- Email: [email protected]
CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles
PTEC-14898; No of Pages 11 Powder Technology xxx (2019) xxx
Contents lists available at ScienceDirect
Powder Technology journal homepage: www.elsevier.com/locate/powtec
CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles Foad Farivar a,⁎, Hu Zhang a,b, Zhao F. Tian c, Anshul Gupte d a
The University of Adelaide, School of Chemical Engineering, Adelaide, SA 5005, Australia Amgen Bioprocessing Centre, Keck Graduate Institute, Claremont, CA 91711, United States of America The University of Adelaide, School of Mechanical Engineering, Adelaide, SA 5005, Australia d Mayne Pharma, Salisbury south, Adelaide, SA 5106, Australia b c
a r t i c l e
i n f o
Article history: Received 9 May 2019 Received in revised form 4 November 2019 Accepted 6 November 2019 Available online xxxx Keywords: Computational fluid dynamics Discrete element method Fluidized bed Multi-sphere method Simulation
a b s t r a c t A hybrid Open-MP (Shared memory) – MPI (Distributed memory) parallel Computational Fluid Dynamics (CFD)Discrete Element Method (DEM) model has been developed for simulating fluidization of cylindrical particles in a fluidized bed. In the DEM model, five, seven, nineteen and sixty overlapping spherical particles in different formation methods including single-row, multi-row and hybrid formations have been used to represent one cylindrical particle. The simulation results of the hybrid representation for mimicking a cylindrical particle are more accurate than other representations; particularly in terms of, particle number density along the bed height and particle orientation distribution in the bed, in comparison with previously reported experimental data. The simulation results show that the preferential orientation of the particles is horizontal in the fluidize bed, but the number of particles with vertical orientation increases when they approach to walls. Only multi-row and hybrid formations could capture this increase near the walls. © 2019 Elsevier B.V. All rights reserved.
1. Introduction Gas-fluidized beds are widely applied in various industries such as chemical, energy, pharmaceutical and environmental industries [1]. Understanding the hydrodynamics and particle movement in the bed is of paramount importance for designing, optimizing and scaling-up in fluidization process. Numerical simulation has played a significant role in revealing the fundamentals in fluidized beds [2,3]. Multi-scale modelling approaches, such as a combination of Discrete Element Method (DEM) and Computational Fluid Dynamics (CFD), has been developed to describe the hydrodynamics of both gas and solid phases in fluidized beds [4,5]. Spherical particles are often adopted for CFD-DEM modelling of fluidized beds for simplicity. However, as most particles in industrial processes are actually non-spherical, recently modelling of nonspherical particles has been receiving an increasing amount of attention [6–8]. Cylindrical particles are common in several industrial processes including polymer suspension, biomass combustion and fluidized beds [9–12]. In pharmaceutical industry, extrusion is a common process in the manufacture of various drug delivery systems such as granules, pellets and tablets which results in cylindrical shape products [13]. Therefore, understanding the aerodynamic behaviour of cylindrical particles is of great importance. Representation of non-spherical particles and ⁎ Corresponding author. E-mail address: [email protected] (F. Farivar).
detection of particle-particle contact are the main challenges of using DEM for non-spherical particles [14]. There are three main approaches for representing a non-spherical particle in DEM, namely multisphere, polyhedral and super quadratic representations [14]. Super quadratic and polyhedral shape representation methods have been used for simulation of various non-spherical particles in fluidized beds [8,15,16]. However, these approaches demand complicated and computationally expensive contact detection algorithms [14]. On the other hand, in the multi-sphere method, a non-spherical body is represented by a certain number of overlapping spheres. This method is easy to be implemented and can be readily extended from the DEM code for spherical particles [17]. Since the multi-sphere approach employs a fast and well-validated sphere-sphere contact detection algorithm, it is the most efficient and versatile in DEM simulation of non-spherical particles [14]. In this representation method, an increase in the number of spheres may lead to more accurate representation of non-spherical surfaces, while this also increases the computational costs [17]. Therefore, one challenge for this multi-sphere approach is to determine the optimal number of spheres [14,17], especially the topology assembled from spheres to capture particle behaviour in the multi-phase system. There are many studies that emphasise the significant effect of the particle shape on fluidized bed dynamics [6,15,16,18,19], but very few on the effect of particle representation [20]. Zhong et al. [21] pioneered in using a multi-sphere approach to investigate particle behaviour of cylindrical particles in fluidized beds.
https://doi.org/10.1016/j.powtec.2019.11.016 0032-5910/© 2019 Elsevier B.V. All rights reserved.
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
2
F. Farivar et al. / Powder Technology xxx (2019) xxx
In their study, 5 overlapping spheres with a diameter of 2.6 mm are used to represent a cylindrical particle with the same diameter and a length of 6 mm. Reasonable accuracy of predictions from this approach has been achieved when comparing their simulation results with experimental data. Ren et al. [22,23] continued to examine the dynamic behaviour of corn-shaped and cylindroid particles in a spouted bed using this multi-sphere representation approach. They used 4 overlapping spheres to represent a 10 mm long cylinder with a diameter of 4 mm. Their simulation results captured the spouting action of non-spherical particles with great accuracy. Recently, Oschmann et al. [20] studied various elongated non-spherical particles in a fluidized bed with a focus on orientation and mixing of cylindrical and cuboidal particle. They compared the results from polyhedral and multi-spherical shape representation and concluded that the accuracy of the particle shape representation has a significant influence on the simulation results [20]. In their study 36 overlapping spheres is used to represent a cylinder of 6 mm length and 6 mm diameter. Despite the popularity of multi-sphere approach for simulation of non-spherical particles in fluidized beds [14,20,22], the effect of number and formation of spheres on the accuracy of simulation results is so far not fully understood. In this study, a detailed investigation on the effect of cylindrical particle representation by the multi-sphere method on particle orientation and bed hydrodynamics in a fluidized bed is presented. The soft-sphere DEM method is used for simulation of particles movement and the MFIX CFD solver for fluid flow. The simulation results, including pressure drop, particle distribution, and particle orientation distribution, are compared to experimental data of Vollmari et al. [16].
!b !s vice versa (M ¼ AM ) [25]: 2
q20 þ q21 −q22 −q23 6 A¼6 4 2ðq1 q2 −q0 q3 Þ
2ðq1 q2 þ q0 q3 Þ q20 −q21 þ q22 −q23
2ðq1 q3 þ q0 q2 Þ
2ðq2 q3 −q0 q1 Þ
2ðq1 q3 −q0 q2 Þ
3
7 2ðq2 q3 þ q0 q1 Þ 7 5 q20 −q21 −q22
þ
ð4Þ
q23
where q0, q1, q2, q3 are components of a unit quaternion which can be represented in terms of Euler angles as [25]: θ ϕþψ q0 ¼ cos þ cos 2 2
ð5Þ
θ ϕ−ψ q1 ¼ sin þ cos 2 2
ð6Þ
θ ϕ−ψ q2 ¼ sin þ sin 2 2
ð7Þ
θ ϕþψ q3 ¼ cos þ sin 2 2
ð8Þ
2.1.1. Contact force To calculate particle-particle and particle-wall contact forces, the linear spring-dashpot model proposed by Cundall and Strack [26] is used. In this model, the normal component of the contact force from the linear spring- dashpot model is calculated as: ! ! ! F ij;n ¼ −kn δn n ij −ηn u ij;n
ð9Þ
2. Mathematical equations 2.1. Solid phase equations In this work a three dimensional DEM code is developed to model the motion of each particle in the system, based on the Newton's second law [24]: ! d u i !c !pf !g mi ¼ Fi þ Fi þ Fi dt ! ! d I i :ω i dt
! ¼ Ti
! ! δn ¼ Ri þ Ri − j r i −r j j ð1Þ
ð2Þ
! ! where mi, u i and ω i are mass, translational and angular velocities of par!c ticle i, respectively. F i is the summation of all contact forces acting on !g !pf particle i. F i and F i are particle-fluid interaction and gravity forces, re! ! spectively. I i is the particle's moment of inertia, and T i is the summation of all torques acting on particle i. Since the moment of inertia changes with the orientation of nonspherical particles, the equations of rotational motion (Eq. (2)) are solved in a body-fixed coordinate system. The axes in this system are the principal axes of inertia, therefore, Eq. (2) in the body-fixed frame can be expressed as [17]:
I0i
!b d ωi !b !b !b þ ω i I 0i ω i ¼ T i dt
where kn and ηn are the normal spring stiffness and the normal damping ! coefficient respectively. n ij is the unit vector from the centre of particle i ! to the centre of particle j, and u ij;n is the normal relative velocity of the contacting particles. The normal overlap δn of spherical particles can be calculated as:
ð3Þ
!b !b where Ii′, ω i and T i are principal moment of inertia, angular velocity in the body-fixed frame and the summation of all torques acting on the particle i in the body- fixed frame, respectively. The rotation matrix A is used to transfer a vector from body-fixed to space fixed frame and
ð10Þ
! ! where r i , r j are position vectors and Ri, Rj are radii of particle i and j, respectively. Normal spring stiffness and the damping coefficient are contact parameters. For a contact between particle i and particle j, these parameters are given as [27]: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RR i j Ri þ R j
kn ¼
4 E 3
E ¼
1−ν 2i 1−ν j þ Ei Ej
ð11Þ
2
2 lnen ηn ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 π2 þ ln en
ð12Þ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mi m j kn mi þ m j
ð13Þ
where E, ν and en are Young's module, the Poisson ratio and coefficient of normal restitution, respectively. Similar to the normal contact force, the tangential force is calculated as [28]: 8 ! ! ! < −μ ! F ij;n t ij if : F ij;t Nμ F ij;n ! F ij;t ¼ ! : −k δ t −η ! t t ij t u ij;t else
ð14Þ
! ! In this equation, when j F ij;t jNμj F ij;n j, sliding occurs and the tangential force is obtained from the Coulomb friction law [28]. In Eq. (14), kt ;
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
F. Farivar et al. / Powder Technology xxx (2019) xxx
! ηt ; u ij;t ; μ and δt are tangential spring stiffness, damping coefficient, tangential relative velocity, friction coefficient and tangential displacement, respectively. The tangential spring stiffness and damping coefficient in Eq. (14) for any contact between particle i and j is given as [27]: kt ¼ 8 G
G ¼
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ri R j δn Ri þ R j
ð15Þ
2−νi 2−ν j þ Gi Gj
2 ln et ηt ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 π 2 þ ln et
ð16Þ
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mi m j kt mi þ m j
ð17Þ
where G and et are the shear module and coefficient of tangential restitution. ! ! The tangential unit vector ( t ij) is directed along u ij;t , which is calcu! lated using the tangential displacement, δt, in case u ij;n is equal to zero [29]. 8 ! u ij;t > > ! ; > > > u ij;t > < ! t ij ¼ δt > ; > > > > jδt j > : 0;
! u ij;t ≠0 ð18Þ
otherwise
2.1.2. Particle- fluid interaction force
!D The main particle-fluid interaction forces are the drag force F i , lift !∇P !L force F i and pressure gradient force F . !pf !D !L !∇P Fi ¼ Fi þ Fi þ F
ð19Þ
Accurate calculation of these forces, especially the drag force is vital for solid – fluid system simulations [30–32]. The equation of drag force acting on each particle is given as [15]: ! ! !D 1 ! ! F i ¼ ρ f C D v − u i A⊥ v − u i ϵ2−X f 2
ð20Þ
! ! where v and u i are fluid and particle velocities, respectively. CD is the drag coefficient, ρf and A⊥ are fluid density and cross-sectional area of the particle normal to the fluid flow, respectively. ϵf is the fluid void fraction and X is given by [33]: " X ¼ 3:7−0:65 exp −
1:5− logRep 2
2 # ð21Þ
where Rep is the particle's Reynolds number. A⊥ in this work is calculated as the cross-sectional area of a cylinder with diameter of D and length of L normal to the fluid velocity (A⊥ = D L sin α). Because of general applicability and high accuracy of the Hӧlzer/ Sommerfeld model for calculation of the drag coefficient for nonspherical particles [9,16], this model is used in this study [34]: CD ¼
considered particle. The crosswise sphericity is the ratio between the cross sectional area of a spherical particle with the equivalent volume and the projected cross sectional area of the considered particle perpendicular to the flow [34]. Lift force is included because of asymmetric pressure distribution on the surface of non-spherical particles [31]. The lift force acts in the direction perpendicular to the drag force and it can be calculated as [35]: !0 ! ! h i z : v−ui !L 1 !0 ! ! ! ! F i ¼ ρ f C L A‖ ! ! z v−ui v−ui v − u i 2
0:2 1 8 1 16 1 3 1 pffiffiffiffi þ pffiffiffiffiffiffiffiffi 3=4 þ 0:42 100:4ð−logΦÞ pffiffiffiffiffiffiffi þ Rep Φ⊥ Rep Φ Φ⊥ Rep Φ
ð22Þ
where Φ and Φ⊥are regular and crosswise sphericities, respectively. The regular sphericity is the ratio between the surface area of a spherical particle with the equivalent volume and the surface area of the
ð23Þ
where A‖ is the particle area, normal to the direction of the lift force and !0 z is the particle's major axis direction. CL is the lift coefficient and it is determined as [31]: CL ¼ sin2 θ cosθ CD
ð24Þ
where, θ is the incident angle between the major axis of the particle and its relative velocity. Finally, the pressure gradient force can be calculated as [17]: !∇P F i ¼ V i ∇P
jδt j≠0
3
ð25Þ
where Vi is the volume of particle i. ∇P is the pressure gradient which is obtained from CFD solver at the particle's location. For non- spherical particles, the pressure centre in most cases does not coincide with the centre of gravity. The distance between the centre of pressure and the centre of gravity for cylindrical particles (Xcp) depends on particle length and orientation, and is obtained from [35]: X cp ¼ 0:25 L 1−e3ð1−L=dÞ cos3 θ
ð26Þ
where, L and θ are the particle's length and incident angle, respectively. As a result, the drag and lift forces induce a rotational movement !L !D due to the additional torques (X cp F and X cp F ) in Eq. (2), on the particles [14]. 2.2. Fluid flow equations The fluid flow is considered as laminar flow and solved using the volume averaged continuity and Navier–Stokes equations as [2]: ∂ ϵf ρf ∂t
! þ ∇:ϵ f ρ f v ¼ 0
! ∂ ϵf ρf v ∂t
ð27Þ
! !! ! ! þ ∇:ϵ f ρ f v v ¼ −ϵ f ∇P− ∇:ϵ f τ f þ ϵ f ρ f g − S p ð28Þ
! The fluid viscus stress tensor, τ f , is given by [4]: 2 ! ! T τ f ¼ − λ f − μ f Þ ∇:ug I−μ f ∇ v þ ð∇! vÞ 3
ð29Þ
where μf is the gas viscosity and I is the unit tensor. λf is the bulk viscosity which is neglected for gases [3]. ! The source term ( S p ) in the momentum equation which is due to momentum exchange with particles in the cell, can be calculated as [17]: kp ! !pf 1 X Sp ¼ F V cell i¼1 i
ð30Þ
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
4
F. Farivar et al. / Powder Technology xxx (2019) xxx
where kp is the number of particles residing in each CFD cell with a volume of Vcell. ϵf in each CFD cell represents the volume fraction of the gas phase and can be estimated as [36]: kp 1 X ϵ f ¼ 1− ϕV V cell i¼1 i i
ð31Þ
where ϕi ∈ [0, 1] is the volume fraction of particle i inside the CFD cell. In this work, a general Lagrangian variable Γp is mapped to an Eulerian cell Γcell and vice versa, by the weighted averaging method as [37]: Γp ¼
kc X
p
f j Γj
ð32Þ
j¼1
Γcell ¼
kp X
cell
f i Γi
ð33Þ
i¼1
where fpj is the volume fraction of the particle residing in cell j and kc is the number of overlapping cells with the particle. fcell i is the volume fraction of the cell occupied by particle i and kp is the number of particles residing in the cell.
Fig. 1. a) Multi-sphere representation of a cylindrical particle using, 5, 7, 60 and 19 spheres. b) The fluidized bed dimensions and the wall zone in the system.
3. Simulation setup and parameters The parameters and dimensions used in this simulation are presented in Table 1. The dimensions of fluid bed and particles are the same as those from the experimental work of Vollmari et al. [16]. As shown in Fig. 1a, the cylindrical particles are represented using a multi-sphere method with 5, 7, 19 and 60 overlapping spheres and named as Cyl-5, Cyl-7, Cyl-19 and Cyl-60, respectively. It is worth noting that all of the multi-sphere representations have the same flow conditions and properties in the simulations, only the number of spheres to represent each non-spherical particle is different. In the simulations, 4185 cylindrical particles are used. This number is calculated from the particle density and a mass of 0.53 kg used in the experiment. The particles are randomly packed and 4 s are allowed to settle down the particles in the bed. Fig. 1b shows the packed particles in the fluidized bed. As shown in this figure, the initial height of the bed in the simulation is roughly 0.1 m. The inlet superficial gas velocity is increased from 0.4 m/s up to 2.4 m/s in an increment of 0.2 m/s and 2 s of simulation is
Table 1 The parameters and dimensions used in the simulation. Model parameters Bed dimensions CFD cell average volume Number of CFD cells Initial pressure Coefficient of restitution (en) Friction coefficient (μn) Rolling friction coefficient (μr) CFD time step DEM time step Particle's properties Cylinder size
0.11 × 0.11 × 0.41 (m) 2.74 × 10−6 (m3) 1920 101,325 (Pa) 0.5 0.4 0.1 1 × 10−4 (s) 1 × 10−5 (s)
Poisson ratio (ν) Density Total number of particles
D = 3.9 (mm) L = 14 (mm) 0.32 764.4 (kg/m3) 4185
Fluid properties Viscosity (μf) Temperature Ambient pressure
1.82 × 10−5 (kg/m.s) 291.15 (K) 101,325 (Pa)
applied for each velocity. The same process is repeated for all of the multi-sphere representations in Fig. 1a. 3.1. Multi-sphere particle representation Representing non-spherical particles in DEM can be done using different methods. Zhong et al. [14] have classified these methods into two groups; composite-particle and single-particle methods. In the composite-particle approach complex shapes are created by combining simple geometries. The most important example is multi- sphere or glued spheres method [14]. This method was first proposed by Gallas and Sokolowski [38] to study the behaviour of granular material with non-spherical particles. In this study, multisphere method is used to represent the cylindrical particles in a gas fluidized bed. There are many possibilities to combine spherical particles with different sizes to form a cylinder-like particle, and there is no systematic method to optimize the combination. We herein aim to provide a general guideline to represent a cylindrical particle using this multi-sphere method for DEM studies. The simplest way to form a cylinder using overlapping spheres is the single-row formation. In this formation a specific number of overlapping spheres with the same diameter as the cylinder is located on the axial direction of the cylinder (along the length). Cyl-5 and Cyl-7 in Table 2 are formed using this combination. In this figure, NL is the number of spheres along the length (axial direction) and NR is the number of spheres along radial direction of the cylinder. In the single-row formation, increasing the number of spheres results in a smoother surface along the length and more accurate representation of the cylinder. To quantify this accuracy, the best parameter is the distance between centres of adjacent spheres (Dl). A higher Dl means a less accurate lengthwise representation and vice versa. Although the single row formation with a small number of spheres decreases the computational cost, this formation poorly represent the flat ends of a cylinder. One possible way to improve this is to use a multi-row formation in which spheres with a smaller size are located in multiple rows along the length of the cylinder. Cyl-60 in Table 2 is an example of this formation. In this formation, to measure the accuracy of radial representation at both ends of the cylinder, similar to the lengthwise representation Dl, Dr is defined as the longest distance between adjacent spheres forming the ends of a cylinder. For single-row
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
F. Farivar et al. / Powder Technology xxx (2019) xxx
5
Table 2 The parameters and dimensions used in the simulation. Cyl-5
Cyl-7
Cyl-60
Cy-19
5×1 (NL × NR) 2.52 × 10−3 3.90 × 10−3
7×1 (NL × NR) 1.68 × 10−3 3.90 × 10−3
10 × 6 (NL × NR) 1.33 × 10−3 1.90 × 10−3
7 + (2 × 6) (NL + 2 × NR) 1.68 × 10−3 0.95 × 10−3
3D representation
2D representation Number of spheres Dl (m) Dr(m)
formations, since there is no sphere to represent the flat ends of a cylinder Dr is the diameter of the cylinder, meaning the maximum error along the radial direction. Apart from multi-row formation, another strategy to improve the representation of the flat ends is to adopt a hybrid formation. In this hybrid formation the main body of the cylinder is represented by single-row formation and the diameter of each sphere is same as the cylinder diameter, but for both ends the multi-row formation is used and the diameter of spheres is half of the cylinder diameter. Cyl-19 in Table 2 is an example of the hybrid formation. Fig. 2 schematically represents Dr and Dl for single- row, multi- row and hybrid formations. The values of Dr and Dl for Cyl-5, Cyl7, Cyl-60 and Cyl-19 are presented in Table 2. Cyl-5 and Cyl-7 have the highest Dr with 0.0039 (m) which means they have the lowest accuracy to represent the flat ends of the cylinder. This number in Cyl-60 and Cyl-19 representations decreased to 0.0019 (m) and 0.00095(m), respectively. Dl is 0.00168 (m) for Cyl-7 and Cyl-19, and 0.00133 (m) for Cyl-60. In Section 5, the simulation results of particle dynamics in a fluidized bed using each representation are compared with experiment data to reveal the effect of shape representation on the simulation results in a multi-particle system.
solvers are run on two separate nodes in parallel and Message Passing Interface (MPI) is used to exchange data between two solvers. Each node has access to 16 shared memory processors and the OpenMP parallelization method is used separately in both CFD and DEM code. After initialization of CFD and DEM, CFD data such as pressure and gas velocity are passed to DEM. In the DEM code after defining the initial particle positions and velocities, a combination of No Binary Search (NBS) and Munjiza algorithm [17] is used to detect contacts for particle-particle and particle-wall. The contact forces are calculated using the DEM soft-sphere model based on overlapped contacts. Then the particle-fluid interaction forces acting on each particle are calculated using the CFD data. Finally, new positions and velocities of particles are obtained by solving the differential equations based on the Newton's second law (Eqs. (1) and (2)). This process repeats 10 times to match with the CFD time step. DEM data such as particles position and particle-fluid interaction forces are sent to the CFD solver, where the volume-averaging method is applied to map the DEM data into CFD cells. After calculation of bed void fraction, the fluid flow equations (Eqs. (27) and (28)) are solved to update fluid flow velocities in each cell. This process is completed at the final time step of CFD.
4. Coupled CFD-DEM algorithm
5. Results and discussions
In this work, an open source MFIX CFD solver [39] is coupled with a three dimensional in-house DEM code. The DEM code is written in FORTRAN 95 programming language and it is compatible with the CFD solver. They are compiled by gfortran compiler on a LINUX system. Fig. 3 shows the general algorithm of the CFD-DEM model. The developed CFD-DEM coupling is parallelized using a combination model for shared memory and distributed memory. CFD and DEM
5.1. Initial fluidization process The snapshots of initial fluidization process at a superficial velocity of 2.4 m/s for Cyl-7 are shown in Fig. 4a. At the initial stage (t = 0.0 s), particles are compacted in the bottom of the bed mostly in horizontal orientation. At t = 0.1 s, gas permeates through the gap between particles to create more volume voidage in the bed. The incoming gas
Fig. 2. Schematic representation of Dr and Dl for single-row, multi-row and hybrid formations.
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
6
F. Farivar et al. / Powder Technology xxx (2019) xxx
Fig. 3. The CFD-DEM model algorithm.
lifts the particles out of the bottom at t = 0.15 s and the bed height changes slightly in comparison with t = 0.1 s. After 0.2 s, the top part of the bed is fluidized, while a dense structure of particle packing remains at the bottom of the fluidized bed. Spherical particles packed in the bed are also simulated for fluidization at the same superficial velocity. Gas penetrates through gaps between spherical particles and lifts the particles at t = 0.1 s, resulting in an increase in the bed voidage in the middle of the bed. At t = 0.2 s, air partially fluidizes particles in the middle of bed and pushes the particles to the sides of the bed, and then the particles return to the bottom part of the fluidized bed, re-forming a packed structure at t = 0.3 s. The difference in the fluidization behaviour between non-spherical and spherical particles may be due to particle-particle and particlewall interactions. To further understand the fluidization process of non-spherical particles, pressure drop, particle density and particle orientation distribution are presented below by averaging the simulation results from 100 time steps (2 s of simulation with an interval of 0.02 s) for each superficial velocity.
drops for Cyl-5 and Cyl-7 are nearly identical at most of gas velocities, while the maximum difference of 26 Pa is found at the velocity of 2.4 m/s. Simulations for Cyl-5 and Cyl-7 slightly under-predict the pressure drop after the fluidization occurs at around 0.8 m/s. The simulation results for Cyl-60 and Cyl-19 have the better match with the experiment data. However, between 0.8 m/s and 0.9 m/s, where the fluidization starts, the simulation cannot capture the peak of the pressure drop. This peak is the result of interlocking effect between cylindrical particles in the bed. Since the static friction is not modelled in the DEM and only simplified constant friction factor is included in the calculations, this model cannot accurately capture this peak of pressure drop at the onset of fluidisation. Fig. 5b shows that the simulation results for the pressure drop of spherical particles with the diameter of 7 mm at different gas superficial velocities is in good agreement to the experiment. In contrast to cylindrical particle, the pressure drop for spherical particles does not drop after the fluidization. This is because spherical particles have smooth surface and no interlocking happens between them.
5.2. Pressure drop
5.3. Particle density distribution
The average pressure drop in the bed for Cyl-5, Cyl-7, Cyl-60 and Cyl19 with different gas superficial velocities is displayed in Fig. 5a. The results for each gas velocity are presented as the average data for 2 s simulation. The simulation results are in good agreement with the experiment in which the pressure at the top of the bed (height = 0.31 m) and below the flow distribution plate is measured by a pressure transducer [16]. As the superficial velocity increases from 0.4 m/s to 0.8 m/s, the pressure drop in the bed increases from around 150 Pa to 370 Pa. At 0.8 m/s, the pressure drop reaches a plateau. This is a typical pressure drop profile with the gas superficial velocity in the fluidized bed and the onset of fluidization occurs at a superficial velocity of 0.8 m/s. For different multi-spheres to represent one cylindrical particle, pressure
In Fig. 6, the average number of particles at the cross-section of different bed heights is used to present the particle spatial distribution in the bed. In this figure, the values of the predicted number density at gas superficial velocities of 1.6 m/s, 2.0 m/s and 2.4 m/s for Cyl-5, Cyl7, Cyl-60 and Cyl-19 are compared to the experiment data. Generally, the predicted data for Cyl-5, Cyl-7, Cyl-60 and Cyl-19 have a similar trend as the experimental results: a relatively higher number of particles in the bottom than on the top of the bed, suggesting more particles are accumulated in the bottom and less particles fluidize at the top. This is also verified from the snapshots of particle fluidization in Fig. 4. As the superficial velocity increases, particles are more dispersed, which is indicated by the average number density of 0.07 for 1.6 m/s down to 0.05 for 2.4 m/s at lower part of the bed in Fig. 6.
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
F. Farivar et al. / Powder Technology xxx (2019) xxx
7
Fig. 4. Paraview snapshots of the fluidized bed at a gas superficial velocity of 2.4 m/s for a) Cyl-7 for non-spherical particles with d = 3.9 mm and L = 14 mm, and b) spherical particle with d = 7 mm.
At a gas superficial velocity of 1.6 m/s (Fig. 6a), simulation results for Cyl-5 have the largest average deviation of 14.4e-3 from the experimental data. This number for Cyl-7 and Cyl-60 reduces to 7.5e-3 and 5.9e-3, respectively. Simulation results for Cyl-19 have the lowest average deviation of 5.8e-3. Fig. 6(b, c) show that all of the cylinder representations except Cyl-5, have a good match with the experimental observations at velocities of 2.0 m/s and 2.4 m/s. Simulation data for Cyl-60 and Cyl-19 have comparable accuracy with an average deviation of 5.26e-3 and 5.19e-3 at 2.0 m/s, respectively. The simulation results for the particles density distribution show that increasing the number of spheres in the single-row formation to decrease the lengthwise representation error, improves the accuracy of the model substantially. On the other hand, using both multi-row and hybrid formations to further improve the multi-sphere representation, reduce the simulation errors and improves the accuracy of the model roughly the same amount. However, in some cases such as at Vg = 2 m/s, this improvement can be marginal. 5.4. Particle orientation distribution Unlike spherical particles, particle orientation plays a significant role on cylindrical particle behaviour in the fluidized bed. In this study, the angle of 0.0° refers to the particle in the horizontal orientation aligning with the longest dimension, while the angle of 90.0° for the vertical orientation. We attempt to compare the particle orientation in the fluidized bed with experimental data. Since the experimental imaging may only capture the particles close to the experimental walls, two
simulation results are obtained: one from averaged particles in the domain (Cyl-5, Cyl-7, Cyl-60 and Cyl-19), and the other from averaged data only in the vicinity of walls (thickness of 0.01 m) with a defined wall zone shown in Fig. 1 (Cyl-5 wall zone, Cyl-7 wall zone, Cyl-60 wall zone and Cyl-19 wall zone). The number density vs the oriental angle at the superficial velocity of 1.6 m/s, 2.0 m/s and 2.4 m/s for Cyl-5, Cyl-7, Cyl-60 and Cyl-19 is presented in Fig. 7. The experimental data are also included in Fig. 7 for comparison. Simulation results for average orientation of the particles in the fluidized bed are presented in Fig. 7 for all particles in the whole bed (red lines with circular marker) as well as particles in the vicinity of the walls (green lines with triangular marker). At a superficial velocity of 1.6 m/s (Fig. 7a, d, g, j), experimental results [16] show that horizontal (between 0° and 10°) and vertical (between 80° and 90°) orientations are dominant for the particles in the bed. Particles in the horizontal orientation are relatively stable while they are likely to orient vertically (parallel to the walls) when they approach to walls. Similar wall effect on particles orientation is also reported in a previous study [8]. The number of particles decreases as the angle increases from horizontal to vertical orientation except at the angle of 90o. The simulation data for the average number of particles in the whole domain or in the wall zone follow a very similar trend as the experimental observations. However, predication for the number of particles at the full vertical orientation (90o) fails for Cyl-5 and Cyl7, while the improvement of the predications is shown in Fig. 7g for Cyl-60 and Fig. 7j for Cyl-19 in the wall zone. The poor representation
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
8
F. Farivar et al. / Powder Technology xxx (2019) xxx
Fig. 5. Pressure drop of the fluidized bed obtained from experiment [16] and simulation results for, (a) Cyl-5, Cyl-7, Cyl-60 and Cyl-19 and (b) spherical particles with d = 7 mm.
of both ends in single-row formations (Cyl-5 and Cyl-7) may lead to unstable vertical orientation (standing position), and the round ends are not able to account for accurate particle-wall interactions. The simulation data averaged in the whole domain indicate that the vertical orientation is the least preferred for particles, opposite to the experimental results. This may be due to the difference in counting particles by the experimental imaging method. The two dimensional nature of experimental image analysis suggests that particles close to the walls are more likely to be detected. An increase in the superficial velocity does not affect the trend of the number density versus the oriental angle for Cyl-5, Cyl-7, Cyl-60 and Cyl-19, while the number density with vertical orientation slightly increases by increasing the superficial velocity, which is supported from both experimental and simulation data. In all multi-sphere representations, simulation results of particles' orientation distribution for hybrid formation (Cyl-19) have the highest accuracy especially at a gas superficial velocity of 2.4 m/s. The simulation results suggest that an increasing number of spheres in the multisphere method using the same formation method can improve the accuracy of simulation results. However, the best representation is not always the one with higher number of spheres. After averaging the particle orientation of all particles on a cross section at a specific bed height, the averaged orientation at different bed heights is shown in Fig. 8. Particles preferably stay in the horizontal orientation at the bottom of the bed. Slightly above the bottom of the bed, the average angle fluctuates around 35° for the velocity of 1.6 m/s and 2.0 m/s, while for 1.0 m/s it remains stable around 20°. This behaviour can be attributed to the effect of particle-particle and particle-wall interactions on the average oriental angle along the vertical direction [40]. By increasing the gas superficial velocity, the number of particle- particle and particle- wall interaction of fluidizing particles increases, which leads to an increase in the particles' average orientation [40].
Fig. 6. Average height of the particles at gas superficial velocities of a) 1.6 (m/s), b) 2.0 (m/ s) and c) 2.4 (m/s).
5.5. The effect of torque induced by drag and lift forces Cylindrical particles experience rotation due to particle-particle and particle-wall interactions as well as drag and lift forces exerted from fluid on particles. The most stable orientation for both stationary and fluidizing particles with a cylindrical shape is horizontal (0°), and the fluidizing particles would have no preferred orientation without any influence of fluid flow on the particle orientation. In order to reveal the effect of torque induced by drag and lift forces, we conducted two separate simulations using Cyl-60: one with the assumption of Xcp = 0 in Eq. 27 and the other including this parameter in every time step (Xcp ≠ 0). The average number of particles at different orientation angles for two scenarios are shown in Fig. 9. At a gas superficial velocity of 1.6 m/s including the torque induced by fluid flow, the number density vs the oriental angle is very similar
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
F. Farivar et al. / Powder Technology xxx (2019) xxx
9
Fig. 7. Particle orientation distribution averaged through all of the particles and particles in the wall zone, for Cyl-5 (a, b, c), Cyl-7 (d, e, f), Cyl-60 (g, h, I) and Cyl-19 (j, k, l) at gas superficial velocities of 1.6 m/s (left column), 2 m/s (middle column) and 2.4 m/s (right column).
to that without the torque, which indicates that the fluid flow torque has no effect on the average orientation. At this velocity, most of the particles retain at the bottom of the bed and they are in the horizontal
orientation. Thereby, fluid flow has a negligible influence on the average orientation of the particles. However, the difference in the particle orientation distribution becomes very significant at high velocities between with and without the torque as shown in Fig. 9b and c. Since particle-fluid interaction forces (drag and lift forces) become more dominant, the torque caused by these forces becomes more significant. Therefore, at higher gas superficial velocities, the torque effect induced by drag and lift forces torque on cylindrical particle orientation should be considered. 6. Simulation run time
Fig. 8. Average orientation along the bed height for gas superficial velocities of 1.0, 1.6 and 2.0 m/s (Cyl-60).
Fig. 10 compares the computing time for 1 s simulation of fluidization process for Cyl_5, Cyl_7, Cyl_60 and Cyl-19. For the Cyl_5 representation, the total number of spheres in the simulation was 20,925 (5 × 4185) and 0.52 h is spent to simulate the process for 1 s. The number of spheres used for Cyl_7 representation was 29,295 and it takes 1.22 h to complete the simulation of the fluidization process. As expected, the simulations using Cyl_60 representation was the most computationally expensive among four representations: about 4.08 h. This number is almost 2 times higher than Cyl_19 run time at 2.6 h, while the accuracy of
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
10
F. Farivar et al. / Powder Technology xxx (2019) xxx
Fig. 9. Comparison of the average particles' orientation distribution for assumptions of Xcp = 0 (drag and lift forces do not cause rotation) and Xcp ≠ 0 (drag and lift forces cause rotation), at gas superficial velocity of, a) 1.6 m/s, b) 2 m/s and c) 2.4 m/s.
Acknowledgements This work was supported by an Australian Government Research Training Program (RTP) Scholarship; and Mayne Pharma International PTY LTD. The authors also acknowledge the supercomputing resources were provided by the Phoenix HPC service at the University of Adelaide. References
Fig. 10. Comparison of the computing time between Cyl_5, Cyl_7, Cyl_60 and Cyl-19 representations for simulation of one second of fluidization process.
the simulation results is almost same regarding particle density distribution and slightly better for Cyl-19 representation in particle orientation distribution in the fluidized bed.
7. Conclusions In this work a CFD-DEM model was developed to reveal the fluidization process of cylinder-like particles in a fluidized bed using different numbers of spheres and formation methods to represent one cylindrical particle by using the multi-sphere method. The simulation results of bed pressure drop using multi-row and hybrid formations (Cyl-60 and Cyl19) have better agreement with experimental data than those from single-row formations (Cyl-5 and Cyl-7). The average axial distribution of the particles is better mimicked by increasing the number of spheres from 5 to 7 in single-row formation. In both multi-row and hybrid formations the simulation results accuracy are slightly better than singlerow formation with 7 spheres. The preferential orientation is horizontal at all gas velocities, but the number of particles with vertical orientation increase when particle approach to walls and this increase is captured with multi-row and hybrid formations. In general, the hybrid formation results has higher accuracy and much lower computational time compared to the multi-row formation in this study. Furthermore, the torque effect induced by drag and lift forces from fluid flow on average particle orientation becomes significant at a high superficial velocity. This work also provides general guidelines for non-spherical particles representation using multi-sphere method and investigates the effect of number of spheres and formations methods on the dynamics of fluidized beds.
[1] C.T. Crowe, Multiphase Flow Handbook, CRC press, 2005. [2] M.A. van der Hoef, M. van Sint Annaland, J.A.M. Kuipers, Computational fluid dynamics for dense gas–solid fluidized beds: a multi-scale modeling strategy, Chem. Eng. Sci. 59 (22−23) (2004) 5157–5165. [3] M.A. van der Hoef, et al., Multiscale modeling of gas-fluidized beds, in: B.M. Guy (Ed.), Advances in Chemical Engineering, Academic Press 2006, pp. 65–149. [4] N.G. Deen, et al., Review of discrete particle modeling of fluidized beds, Chem. Eng. Sci. 62 (1–2) (2007) 28–44. [5] H.P. Zhu, et al., Discrete particle simulation of particulate systems: a review of major applications and findings, Chem. Eng. Sci. 63 (23) (2008) 5728–5770. [6] Z.Y. Zhou, et al., Discrete particle simulation of gas fluidization of ellipsoidal particles, Chem. Eng. Sci. 66 (23) (2011) 6128–6145. [7] J. Cai, Z. Yuan, X. Zhao, Study on two-way coupling of gas–solid two-phase flow of cylindrical particles, Powder Technol. 300 (2016) 136–145. [8] H. Ma, Y. Zhao, CFD-DEM investigation of the fluidization of binary mixtures containing rod-like particles and spherical particles in a fluidized bed, Powder Technol. 336 (2018) 533–545. [9] H. Ma, L. Xu, Y. Zhao, CFD-DEM simulation of fluidization of rod-like particles in a fluidized bed, Powder Technol. 314 (2017) 355–366. [10] M. Mohseni, B. Peters, M. Baniasadi, Conversion analysis of a cylindrical biomass particle with a DEM-CFD coupling approach, Case Stud. Therm. Eng. 10 (2017) 343–356. [11] L. Boer, et al., Experimental study on orientation and de-mixing phenomena of elongated particles in gas-fluidized beds, Powder Technol. 329 (2018) 332–344. [12] X. Chen, W. Zhong, T.J. Heindel, Orientation of cylindrical particles in a fluidized bed based on stereo X-ray particle tracking velocimetry (XPTV), Chem. Eng. Sci. 203 (2019) 104–112. [13] J. Thiry, F. Krier, B. Evrard, A review of pharmaceutical extrusion: critical process parameters and scaling-up, Int. J. Pharm. 479 (1) (2015) 227–240. [14] W. Zhong, et al., DEM/CFD-DEM modelling of non-spherical particulate systems: theoretical developments and applications, Powder Technol. 302 (2016) 108–152. [15] J.E. Hilton, L.R. Mason, P.W. Cleary, Dynamics of gas–solid fluidised beds with nonspherical particle geometry, Chem. Eng. Sci. 65 (5) (2010) 1584–1596. [16] K. Vollmari, R. Jasevičius, H. Kruggel-Emden, Experimental and numerical study of fluidization and pressure drop of spherical and non-spherical particles in a model scale fluidized bed, Powder Technol. 291 (2016) 506–521. [17] H.R. Norouzi, et al., Coupled CFD-DEM Modeling: Formulation, Implementation and Application to Multiphase Flows, John Wiley & Sons, 2016. [18] B. Liu, et al., Fluidization of non-spherical particles: sphericity, Zingg factor and other fluidization parameters, Particuology 6 (2) (2008) 125–129. [19] G. Lu, J. Third, C. Müller, Discrete element models for non-spherical particle systems: from theoretical developments to applications, Chem. Eng. Sci. 127 (2015) 425–465. [20] T. Oschmann, J. Hold, H. Kruggel-Emden, Numerical investigation of mixing and orientation of non-spherical particles in a model type fluidized bed, Powder Technol. 258 (2014) 304–323. [21] W.Q. Zhong, et al., Discrete element method simulation of cylinder-shaped particle flow in a gas-solid fluidized bed, Chem. Eng. Technol. 32 (3) (2009) 386–391. [22] B. Ren, et al., CFD-DEM simulation of spouting of corn-shaped particles, Particuology 10 (5) (2012) 562–572. [23] B. Ren, et al., Numerical simulation of spouting of cylindroid particles in a spouted bed, Can. J. Chem. Eng. 92 (5) (2014) 928–934. [24] H.P. Zhu, et al., Discrete particle simulation of particulate systems: theoretical developments, Chem. Eng. Sci. 62 (13) (2007) 3378–3396. [25] H.-G.M.J. Chen, Understanding the Discrete Element Method: Simulation of NonSpherical Particles for Granular and Multi-body Systems, John Wiley & Sons, 2014 448. [26] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Géotechnique 29 (1) (1979) 47–65.
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
F. Farivar et al. / Powder Technology xxx (2019) xxx [27] L. Fries, et al., DEM–CFD modeling of a fluidized bed spray granulator, Chem. Eng. Sci. 66 (11) (2011) 2340–2355. [28] B.P.B. Hoomans, Granular Dynamics of Gas-Solid Two-Phase Flows, Universiteit Twente, 2000. [29] A. Džiugys, B. Peters, An approach to simulate the motion of spherical and nonspherical fuel particles in combustion chambers, Granul. Matter 3 (4) (2001) 231–266. [30] Y. Chen, J.R. Third, C.R. Müller, A drag force correlation for approximately cubic particles constructed from identical spheres, Chem. Eng. Sci. 123 (2015) 146–154. [31] M. Zastawny, et al., Derivation of drag and lift force and torque coefficients for nonspherical particles in flows, Int. J. Multiphase Flow 39 (2012) 227–239. [32] Y.Q. Feng, A.B. Yu, Assessment of model formulations in the discrete particle simulation of gas−solid flow, Ind. Eng. Chem. Res. 43 (26) (2004) 8378–8390. [33] R. Di Felice, The voidage function for fluid-particle interaction systems, Int. J. Multiphase Flow 20 (1) (1994) 153–159. [34] A. Hölzer, M. Sommerfeld, New simple correlation formula for the drag coefficient of non-spherical particles, Powder Technol. 184 (3) (2008) 361–365.
11
[35] C. Yin, et al., Modelling the motion of cylindrical particles in a nonuniform flow, Chem. Eng. Sci. 58 (15) (2003) 3489–3498. [36] B.P.B. Hoomans, et al., Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidised bed: a hard-sphere approach, Chem. Eng. Sci. 51 (1) (1996) 99–118. [37] M.v.S. Annaland, N.G. Deen, J.A.M. Kuipers, Chapter 23 Multi-level computational fluid dynamics models for the description of particle mixing and granulation in fluidized beds, in: M.J.H.A.D. Salman, J.P.K. Seville (Eds.), Handbook of Powder Technology, Elsevier Science B.V. 2007, pp. 1071–1107. [38] J.A.C. Gallas, S. Sokolowski, Grain non-sphericity effects on the angle of repose of granular material, Int. J. Mod. Phys. B 07 (09n10) (1993) 2037–2046. [39] H. Liu, D.K. Tafti, T. Li, Hybrid parallelism in MFIX CFD-DEM using OpenMP, Powder Technol. 259 (2014) 22–29. [40] F. Farivar, et al., Capturing particle-particle interactions for cylindrical fibrous particles in different flow regimes, Powder Technol. 330 (2018) 418–424.
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
Contents lists available at ScienceDirect
Powder Technology journal homepage: www.elsevier.com/locate/powtec
CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles Foad Farivar a,⁎, Hu Zhang a,b, Zhao F. Tian c, Anshul Gupte d a
The University of Adelaide, School of Chemical Engineering, Adelaide, SA 5005, Australia Amgen Bioprocessing Centre, Keck Graduate Institute, Claremont, CA 91711, United States of America The University of Adelaide, School of Mechanical Engineering, Adelaide, SA 5005, Australia d Mayne Pharma, Salisbury south, Adelaide, SA 5106, Australia b c
a r t i c l e
i n f o
Article history: Received 9 May 2019 Received in revised form 4 November 2019 Accepted 6 November 2019 Available online xxxx Keywords: Computational fluid dynamics Discrete element method Fluidized bed Multi-sphere method Simulation
a b s t r a c t A hybrid Open-MP (Shared memory) – MPI (Distributed memory) parallel Computational Fluid Dynamics (CFD)Discrete Element Method (DEM) model has been developed for simulating fluidization of cylindrical particles in a fluidized bed. In the DEM model, five, seven, nineteen and sixty overlapping spherical particles in different formation methods including single-row, multi-row and hybrid formations have been used to represent one cylindrical particle. The simulation results of the hybrid representation for mimicking a cylindrical particle are more accurate than other representations; particularly in terms of, particle number density along the bed height and particle orientation distribution in the bed, in comparison with previously reported experimental data. The simulation results show that the preferential orientation of the particles is horizontal in the fluidize bed, but the number of particles with vertical orientation increases when they approach to walls. Only multi-row and hybrid formations could capture this increase near the walls. © 2019 Elsevier B.V. All rights reserved.
1. Introduction Gas-fluidized beds are widely applied in various industries such as chemical, energy, pharmaceutical and environmental industries [1]. Understanding the hydrodynamics and particle movement in the bed is of paramount importance for designing, optimizing and scaling-up in fluidization process. Numerical simulation has played a significant role in revealing the fundamentals in fluidized beds [2,3]. Multi-scale modelling approaches, such as a combination of Discrete Element Method (DEM) and Computational Fluid Dynamics (CFD), has been developed to describe the hydrodynamics of both gas and solid phases in fluidized beds [4,5]. Spherical particles are often adopted for CFD-DEM modelling of fluidized beds for simplicity. However, as most particles in industrial processes are actually non-spherical, recently modelling of nonspherical particles has been receiving an increasing amount of attention [6–8]. Cylindrical particles are common in several industrial processes including polymer suspension, biomass combustion and fluidized beds [9–12]. In pharmaceutical industry, extrusion is a common process in the manufacture of various drug delivery systems such as granules, pellets and tablets which results in cylindrical shape products [13]. Therefore, understanding the aerodynamic behaviour of cylindrical particles is of great importance. Representation of non-spherical particles and ⁎ Corresponding author. E-mail address: [email protected] (F. Farivar).
detection of particle-particle contact are the main challenges of using DEM for non-spherical particles [14]. There are three main approaches for representing a non-spherical particle in DEM, namely multisphere, polyhedral and super quadratic representations [14]. Super quadratic and polyhedral shape representation methods have been used for simulation of various non-spherical particles in fluidized beds [8,15,16]. However, these approaches demand complicated and computationally expensive contact detection algorithms [14]. On the other hand, in the multi-sphere method, a non-spherical body is represented by a certain number of overlapping spheres. This method is easy to be implemented and can be readily extended from the DEM code for spherical particles [17]. Since the multi-sphere approach employs a fast and well-validated sphere-sphere contact detection algorithm, it is the most efficient and versatile in DEM simulation of non-spherical particles [14]. In this representation method, an increase in the number of spheres may lead to more accurate representation of non-spherical surfaces, while this also increases the computational costs [17]. Therefore, one challenge for this multi-sphere approach is to determine the optimal number of spheres [14,17], especially the topology assembled from spheres to capture particle behaviour in the multi-phase system. There are many studies that emphasise the significant effect of the particle shape on fluidized bed dynamics [6,15,16,18,19], but very few on the effect of particle representation [20]. Zhong et al. [21] pioneered in using a multi-sphere approach to investigate particle behaviour of cylindrical particles in fluidized beds.
https://doi.org/10.1016/j.powtec.2019.11.016 0032-5910/© 2019 Elsevier B.V. All rights reserved.
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
2
F. Farivar et al. / Powder Technology xxx (2019) xxx
In their study, 5 overlapping spheres with a diameter of 2.6 mm are used to represent a cylindrical particle with the same diameter and a length of 6 mm. Reasonable accuracy of predictions from this approach has been achieved when comparing their simulation results with experimental data. Ren et al. [22,23] continued to examine the dynamic behaviour of corn-shaped and cylindroid particles in a spouted bed using this multi-sphere representation approach. They used 4 overlapping spheres to represent a 10 mm long cylinder with a diameter of 4 mm. Their simulation results captured the spouting action of non-spherical particles with great accuracy. Recently, Oschmann et al. [20] studied various elongated non-spherical particles in a fluidized bed with a focus on orientation and mixing of cylindrical and cuboidal particle. They compared the results from polyhedral and multi-spherical shape representation and concluded that the accuracy of the particle shape representation has a significant influence on the simulation results [20]. In their study 36 overlapping spheres is used to represent a cylinder of 6 mm length and 6 mm diameter. Despite the popularity of multi-sphere approach for simulation of non-spherical particles in fluidized beds [14,20,22], the effect of number and formation of spheres on the accuracy of simulation results is so far not fully understood. In this study, a detailed investigation on the effect of cylindrical particle representation by the multi-sphere method on particle orientation and bed hydrodynamics in a fluidized bed is presented. The soft-sphere DEM method is used for simulation of particles movement and the MFIX CFD solver for fluid flow. The simulation results, including pressure drop, particle distribution, and particle orientation distribution, are compared to experimental data of Vollmari et al. [16].
!b !s vice versa (M ¼ AM ) [25]: 2
q20 þ q21 −q22 −q23 6 A¼6 4 2ðq1 q2 −q0 q3 Þ
2ðq1 q2 þ q0 q3 Þ q20 −q21 þ q22 −q23
2ðq1 q3 þ q0 q2 Þ
2ðq2 q3 −q0 q1 Þ
2ðq1 q3 −q0 q2 Þ
3
7 2ðq2 q3 þ q0 q1 Þ 7 5 q20 −q21 −q22
þ
ð4Þ
q23
where q0, q1, q2, q3 are components of a unit quaternion which can be represented in terms of Euler angles as [25]: θ ϕþψ q0 ¼ cos þ cos 2 2
ð5Þ
θ ϕ−ψ q1 ¼ sin þ cos 2 2
ð6Þ
θ ϕ−ψ q2 ¼ sin þ sin 2 2
ð7Þ
θ ϕþψ q3 ¼ cos þ sin 2 2
ð8Þ
2.1.1. Contact force To calculate particle-particle and particle-wall contact forces, the linear spring-dashpot model proposed by Cundall and Strack [26] is used. In this model, the normal component of the contact force from the linear spring- dashpot model is calculated as: ! ! ! F ij;n ¼ −kn δn n ij −ηn u ij;n
ð9Þ
2. Mathematical equations 2.1. Solid phase equations In this work a three dimensional DEM code is developed to model the motion of each particle in the system, based on the Newton's second law [24]: ! d u i !c !pf !g mi ¼ Fi þ Fi þ Fi dt ! ! d I i :ω i dt
! ¼ Ti
! ! δn ¼ Ri þ Ri − j r i −r j j ð1Þ
ð2Þ
! ! where mi, u i and ω i are mass, translational and angular velocities of par!c ticle i, respectively. F i is the summation of all contact forces acting on !g !pf particle i. F i and F i are particle-fluid interaction and gravity forces, re! ! spectively. I i is the particle's moment of inertia, and T i is the summation of all torques acting on particle i. Since the moment of inertia changes with the orientation of nonspherical particles, the equations of rotational motion (Eq. (2)) are solved in a body-fixed coordinate system. The axes in this system are the principal axes of inertia, therefore, Eq. (2) in the body-fixed frame can be expressed as [17]:
I0i
!b d ωi !b !b !b þ ω i I 0i ω i ¼ T i dt
where kn and ηn are the normal spring stiffness and the normal damping ! coefficient respectively. n ij is the unit vector from the centre of particle i ! to the centre of particle j, and u ij;n is the normal relative velocity of the contacting particles. The normal overlap δn of spherical particles can be calculated as:
ð3Þ
!b !b where Ii′, ω i and T i are principal moment of inertia, angular velocity in the body-fixed frame and the summation of all torques acting on the particle i in the body- fixed frame, respectively. The rotation matrix A is used to transfer a vector from body-fixed to space fixed frame and
ð10Þ
! ! where r i , r j are position vectors and Ri, Rj are radii of particle i and j, respectively. Normal spring stiffness and the damping coefficient are contact parameters. For a contact between particle i and particle j, these parameters are given as [27]: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RR i j Ri þ R j
kn ¼
4 E 3
E ¼
1−ν 2i 1−ν j þ Ei Ej
ð11Þ
2
2 lnen ηn ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 π2 þ ln en
ð12Þ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mi m j kn mi þ m j
ð13Þ
where E, ν and en are Young's module, the Poisson ratio and coefficient of normal restitution, respectively. Similar to the normal contact force, the tangential force is calculated as [28]: 8 ! ! ! < −μ ! F ij;n t ij if : F ij;t Nμ F ij;n ! F ij;t ¼ ! : −k δ t −η ! t t ij t u ij;t else
ð14Þ
! ! In this equation, when j F ij;t jNμj F ij;n j, sliding occurs and the tangential force is obtained from the Coulomb friction law [28]. In Eq. (14), kt ;
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
F. Farivar et al. / Powder Technology xxx (2019) xxx
! ηt ; u ij;t ; μ and δt are tangential spring stiffness, damping coefficient, tangential relative velocity, friction coefficient and tangential displacement, respectively. The tangential spring stiffness and damping coefficient in Eq. (14) for any contact between particle i and j is given as [27]: kt ¼ 8 G
G ¼
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ri R j δn Ri þ R j
ð15Þ
2−νi 2−ν j þ Gi Gj
2 ln et ηt ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 π 2 þ ln et
ð16Þ
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mi m j kt mi þ m j
ð17Þ
where G and et are the shear module and coefficient of tangential restitution. ! ! The tangential unit vector ( t ij) is directed along u ij;t , which is calcu! lated using the tangential displacement, δt, in case u ij;n is equal to zero [29]. 8 ! u ij;t > > ! ; > > > u ij;t > < ! t ij ¼ δt > ; > > > > jδt j > : 0;
! u ij;t ≠0 ð18Þ
otherwise
2.1.2. Particle- fluid interaction force
!D The main particle-fluid interaction forces are the drag force F i , lift !∇P !L force F i and pressure gradient force F . !pf !D !L !∇P Fi ¼ Fi þ Fi þ F
ð19Þ
Accurate calculation of these forces, especially the drag force is vital for solid – fluid system simulations [30–32]. The equation of drag force acting on each particle is given as [15]: ! ! !D 1 ! ! F i ¼ ρ f C D v − u i A⊥ v − u i ϵ2−X f 2
ð20Þ
! ! where v and u i are fluid and particle velocities, respectively. CD is the drag coefficient, ρf and A⊥ are fluid density and cross-sectional area of the particle normal to the fluid flow, respectively. ϵf is the fluid void fraction and X is given by [33]: " X ¼ 3:7−0:65 exp −
1:5− logRep 2
2 # ð21Þ
where Rep is the particle's Reynolds number. A⊥ in this work is calculated as the cross-sectional area of a cylinder with diameter of D and length of L normal to the fluid velocity (A⊥ = D L sin α). Because of general applicability and high accuracy of the Hӧlzer/ Sommerfeld model for calculation of the drag coefficient for nonspherical particles [9,16], this model is used in this study [34]: CD ¼
considered particle. The crosswise sphericity is the ratio between the cross sectional area of a spherical particle with the equivalent volume and the projected cross sectional area of the considered particle perpendicular to the flow [34]. Lift force is included because of asymmetric pressure distribution on the surface of non-spherical particles [31]. The lift force acts in the direction perpendicular to the drag force and it can be calculated as [35]: !0 ! ! h i z : v−ui !L 1 !0 ! ! ! ! F i ¼ ρ f C L A‖ ! ! z v−ui v−ui v − u i 2
0:2 1 8 1 16 1 3 1 pffiffiffiffi þ pffiffiffiffiffiffiffiffi 3=4 þ 0:42 100:4ð−logΦÞ pffiffiffiffiffiffiffi þ Rep Φ⊥ Rep Φ Φ⊥ Rep Φ
ð22Þ
where Φ and Φ⊥are regular and crosswise sphericities, respectively. The regular sphericity is the ratio between the surface area of a spherical particle with the equivalent volume and the surface area of the
ð23Þ
where A‖ is the particle area, normal to the direction of the lift force and !0 z is the particle's major axis direction. CL is the lift coefficient and it is determined as [31]: CL ¼ sin2 θ cosθ CD
ð24Þ
where, θ is the incident angle between the major axis of the particle and its relative velocity. Finally, the pressure gradient force can be calculated as [17]: !∇P F i ¼ V i ∇P
jδt j≠0
3
ð25Þ
where Vi is the volume of particle i. ∇P is the pressure gradient which is obtained from CFD solver at the particle's location. For non- spherical particles, the pressure centre in most cases does not coincide with the centre of gravity. The distance between the centre of pressure and the centre of gravity for cylindrical particles (Xcp) depends on particle length and orientation, and is obtained from [35]: X cp ¼ 0:25 L 1−e3ð1−L=dÞ cos3 θ
ð26Þ
where, L and θ are the particle's length and incident angle, respectively. As a result, the drag and lift forces induce a rotational movement !L !D due to the additional torques (X cp F and X cp F ) in Eq. (2), on the particles [14]. 2.2. Fluid flow equations The fluid flow is considered as laminar flow and solved using the volume averaged continuity and Navier–Stokes equations as [2]: ∂ ϵf ρf ∂t
! þ ∇:ϵ f ρ f v ¼ 0
! ∂ ϵf ρf v ∂t
ð27Þ
! !! ! ! þ ∇:ϵ f ρ f v v ¼ −ϵ f ∇P− ∇:ϵ f τ f þ ϵ f ρ f g − S p ð28Þ
! The fluid viscus stress tensor, τ f , is given by [4]: 2 ! ! T τ f ¼ − λ f − μ f Þ ∇:ug I−μ f ∇ v þ ð∇! vÞ 3
ð29Þ
where μf is the gas viscosity and I is the unit tensor. λf is the bulk viscosity which is neglected for gases [3]. ! The source term ( S p ) in the momentum equation which is due to momentum exchange with particles in the cell, can be calculated as [17]: kp ! !pf 1 X Sp ¼ F V cell i¼1 i
ð30Þ
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
4
F. Farivar et al. / Powder Technology xxx (2019) xxx
where kp is the number of particles residing in each CFD cell with a volume of Vcell. ϵf in each CFD cell represents the volume fraction of the gas phase and can be estimated as [36]: kp 1 X ϵ f ¼ 1− ϕV V cell i¼1 i i
ð31Þ
where ϕi ∈ [0, 1] is the volume fraction of particle i inside the CFD cell. In this work, a general Lagrangian variable Γp is mapped to an Eulerian cell Γcell and vice versa, by the weighted averaging method as [37]: Γp ¼
kc X
p
f j Γj
ð32Þ
j¼1
Γcell ¼
kp X
cell
f i Γi
ð33Þ
i¼1
where fpj is the volume fraction of the particle residing in cell j and kc is the number of overlapping cells with the particle. fcell i is the volume fraction of the cell occupied by particle i and kp is the number of particles residing in the cell.
Fig. 1. a) Multi-sphere representation of a cylindrical particle using, 5, 7, 60 and 19 spheres. b) The fluidized bed dimensions and the wall zone in the system.
3. Simulation setup and parameters The parameters and dimensions used in this simulation are presented in Table 1. The dimensions of fluid bed and particles are the same as those from the experimental work of Vollmari et al. [16]. As shown in Fig. 1a, the cylindrical particles are represented using a multi-sphere method with 5, 7, 19 and 60 overlapping spheres and named as Cyl-5, Cyl-7, Cyl-19 and Cyl-60, respectively. It is worth noting that all of the multi-sphere representations have the same flow conditions and properties in the simulations, only the number of spheres to represent each non-spherical particle is different. In the simulations, 4185 cylindrical particles are used. This number is calculated from the particle density and a mass of 0.53 kg used in the experiment. The particles are randomly packed and 4 s are allowed to settle down the particles in the bed. Fig. 1b shows the packed particles in the fluidized bed. As shown in this figure, the initial height of the bed in the simulation is roughly 0.1 m. The inlet superficial gas velocity is increased from 0.4 m/s up to 2.4 m/s in an increment of 0.2 m/s and 2 s of simulation is
Table 1 The parameters and dimensions used in the simulation. Model parameters Bed dimensions CFD cell average volume Number of CFD cells Initial pressure Coefficient of restitution (en) Friction coefficient (μn) Rolling friction coefficient (μr) CFD time step DEM time step Particle's properties Cylinder size
0.11 × 0.11 × 0.41 (m) 2.74 × 10−6 (m3) 1920 101,325 (Pa) 0.5 0.4 0.1 1 × 10−4 (s) 1 × 10−5 (s)
Poisson ratio (ν) Density Total number of particles
D = 3.9 (mm) L = 14 (mm) 0.32 764.4 (kg/m3) 4185
Fluid properties Viscosity (μf) Temperature Ambient pressure
1.82 × 10−5 (kg/m.s) 291.15 (K) 101,325 (Pa)
applied for each velocity. The same process is repeated for all of the multi-sphere representations in Fig. 1a. 3.1. Multi-sphere particle representation Representing non-spherical particles in DEM can be done using different methods. Zhong et al. [14] have classified these methods into two groups; composite-particle and single-particle methods. In the composite-particle approach complex shapes are created by combining simple geometries. The most important example is multi- sphere or glued spheres method [14]. This method was first proposed by Gallas and Sokolowski [38] to study the behaviour of granular material with non-spherical particles. In this study, multisphere method is used to represent the cylindrical particles in a gas fluidized bed. There are many possibilities to combine spherical particles with different sizes to form a cylinder-like particle, and there is no systematic method to optimize the combination. We herein aim to provide a general guideline to represent a cylindrical particle using this multi-sphere method for DEM studies. The simplest way to form a cylinder using overlapping spheres is the single-row formation. In this formation a specific number of overlapping spheres with the same diameter as the cylinder is located on the axial direction of the cylinder (along the length). Cyl-5 and Cyl-7 in Table 2 are formed using this combination. In this figure, NL is the number of spheres along the length (axial direction) and NR is the number of spheres along radial direction of the cylinder. In the single-row formation, increasing the number of spheres results in a smoother surface along the length and more accurate representation of the cylinder. To quantify this accuracy, the best parameter is the distance between centres of adjacent spheres (Dl). A higher Dl means a less accurate lengthwise representation and vice versa. Although the single row formation with a small number of spheres decreases the computational cost, this formation poorly represent the flat ends of a cylinder. One possible way to improve this is to use a multi-row formation in which spheres with a smaller size are located in multiple rows along the length of the cylinder. Cyl-60 in Table 2 is an example of this formation. In this formation, to measure the accuracy of radial representation at both ends of the cylinder, similar to the lengthwise representation Dl, Dr is defined as the longest distance between adjacent spheres forming the ends of a cylinder. For single-row
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
F. Farivar et al. / Powder Technology xxx (2019) xxx
5
Table 2 The parameters and dimensions used in the simulation. Cyl-5
Cyl-7
Cyl-60
Cy-19
5×1 (NL × NR) 2.52 × 10−3 3.90 × 10−3
7×1 (NL × NR) 1.68 × 10−3 3.90 × 10−3
10 × 6 (NL × NR) 1.33 × 10−3 1.90 × 10−3
7 + (2 × 6) (NL + 2 × NR) 1.68 × 10−3 0.95 × 10−3
3D representation
2D representation Number of spheres Dl (m) Dr(m)
formations, since there is no sphere to represent the flat ends of a cylinder Dr is the diameter of the cylinder, meaning the maximum error along the radial direction. Apart from multi-row formation, another strategy to improve the representation of the flat ends is to adopt a hybrid formation. In this hybrid formation the main body of the cylinder is represented by single-row formation and the diameter of each sphere is same as the cylinder diameter, but for both ends the multi-row formation is used and the diameter of spheres is half of the cylinder diameter. Cyl-19 in Table 2 is an example of the hybrid formation. Fig. 2 schematically represents Dr and Dl for single- row, multi- row and hybrid formations. The values of Dr and Dl for Cyl-5, Cyl7, Cyl-60 and Cyl-19 are presented in Table 2. Cyl-5 and Cyl-7 have the highest Dr with 0.0039 (m) which means they have the lowest accuracy to represent the flat ends of the cylinder. This number in Cyl-60 and Cyl-19 representations decreased to 0.0019 (m) and 0.00095(m), respectively. Dl is 0.00168 (m) for Cyl-7 and Cyl-19, and 0.00133 (m) for Cyl-60. In Section 5, the simulation results of particle dynamics in a fluidized bed using each representation are compared with experiment data to reveal the effect of shape representation on the simulation results in a multi-particle system.
solvers are run on two separate nodes in parallel and Message Passing Interface (MPI) is used to exchange data between two solvers. Each node has access to 16 shared memory processors and the OpenMP parallelization method is used separately in both CFD and DEM code. After initialization of CFD and DEM, CFD data such as pressure and gas velocity are passed to DEM. In the DEM code after defining the initial particle positions and velocities, a combination of No Binary Search (NBS) and Munjiza algorithm [17] is used to detect contacts for particle-particle and particle-wall. The contact forces are calculated using the DEM soft-sphere model based on overlapped contacts. Then the particle-fluid interaction forces acting on each particle are calculated using the CFD data. Finally, new positions and velocities of particles are obtained by solving the differential equations based on the Newton's second law (Eqs. (1) and (2)). This process repeats 10 times to match with the CFD time step. DEM data such as particles position and particle-fluid interaction forces are sent to the CFD solver, where the volume-averaging method is applied to map the DEM data into CFD cells. After calculation of bed void fraction, the fluid flow equations (Eqs. (27) and (28)) are solved to update fluid flow velocities in each cell. This process is completed at the final time step of CFD.
4. Coupled CFD-DEM algorithm
5. Results and discussions
In this work, an open source MFIX CFD solver [39] is coupled with a three dimensional in-house DEM code. The DEM code is written in FORTRAN 95 programming language and it is compatible with the CFD solver. They are compiled by gfortran compiler on a LINUX system. Fig. 3 shows the general algorithm of the CFD-DEM model. The developed CFD-DEM coupling is parallelized using a combination model for shared memory and distributed memory. CFD and DEM
5.1. Initial fluidization process The snapshots of initial fluidization process at a superficial velocity of 2.4 m/s for Cyl-7 are shown in Fig. 4a. At the initial stage (t = 0.0 s), particles are compacted in the bottom of the bed mostly in horizontal orientation. At t = 0.1 s, gas permeates through the gap between particles to create more volume voidage in the bed. The incoming gas
Fig. 2. Schematic representation of Dr and Dl for single-row, multi-row and hybrid formations.
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
6
F. Farivar et al. / Powder Technology xxx (2019) xxx
Fig. 3. The CFD-DEM model algorithm.
lifts the particles out of the bottom at t = 0.15 s and the bed height changes slightly in comparison with t = 0.1 s. After 0.2 s, the top part of the bed is fluidized, while a dense structure of particle packing remains at the bottom of the fluidized bed. Spherical particles packed in the bed are also simulated for fluidization at the same superficial velocity. Gas penetrates through gaps between spherical particles and lifts the particles at t = 0.1 s, resulting in an increase in the bed voidage in the middle of the bed. At t = 0.2 s, air partially fluidizes particles in the middle of bed and pushes the particles to the sides of the bed, and then the particles return to the bottom part of the fluidized bed, re-forming a packed structure at t = 0.3 s. The difference in the fluidization behaviour between non-spherical and spherical particles may be due to particle-particle and particlewall interactions. To further understand the fluidization process of non-spherical particles, pressure drop, particle density and particle orientation distribution are presented below by averaging the simulation results from 100 time steps (2 s of simulation with an interval of 0.02 s) for each superficial velocity.
drops for Cyl-5 and Cyl-7 are nearly identical at most of gas velocities, while the maximum difference of 26 Pa is found at the velocity of 2.4 m/s. Simulations for Cyl-5 and Cyl-7 slightly under-predict the pressure drop after the fluidization occurs at around 0.8 m/s. The simulation results for Cyl-60 and Cyl-19 have the better match with the experiment data. However, between 0.8 m/s and 0.9 m/s, where the fluidization starts, the simulation cannot capture the peak of the pressure drop. This peak is the result of interlocking effect between cylindrical particles in the bed. Since the static friction is not modelled in the DEM and only simplified constant friction factor is included in the calculations, this model cannot accurately capture this peak of pressure drop at the onset of fluidisation. Fig. 5b shows that the simulation results for the pressure drop of spherical particles with the diameter of 7 mm at different gas superficial velocities is in good agreement to the experiment. In contrast to cylindrical particle, the pressure drop for spherical particles does not drop after the fluidization. This is because spherical particles have smooth surface and no interlocking happens between them.
5.2. Pressure drop
5.3. Particle density distribution
The average pressure drop in the bed for Cyl-5, Cyl-7, Cyl-60 and Cyl19 with different gas superficial velocities is displayed in Fig. 5a. The results for each gas velocity are presented as the average data for 2 s simulation. The simulation results are in good agreement with the experiment in which the pressure at the top of the bed (height = 0.31 m) and below the flow distribution plate is measured by a pressure transducer [16]. As the superficial velocity increases from 0.4 m/s to 0.8 m/s, the pressure drop in the bed increases from around 150 Pa to 370 Pa. At 0.8 m/s, the pressure drop reaches a plateau. This is a typical pressure drop profile with the gas superficial velocity in the fluidized bed and the onset of fluidization occurs at a superficial velocity of 0.8 m/s. For different multi-spheres to represent one cylindrical particle, pressure
In Fig. 6, the average number of particles at the cross-section of different bed heights is used to present the particle spatial distribution in the bed. In this figure, the values of the predicted number density at gas superficial velocities of 1.6 m/s, 2.0 m/s and 2.4 m/s for Cyl-5, Cyl7, Cyl-60 and Cyl-19 are compared to the experiment data. Generally, the predicted data for Cyl-5, Cyl-7, Cyl-60 and Cyl-19 have a similar trend as the experimental results: a relatively higher number of particles in the bottom than on the top of the bed, suggesting more particles are accumulated in the bottom and less particles fluidize at the top. This is also verified from the snapshots of particle fluidization in Fig. 4. As the superficial velocity increases, particles are more dispersed, which is indicated by the average number density of 0.07 for 1.6 m/s down to 0.05 for 2.4 m/s at lower part of the bed in Fig. 6.
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
F. Farivar et al. / Powder Technology xxx (2019) xxx
7
Fig. 4. Paraview snapshots of the fluidized bed at a gas superficial velocity of 2.4 m/s for a) Cyl-7 for non-spherical particles with d = 3.9 mm and L = 14 mm, and b) spherical particle with d = 7 mm.
At a gas superficial velocity of 1.6 m/s (Fig. 6a), simulation results for Cyl-5 have the largest average deviation of 14.4e-3 from the experimental data. This number for Cyl-7 and Cyl-60 reduces to 7.5e-3 and 5.9e-3, respectively. Simulation results for Cyl-19 have the lowest average deviation of 5.8e-3. Fig. 6(b, c) show that all of the cylinder representations except Cyl-5, have a good match with the experimental observations at velocities of 2.0 m/s and 2.4 m/s. Simulation data for Cyl-60 and Cyl-19 have comparable accuracy with an average deviation of 5.26e-3 and 5.19e-3 at 2.0 m/s, respectively. The simulation results for the particles density distribution show that increasing the number of spheres in the single-row formation to decrease the lengthwise representation error, improves the accuracy of the model substantially. On the other hand, using both multi-row and hybrid formations to further improve the multi-sphere representation, reduce the simulation errors and improves the accuracy of the model roughly the same amount. However, in some cases such as at Vg = 2 m/s, this improvement can be marginal. 5.4. Particle orientation distribution Unlike spherical particles, particle orientation plays a significant role on cylindrical particle behaviour in the fluidized bed. In this study, the angle of 0.0° refers to the particle in the horizontal orientation aligning with the longest dimension, while the angle of 90.0° for the vertical orientation. We attempt to compare the particle orientation in the fluidized bed with experimental data. Since the experimental imaging may only capture the particles close to the experimental walls, two
simulation results are obtained: one from averaged particles in the domain (Cyl-5, Cyl-7, Cyl-60 and Cyl-19), and the other from averaged data only in the vicinity of walls (thickness of 0.01 m) with a defined wall zone shown in Fig. 1 (Cyl-5 wall zone, Cyl-7 wall zone, Cyl-60 wall zone and Cyl-19 wall zone). The number density vs the oriental angle at the superficial velocity of 1.6 m/s, 2.0 m/s and 2.4 m/s for Cyl-5, Cyl-7, Cyl-60 and Cyl-19 is presented in Fig. 7. The experimental data are also included in Fig. 7 for comparison. Simulation results for average orientation of the particles in the fluidized bed are presented in Fig. 7 for all particles in the whole bed (red lines with circular marker) as well as particles in the vicinity of the walls (green lines with triangular marker). At a superficial velocity of 1.6 m/s (Fig. 7a, d, g, j), experimental results [16] show that horizontal (between 0° and 10°) and vertical (between 80° and 90°) orientations are dominant for the particles in the bed. Particles in the horizontal orientation are relatively stable while they are likely to orient vertically (parallel to the walls) when they approach to walls. Similar wall effect on particles orientation is also reported in a previous study [8]. The number of particles decreases as the angle increases from horizontal to vertical orientation except at the angle of 90o. The simulation data for the average number of particles in the whole domain or in the wall zone follow a very similar trend as the experimental observations. However, predication for the number of particles at the full vertical orientation (90o) fails for Cyl-5 and Cyl7, while the improvement of the predications is shown in Fig. 7g for Cyl-60 and Fig. 7j for Cyl-19 in the wall zone. The poor representation
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
8
F. Farivar et al. / Powder Technology xxx (2019) xxx
Fig. 5. Pressure drop of the fluidized bed obtained from experiment [16] and simulation results for, (a) Cyl-5, Cyl-7, Cyl-60 and Cyl-19 and (b) spherical particles with d = 7 mm.
of both ends in single-row formations (Cyl-5 and Cyl-7) may lead to unstable vertical orientation (standing position), and the round ends are not able to account for accurate particle-wall interactions. The simulation data averaged in the whole domain indicate that the vertical orientation is the least preferred for particles, opposite to the experimental results. This may be due to the difference in counting particles by the experimental imaging method. The two dimensional nature of experimental image analysis suggests that particles close to the walls are more likely to be detected. An increase in the superficial velocity does not affect the trend of the number density versus the oriental angle for Cyl-5, Cyl-7, Cyl-60 and Cyl-19, while the number density with vertical orientation slightly increases by increasing the superficial velocity, which is supported from both experimental and simulation data. In all multi-sphere representations, simulation results of particles' orientation distribution for hybrid formation (Cyl-19) have the highest accuracy especially at a gas superficial velocity of 2.4 m/s. The simulation results suggest that an increasing number of spheres in the multisphere method using the same formation method can improve the accuracy of simulation results. However, the best representation is not always the one with higher number of spheres. After averaging the particle orientation of all particles on a cross section at a specific bed height, the averaged orientation at different bed heights is shown in Fig. 8. Particles preferably stay in the horizontal orientation at the bottom of the bed. Slightly above the bottom of the bed, the average angle fluctuates around 35° for the velocity of 1.6 m/s and 2.0 m/s, while for 1.0 m/s it remains stable around 20°. This behaviour can be attributed to the effect of particle-particle and particle-wall interactions on the average oriental angle along the vertical direction [40]. By increasing the gas superficial velocity, the number of particle- particle and particle- wall interaction of fluidizing particles increases, which leads to an increase in the particles' average orientation [40].
Fig. 6. Average height of the particles at gas superficial velocities of a) 1.6 (m/s), b) 2.0 (m/ s) and c) 2.4 (m/s).
5.5. The effect of torque induced by drag and lift forces Cylindrical particles experience rotation due to particle-particle and particle-wall interactions as well as drag and lift forces exerted from fluid on particles. The most stable orientation for both stationary and fluidizing particles with a cylindrical shape is horizontal (0°), and the fluidizing particles would have no preferred orientation without any influence of fluid flow on the particle orientation. In order to reveal the effect of torque induced by drag and lift forces, we conducted two separate simulations using Cyl-60: one with the assumption of Xcp = 0 in Eq. 27 and the other including this parameter in every time step (Xcp ≠ 0). The average number of particles at different orientation angles for two scenarios are shown in Fig. 9. At a gas superficial velocity of 1.6 m/s including the torque induced by fluid flow, the number density vs the oriental angle is very similar
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
F. Farivar et al. / Powder Technology xxx (2019) xxx
9
Fig. 7. Particle orientation distribution averaged through all of the particles and particles in the wall zone, for Cyl-5 (a, b, c), Cyl-7 (d, e, f), Cyl-60 (g, h, I) and Cyl-19 (j, k, l) at gas superficial velocities of 1.6 m/s (left column), 2 m/s (middle column) and 2.4 m/s (right column).
to that without the torque, which indicates that the fluid flow torque has no effect on the average orientation. At this velocity, most of the particles retain at the bottom of the bed and they are in the horizontal
orientation. Thereby, fluid flow has a negligible influence on the average orientation of the particles. However, the difference in the particle orientation distribution becomes very significant at high velocities between with and without the torque as shown in Fig. 9b and c. Since particle-fluid interaction forces (drag and lift forces) become more dominant, the torque caused by these forces becomes more significant. Therefore, at higher gas superficial velocities, the torque effect induced by drag and lift forces torque on cylindrical particle orientation should be considered. 6. Simulation run time
Fig. 8. Average orientation along the bed height for gas superficial velocities of 1.0, 1.6 and 2.0 m/s (Cyl-60).
Fig. 10 compares the computing time for 1 s simulation of fluidization process for Cyl_5, Cyl_7, Cyl_60 and Cyl-19. For the Cyl_5 representation, the total number of spheres in the simulation was 20,925 (5 × 4185) and 0.52 h is spent to simulate the process for 1 s. The number of spheres used for Cyl_7 representation was 29,295 and it takes 1.22 h to complete the simulation of the fluidization process. As expected, the simulations using Cyl_60 representation was the most computationally expensive among four representations: about 4.08 h. This number is almost 2 times higher than Cyl_19 run time at 2.6 h, while the accuracy of
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
10
F. Farivar et al. / Powder Technology xxx (2019) xxx
Fig. 9. Comparison of the average particles' orientation distribution for assumptions of Xcp = 0 (drag and lift forces do not cause rotation) and Xcp ≠ 0 (drag and lift forces cause rotation), at gas superficial velocity of, a) 1.6 m/s, b) 2 m/s and c) 2.4 m/s.
Acknowledgements This work was supported by an Australian Government Research Training Program (RTP) Scholarship; and Mayne Pharma International PTY LTD. The authors also acknowledge the supercomputing resources were provided by the Phoenix HPC service at the University of Adelaide. References
Fig. 10. Comparison of the computing time between Cyl_5, Cyl_7, Cyl_60 and Cyl-19 representations for simulation of one second of fluidization process.
the simulation results is almost same regarding particle density distribution and slightly better for Cyl-19 representation in particle orientation distribution in the fluidized bed.
7. Conclusions In this work a CFD-DEM model was developed to reveal the fluidization process of cylinder-like particles in a fluidized bed using different numbers of spheres and formation methods to represent one cylindrical particle by using the multi-sphere method. The simulation results of bed pressure drop using multi-row and hybrid formations (Cyl-60 and Cyl19) have better agreement with experimental data than those from single-row formations (Cyl-5 and Cyl-7). The average axial distribution of the particles is better mimicked by increasing the number of spheres from 5 to 7 in single-row formation. In both multi-row and hybrid formations the simulation results accuracy are slightly better than singlerow formation with 7 spheres. The preferential orientation is horizontal at all gas velocities, but the number of particles with vertical orientation increase when particle approach to walls and this increase is captured with multi-row and hybrid formations. In general, the hybrid formation results has higher accuracy and much lower computational time compared to the multi-row formation in this study. Furthermore, the torque effect induced by drag and lift forces from fluid flow on average particle orientation becomes significant at a high superficial velocity. This work also provides general guidelines for non-spherical particles representation using multi-sphere method and investigates the effect of number of spheres and formations methods on the dynamics of fluidized beds.
[1] C.T. Crowe, Multiphase Flow Handbook, CRC press, 2005. [2] M.A. van der Hoef, M. van Sint Annaland, J.A.M. Kuipers, Computational fluid dynamics for dense gas–solid fluidized beds: a multi-scale modeling strategy, Chem. Eng. Sci. 59 (22−23) (2004) 5157–5165. [3] M.A. van der Hoef, et al., Multiscale modeling of gas-fluidized beds, in: B.M. Guy (Ed.), Advances in Chemical Engineering, Academic Press 2006, pp. 65–149. [4] N.G. Deen, et al., Review of discrete particle modeling of fluidized beds, Chem. Eng. Sci. 62 (1–2) (2007) 28–44. [5] H.P. Zhu, et al., Discrete particle simulation of particulate systems: a review of major applications and findings, Chem. Eng. Sci. 63 (23) (2008) 5728–5770. [6] Z.Y. Zhou, et al., Discrete particle simulation of gas fluidization of ellipsoidal particles, Chem. Eng. Sci. 66 (23) (2011) 6128–6145. [7] J. Cai, Z. Yuan, X. Zhao, Study on two-way coupling of gas–solid two-phase flow of cylindrical particles, Powder Technol. 300 (2016) 136–145. [8] H. Ma, Y. Zhao, CFD-DEM investigation of the fluidization of binary mixtures containing rod-like particles and spherical particles in a fluidized bed, Powder Technol. 336 (2018) 533–545. [9] H. Ma, L. Xu, Y. Zhao, CFD-DEM simulation of fluidization of rod-like particles in a fluidized bed, Powder Technol. 314 (2017) 355–366. [10] M. Mohseni, B. Peters, M. Baniasadi, Conversion analysis of a cylindrical biomass particle with a DEM-CFD coupling approach, Case Stud. Therm. Eng. 10 (2017) 343–356. [11] L. Boer, et al., Experimental study on orientation and de-mixing phenomena of elongated particles in gas-fluidized beds, Powder Technol. 329 (2018) 332–344. [12] X. Chen, W. Zhong, T.J. Heindel, Orientation of cylindrical particles in a fluidized bed based on stereo X-ray particle tracking velocimetry (XPTV), Chem. Eng. Sci. 203 (2019) 104–112. [13] J. Thiry, F. Krier, B. Evrard, A review of pharmaceutical extrusion: critical process parameters and scaling-up, Int. J. Pharm. 479 (1) (2015) 227–240. [14] W. Zhong, et al., DEM/CFD-DEM modelling of non-spherical particulate systems: theoretical developments and applications, Powder Technol. 302 (2016) 108–152. [15] J.E. Hilton, L.R. Mason, P.W. Cleary, Dynamics of gas–solid fluidised beds with nonspherical particle geometry, Chem. Eng. Sci. 65 (5) (2010) 1584–1596. [16] K. Vollmari, R. Jasevičius, H. Kruggel-Emden, Experimental and numerical study of fluidization and pressure drop of spherical and non-spherical particles in a model scale fluidized bed, Powder Technol. 291 (2016) 506–521. [17] H.R. Norouzi, et al., Coupled CFD-DEM Modeling: Formulation, Implementation and Application to Multiphase Flows, John Wiley & Sons, 2016. [18] B. Liu, et al., Fluidization of non-spherical particles: sphericity, Zingg factor and other fluidization parameters, Particuology 6 (2) (2008) 125–129. [19] G. Lu, J. Third, C. Müller, Discrete element models for non-spherical particle systems: from theoretical developments to applications, Chem. Eng. Sci. 127 (2015) 425–465. [20] T. Oschmann, J. Hold, H. Kruggel-Emden, Numerical investigation of mixing and orientation of non-spherical particles in a model type fluidized bed, Powder Technol. 258 (2014) 304–323. [21] W.Q. Zhong, et al., Discrete element method simulation of cylinder-shaped particle flow in a gas-solid fluidized bed, Chem. Eng. Technol. 32 (3) (2009) 386–391. [22] B. Ren, et al., CFD-DEM simulation of spouting of corn-shaped particles, Particuology 10 (5) (2012) 562–572. [23] B. Ren, et al., Numerical simulation of spouting of cylindroid particles in a spouted bed, Can. J. Chem. Eng. 92 (5) (2014) 928–934. [24] H.P. Zhu, et al., Discrete particle simulation of particulate systems: theoretical developments, Chem. Eng. Sci. 62 (13) (2007) 3378–3396. [25] H.-G.M.J. Chen, Understanding the Discrete Element Method: Simulation of NonSpherical Particles for Granular and Multi-body Systems, John Wiley & Sons, 2014 448. [26] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Géotechnique 29 (1) (1979) 47–65.
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016
F. Farivar et al. / Powder Technology xxx (2019) xxx [27] L. Fries, et al., DEM–CFD modeling of a fluidized bed spray granulator, Chem. Eng. Sci. 66 (11) (2011) 2340–2355. [28] B.P.B. Hoomans, Granular Dynamics of Gas-Solid Two-Phase Flows, Universiteit Twente, 2000. [29] A. Džiugys, B. Peters, An approach to simulate the motion of spherical and nonspherical fuel particles in combustion chambers, Granul. Matter 3 (4) (2001) 231–266. [30] Y. Chen, J.R. Third, C.R. Müller, A drag force correlation for approximately cubic particles constructed from identical spheres, Chem. Eng. Sci. 123 (2015) 146–154. [31] M. Zastawny, et al., Derivation of drag and lift force and torque coefficients for nonspherical particles in flows, Int. J. Multiphase Flow 39 (2012) 227–239. [32] Y.Q. Feng, A.B. Yu, Assessment of model formulations in the discrete particle simulation of gas−solid flow, Ind. Eng. Chem. Res. 43 (26) (2004) 8378–8390. [33] R. Di Felice, The voidage function for fluid-particle interaction systems, Int. J. Multiphase Flow 20 (1) (1994) 153–159. [34] A. Hölzer, M. Sommerfeld, New simple correlation formula for the drag coefficient of non-spherical particles, Powder Technol. 184 (3) (2008) 361–365.
11
[35] C. Yin, et al., Modelling the motion of cylindrical particles in a nonuniform flow, Chem. Eng. Sci. 58 (15) (2003) 3489–3498. [36] B.P.B. Hoomans, et al., Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidised bed: a hard-sphere approach, Chem. Eng. Sci. 51 (1) (1996) 99–118. [37] M.v.S. Annaland, N.G. Deen, J.A.M. Kuipers, Chapter 23 Multi-level computational fluid dynamics models for the description of particle mixing and granulation in fluidized beds, in: M.J.H.A.D. Salman, J.P.K. Seville (Eds.), Handbook of Powder Technology, Elsevier Science B.V. 2007, pp. 1071–1107. [38] J.A.C. Gallas, S. Sokolowski, Grain non-sphericity effects on the angle of repose of granular material, Int. J. Mod. Phys. B 07 (09n10) (1993) 2037–2046. [39] H. Liu, D.K. Tafti, T. Li, Hybrid parallelism in MFIX CFD-DEM using OpenMP, Powder Technol. 259 (2014) 22–29. [40] F. Farivar, et al., Capturing particle-particle interactions for cylindrical fibrous particles in different flow regimes, Powder Technol. 330 (2018) 418–424.
Please cite this article as: F. Farivar, H. Zhang, Z.F. Tian, et al., CFD-DEM simulation of fluidization of multisphere- modelled cylindrical particles, Powder Technol., https://doi.org/10.1016/j.powtec.2019.11.016