CFD-DEM Simulation of Fluidization of Polyhedral Particles in a Fluidized Bed

: Fluidization of non-spherical particles is a common process in energy industries and chemical engineering. Understanding the ﬂuidization of non-spherical particles is important to guide relevant processes. There already have been numerous studies which investigate the behaviors of different non-spherical particles during ﬂuidization, but the investigations of the ﬂuidization of polyhedral particles do not receive much attention. In this study, the investigation of the ﬂuidization of polyhedral particles described by the polyhedron approach is conducted with a numerical CFD-DEM method. Experiments of the ﬂuidization of three kinds of polyhedral particles are conducted under the same condition with corresponding simulations to validate the accuracy of our CFD-DEM model. The results indicate that our CFD-DEM model with the polyhedron approach can predict the behaviors of polyhedral particles with reasonable accuracy. Fluidization behaviors of different polyhedral particles are also investigated in this study. Compared to spherical particles, the motion of polyhedral particles is stronger, and mixing degree is higher under the same ﬂuidization gas velocity.


Introduction
Fluidization is widely employed in different processes in energy industries and chemical engineering, such as mixing, biomass gasification, adsorption and pharmacy [1][2][3]. Obviously, particle shape has a great influence on fluidization behaviors [4,5]. Therefore, understanding the properties of fluidizations with different kinds of particles is important for relevant processes. Two ways are commonly applied to investigate the behaviors of fluidization, including the experimental method and the numerical technique. Compared to the experimental method, more detailed information can be obtained easily with the numerical technique. Therefore, the numerical technique has been widely applied in studies of fluidization. Among all numerical methods, the combined computational fluid dynamics (CFD)-discrete element method (DEM) has been one of the most important approaches of investigating the properties of fluidizations [6][7][8][9].
DEM has been widely used in the investigations of granular flow, which was first proposed by Cundall and Strack [10]. Tsuji et al. [11] first proposed CFD-DEM to investigate the fluidization of spherical particles in a two-dimensional fluidized bed. Subsequently, plenty investigations have verified the validity of the CFD-DEM method, and CFD-DEM has been widely accepted while investigating the fluidization of particles. However, many studies were focused on the fluidization of spherical particles in the past CFD-DEM investigations, but real particles are non-spherical mostly. It has been proved that particle shape is one of the key factors in the simulations of the investigations of granular systems [5,12]. Therefore, it is important to represent non-spherical particles accurately during the simulations. In order to describe non-spherical particles, different shape representation approaches have been proposed. With some approaches, non-spherical particles can be described by just one single particle, such as the ellipsoid approach [13][14][15], super-ellipsoid approach [6,[16][17][18] or polyhedron approach [19][20][21][22][23]. Non-spherical particles can also be described as a composite of some component particles with different numbers, such as the multi-sphere approach [24][25][26] or multi-super-ellipsoid approach [27,28].
For the CFD-DEM simulation of the fluidization of non-spherical particles, the calculation of the drag force is a technological problem. The orientation of a non-spherical particle would affect drag force greatly because of its asymmetric shape. So far, different models have been submitted to calculate the drag force which consider the influence of the shape and orientation of non-spherical particles [29][30][31][32]. In light of the previous research, Hölzer and Sommerfeld [29] put forward a drag force model of non-spherical particles, which has been widely applied in investigations of the fluidization of non-spherical particles [6,[33][34][35].
Besides the drag force, effective contact detection is also important for simulating the motion of non-spherical particles in a fluidized bed. Up to now, there have been different methods of solving the contact detection between different non-spherical particles. For the polyhedral particles which we focus on in this study, different methods have also been proposed in past research, such as the common plane method (CP) [36], the fast common plane method (FCP) [37], the shortest link method (SLM) [38], or the orientation discretization database solution [39]. One of the well-known contact methods between polyhedral particles was proposed on account of the theory of common plane method [36]. Based on the common plane method, the fast common plane method [37] and the shortest link method [38] were developed. Boon et al. [40] presented an method which could solve the contact detection between polygonal particles in 2D or polyhedral convex particles in 3D, and the contact detection was calculated using standard convex optimization procedures. Feng et al. [41] proposed a method which was based on the energy-conservation principle for elastic contact. The contact between arbitrarily shaped particles can be solved with this method.
Up to now, with the contact detection methods mentioned above, a number of DEM investigations involving polygonal and polyhedral particles have been carried out. Fraige et al. [42] investigated the flow behaviors of spherical particles and cubic particles in a rectangular hopper, and the simulations were compared with corresponding experiments. Mack et al. [43] investigated the polyhedral particle flow in a small slice hopper with experiments and simulations. Liu et al. [44] developed a bond and fracture model, and the breaking process of brittle materials was investigated with this method. Govender et al. [45] investigated the behaviors of polyhedral particles in the hoppers with different angles. The number of polyhedral particles could be 1 million via the BlazeDEM3D-GPU code. Xie et al. [23] developed a composite contact method between the polyhedron approach and super-ellipsoid approach. With this method, the contact between polyhedrons and super-ellipsoids could be solved by transforming the super-ellipsoids into polyhedrons. The accuracy of this contact method was validated by experiments and corresponding simulations, and the computation efficiency of this method was also tested. As for the investigations of fluidized beds, Oschmann et al. [46,47] investigated the influence of particle shape on mixing behaviors in a fluidized bed including cuboids. They found particle shape strongly influenced the mixing behaviors. The pressure drops of 13 kinds of particles were investigated by Vollmari et al. [48] experimentally and numerically, which included cubes and plates with different sizes. The CFD-DEM simulation results of most kinds of particles showed a good agreement when comparing with corresponding experiments. The mixing behaviors of four differently shaped particles were also investigated by Vollmari et al. [7] experimentally and numerically, and the particles were spheres, cubes, long cuboids and plates. They found that a faster mixing degree could be obtained by adding plates into the mixture in most cases.
The fluidization of polyhedral particles is a common behavior in many industrial processes, such as hydraulic engineering geological exploration [49], the wear of centrifugal pumps caused by suspended polyhedral particles [50], or the erosion in pneumatic conveying systems [51]. In existing research with the polyhedron approach, the investigations of the fluidization of polyhedral particles with the CFD-DEM method do not receive much attention. The shapes of the polyhedral particles investigated during the fluidization simulations are almost solely cuboids. As such, in this study, a CFD-DEM coupling method with the polyhedron approach is proposed to investigate the fluidization behaviors of polyhedral particles. The particles proposed in this study can be arbitrary convex shape described by the polyhedron approach. The Hölzer/Sommerfeld drag force model is employed to calculate the drag coefficient of polyhedral particles. To demonstrate our CFD-DEM model, the fluidization of three kinds of polyhedral particles including one kind of pentahedron and two kinds of hexahedrons are adopted in this paper experimentally and numerically. Then, the fluidization characteristics including the particle flow behaviors and particle mixing behaviors are investigated.

Particle Shape Representation Approach
The spherical approach and polyhedron approach are applied in this study to describe spherical particles and polyhedral particles, respectively. Polyhedral particles with triangular facets can be obtained by the import of their STL (Stereolithography) file. For a polyhedral particle with an arbitrary shape, the particle mass first needs to be detected. As mentioned above, the surfaces of the polyhedral particle are divided into many triangular facets. In order to detect the mass of the polyhedral particle, we choose a point in the particle randomly. With this point and the triangular facets, several tetrahedron elements can be obtained. The mass m i of the tetrahedron element i can be acquired, and then the mass m of the polyhedral particle can be calculated by: where N is the number of the tetrahedron elements. The particle's mass center needs to be detected before simulation. The mass m i and the mass center (x i , y i , z i ) of the tetrahedron element i can be obtained by the four vertexes of the tetrahedron, so the mass center (x 0 , y 0 , z 0 ) of the polyhedral particle can be obtained by: After the determination of mass center, the component tetrahedron elements of the polyhedral particle can be obtained, which are composed of the mass center and triangular facets. Taking a sphere as an example, a sphere with triangular facets can be seen in Figure 1. With these triangular facets and the mass center O, a sphere can be the composition of several component tetrahedron elements. Because of the complex shapes of polyhedral particles, a local coordinate system and a global coordinate system need to be employed to describe polyhedral particles accurately, which can have arbitrary positions and orientations. For example, in order to obtain the coordinate (x j ) of the vertex j in the global coordinate system, the vector (r j ) from the polyhedral particle's mass center to the vertex j first needs to be detected in the local coordinate system. Then, with the coordinate (x m ) of the particle's mass center in the global coordinate system and the transformation matrix A, the coordinate x j can be calculated by: where A is the transformation matrix from the local coordinate system to the global coordinate system, which can be obtained from the Euler angles [2].
the coordinate (xm) of the particle's mass center in the global co transformation matrix A, the coordinate xj can be calculated by: where A is the transformation matrix from the local coordin coordinate system, which can be obtained from the Euler angles

Particle Motion Equations
The motion of polyhedral particles includes translationa both of which obey Newton's laws of motion. The motion equ as: where m is mass and I denotes the inertia tensor of polyhedral particle translational velocity and particle rotational velocity gravitational acceleration. Fc, Fd and Fb are the contact force, buoyancy, respectively. The standard soft-sphere linear spring-d to calculate the contact force, and the contact torque Tc ca parameters [16,27]. The inertia tensor I of a polyhedral particle in the glob time-varying because of the asymmetric structure and particle inertia tensor I must first be calculated to decide the rotationa particle. Similar to the calculation of mass m, the inertia tensor I

Particle Motion Equations
The motion of polyhedral particles includes translational and rotational motion, both of which obey Newton's laws of motion. The motion equations can be formulated as: where m is mass and I denotes the inertia tensor of polyhedral particle. v and ω are the particle translational velocity and particle rotational velocity, respectively. g is the gravitational acceleration. F c , F d and F b are the contact force, the drag force and the buoyancy, respectively. The standard soft-sphere linear spring-dashpot model is applied to calculate the contact force, and the contact torque T c can be obtained by these parameters [16,27]. The inertia tensor I of a polyhedral particle in the global coordinate system is timevarying because of the asymmetric structure and particle motion. As a result, the inertia tensor I must first be calculated to decide the rotational motion of a polyhedral particle. Similar to the calculation of mass m, the inertia tensor I of a polyhedral particle in the local coordinate system can be calculated according to its component tetrahedron elements. The moment of inertia I xi , I yi , and I zi of the tetrahedron element i can be calculated by its density and volume, and the coordinates of the vertexes. Then, the moment of inertia in the local coordinate system can be written as:  (10) where N is the total number of the tetrahedron elements.
What should be noted is that the inertia tensor I in the local coordinate system can transform to the inertia tensor I in the global coordinate system using the inverse matrix A.

Contact Detection
The particles applied in this study are elastic and the contact interaction model of Feng et al. [23,41] is applied for contact detection with the polyhedron approach. The geometric features of a polyhedral particle include the vertexes, the edges connecting the vertexes, and the surfaces made up of edges. Because of the shape and the orientation of polyhedral particle, different contact types can be detected during the simulation, such as the contact between vertexes, the contact between vertex and surface, or the contact between some other geometric features of contacting polyhedral particles. A new polyhedron can be obtained in either contact type, which is the overlap of contacting polyhedral particles. This new polyhedron is called the overlap polyhedron in this study. Two simple polyhedral particles in contact with a typical contact type are chosen as an example, which is shown in Figure 2. The other contact types can be solved in a similar method. As Figure 2 shows, the overlap between two polyhedral particles forms a new polyhedron, P 1 and P 2 are the vertexes of Polyhedron 1 and Polyhedron 2, respectively. Naturally, the geometric features of the overlap polyhedron also include vertexes, edges, and facets, and all these features belong to the two contacting polyhedral particles. As shown in Figure 2, the features belonging to Polyhedron 1 are colored by red, and the other features belonging to Polyhedron 2 are colored by blue. In the overlap polyhedron, the outward normal unit vectors and areas of the facets belonging to Polyhedron 1 are n i and A i (i = 1, 2 in Figure 2), respectively. With these parameters, the normal unit vector of the overlap polyhedron n, which is normal to the contact plane, can be obtained by the sum of the unit normal vectors of the faces belonging to Polyhedron 1:

Contact Detection
The particles applied in this study are elastic and the co Feng et al. [23,41] is applied for contact detection with the p geometric features of a polyhedral particle include the vertexes vertexes, and the surfaces made up of edges. Because of the sh polyhedral particle, different contact types can be detected dur the contact between vertexes, the contact between vertex an between some other geometric features of contacting poly polyhedron can be obtained in either contact type, which is polyhedral particles. This new polyhedron is called the overlap Two simple polyhedral particles in contact with a typical con example, which is shown in Figure 2. The other contact types method. As Figure 2 shows, the overlap between two polyhed polyhedron, P1 and P2 are the vertexes of Polyhedron 1 and P Naturally, the geometric features of the overlap polyhedron al and facets, and all these features belong to the two contactin shown in Figure 2, the features belonging to Polyhedron 1 ar other features belonging to Polyhedron 2 are colored by blue. I the outward normal unit vectors and areas of the facets belong and Ai (i = 1, 2 in Figure 2), respectively. With these parameters the overlap polyhedron n, which is normal to the contact plan sum of the unit normal vectors of the faces belonging to Polyhe  With the same method, n , the sum of unit normal vectors of the facets belonging to Polyhedron 2, can also be calculated. As for the other contact types, the differences between them are that the overlap polyhedrons can be different. Different overlap polyhedrons correspond to different contact types. The determination of the geometric features involved in the contacted polyhedron is important to define the unit vectors n and n , which can be vertexes, edges, and facets with different amounts for both contacting polyhedral particles. The volume of the overlap polyhedron and the unit normal vector n i of each facet can be obtained with these geometric features.
The contact plane is perpendicular to the normal unit vector n and passes through the mass center of the overlap polyhedron, and the center of the overlap polyhedron is defined as the contact point. The combination method of the deepest-point method and polyhedron approach is applied to obtain the overlap δ. The deepest point in the overlap polyhedron is defined as the point of which the distance between this point and the contact plane has the maximum value on one side of the contact plane. After the detection of two deepest points in the overlap polyhedron on both sides of the contact plane, the overlap δ can be obtained by the distance between these two deepest points along the contact normal direction.

Fluid Motion Equations
The Navier-Stokes formulas are employed to describe the gas flow of a coupled gassolid two-phase system. The equations of continuity and momentum conservation can be formulated as: where ρ, u, p, and µ denote the gas density, gas velocity, gas pressure, and gas viscosity, respectively. x is the coordinates. F s represents the particle-fluid interaction force. The voidage fraction in the computational CFD cell is represented by ε, which can be stated as: where n is the number of particles located in the CFD cell, V p,i and V cell are the volumes of particle i and the corresponding CFD cell, respectively. As for the method of the calculation of the voidage fraction ε, some DEM cells will be associated with the CFD cell firstly, and only the particles residing in these DEM cells need to be considered, which can be seen in Figure 3a. The position of the mass center of the polyhedral particle is detected to determine if the polyhedral particle is in a computational CFD cell traditionally. Nevertheless, the accuracy of this algorithm is not high enough. With this method, if the particle mass center is in the CFD cell, the whole volume of the polyhedral particle will be applied to calculate the voidage fraction. As shown in Figure 3a, some particles are fully involved to calculate the voidage fraction when the mass centers of the particles are in the CFD cell. However, only some of them belong to this CFD cell. Additionally, some particles are not involved in the calculation of the voidage fraction when the mass centers of the particles are not in the CFD cell. However, some of them belong to this CFD cell. To solve this problem, we improved this search method. As mentioned above, an arbitrary polyhedral particle can be composed of several component tetrahedron elements with different amounts. Then, the mass center of each tetrahedron element is detected to decide if the tetrahedron element is in this CFD cell. By detecting the position of each component tetrahedron element rather than the whole particle, the calculation of the voidage fraction can be more accurate. The improved method of searching tetrahedron elements can be seen in Figure 3b. With this method, the voidage fraction can be calculated with part of the particle which is not fully involved in the CFD cell, and the computational accuracy can be improved. A similar method was used in our previous study with a good result [6]. of each tetrahedron element is detected to decide if the tetrahedron element is in this CFD cell. By detecting the position of each component tetrahedron element rather than the whole particle, the calculation of the voidage fraction can be more accurate. The improved method of searching tetrahedron elements can be seen in Figure 3b. With this method, the voidage fraction can be calculated with part of the particle which is not fully involved in the CFD cell, and the computational accuracy can be improved. A similar method was used in our previous study with a good result [6].

Gas-Particle Interaction Force
Several types of interaction force are involved in the gas-particle system, and only buoyancy and drag force are emphasized in this study on account of the modest rotational and translational motion of polyhedral particles in the fluidized bed. The lift forces can be neglected because of the low gas velocity. Ignoring the influence of the circumambient particles, the formula of the drag force imposed on an individual particle can be written as: (15) where ε and CD represent the voidage fraction and drag coefficient, respectively. u and v stand for the gas velocity and the particle velocity, respectively. A is the cross-sectional area of the sphere of which the volume is equal to the considered polyhedral particle. In fact, the ambient particles will affect the drag force. As a result, Equation (15) was modified by Di Felice [52]:

Gas-Particle Interaction Force
Several types of interaction force are involved in the gas-particle system, and only buoyancy and drag force are emphasized in this study on account of the modest rotational and translational motion of polyhedral particles in the fluidized bed. The lift forces can be neglected because of the low gas velocity. Ignoring the influence of the circumambient particles, the formula of the drag force imposed on an individual particle can be written as: where ε and C D represent the voidage fraction and drag coefficient, respectively. u and v stand for the gas velocity and the particle velocity, respectively. A is the cross-sectional area of the sphere of which the volume is equal to the considered polyhedral particle. In fact, the ambient particles will affect the drag force. As a result, Equation (15) was modified by Di Felice [52]: where Re p is the particle Reynolds number and d is the diameter of the sphere of which the volume is equal to the considered polyhedral particle.
Because of the complex shapes and non-symmetry of polyhedral particles, the calculation of drag force is much more complicated for polyhedral particles compared to spherical particles. Different drag force models have been devised for non-spherical particles. The Hölzer/Sommerfeld drag force model [29] is employed in this study to calculate the drag coefficient of polyhedral particles, which has been extensively applied in the simulation of the fluidization with non-spherical particles [6,[33][34][35]. The Hölzer/Sommerfeld drag force model [29] can be formulated as: φ denotes the regular sphericity, which is the surface area ratio of the volume equivalent sphere to the corresponding particle. φ ⊥ represents the crosswise sphericity, which can be acquired by the area ratio of the cross section of the volume equivalent sphere to the projective plane of the corresponding particle against the direction of the fluid flow.
What should be noted is that in this study, the traditional drag force model is applied for spherical particles, which can be formulated as: The equation of the calculation of the buoyancy F b can be calculated as: where V p is the volume of the corresponding particle.
After the calculations of the forces mentioned above, the gas-particle interaction force F s in a CFD cell can be obtained as: where F d,i and F b,i represent the drag force and the buoyancy imposing on particle i, respectively. n is the number of the particles in the corresponding CFD cell.

Solution and Simulation Conditions
Four kinds of particles are adopted to investigate the fluidization of polyhedral particles in this study, which include spherical particles and three kinds of polyhedral particles. The governing equations are solved by the SIMPLEC method and the QUICK differencing scheme, and the fluid motion is solved by the pressure-based implicit integration method, which is solved by the software Fluent. As for the particle motion, it is solved by the linear spring-dashpot contact model and the explicit time integration method, which is realized through a self-developed DEM code (DEMSLab), and a coupling scheme is used between CFD and DEM. More detailed information of the CFD-DEM coupling method adopted in this study can be consulted in our previous studies [3,6,33].
The schematic diagram of the fluidized bed is shown in Figure 4. The sizes of the fluidized bed are 0.2 m, 0.03 m, and 1 m in width, thickness, and height, respectively. The fluid, which is air in this study, passes in from the bottom of the vessel along the direction of z-axis. The four side walls of the vessel are set to be stationary no-slip boundaries, which are the front, rear, right, and left walls. The size of the CFD mesh is 2 cm, 3 cm, and 2 cm in the x, y, and z directions, respectively. The polyhedral particles used in this study include one kind of pentahedron and two kinds of hexahedrons. The shape of the pentahedron is a triangular prism, of which the triangle faces are equilateral triangles, and the shapes of two hexahedrons are cuboids. These polyhedral particles are modeled by the polyhedron approach and the spherical particles are described by the spherical approach. The geometrical shapes of the polyhedral particles can be seen in Figure 5, which are named Polyhedron 1, Polyhedron 2, and Polyhedron 3, respectively. Their detailed geometric parameters are summarized in Table 1. pentahedron is a triangular prism, of which the triangle faces are equilateral and the shapes of two hexahedrons are cuboids. These polyhedral particles are by the polyhedron approach and the spherical particles are described by the approach. The geometrical shapes of the polyhedral particles can be seen in which are named Polyhedron 1, Polyhedron 2, and Polyhedron 3, respective detailed geometric parameters are summarized in Table 1.      Particles are first dropped into the fluidized bed under gravity. All particles are packed with random locations and orientations, and the initial kinetic energies are all zero, allowing free fall and random packing. A total of 3600 particles are applied for each case. The four kinds of particles have the same densities (480 kg/m 3 ), and their equivalent diameters are 6.8 mm. Both the frictional coefficient and the restitution coefficient between particles are 0.5. The restitution coefficient between particle and wall is 0.5 and the frictional coefficient is 0.4. The air flow into the model from the end of the vessel and the superficial velocity starts from 0.4 m/s and ends at 2.0 m/s, of which the step is 0.1 m/s. What should be noted is that the superficial velocity is fixed during each simulation. The time step is 8 × 10 −6 s and the detailed parameters are listed in Table 2.

Validation of the CFD-DEM Model
Corresponding fluidized bed experiments of polyhedral particles are carried out to validate our CFD-DEM model. The sizes of the equipment are same as the simulations, of which the dimensions are 0.2 × 0.03 × 1.0 m (width × thickness × height). The walls of the vessel are made of the transparent polymethyl methacrylate and a high-speed camera is placed in front of the vessel to capture the behaviors of the polyhedral particles during experiments. The gas flow rate is measured with a vortex shedding flowmeter. In order to assure an even gas flow, a porous gas distributor with 20% gas passage is installed at the bottom of the vessel. The pressure drop is measured by two pressure tapping points, one point is set near the outlet and the other one is slightly under the gas distributor. Corresponding simulations are carried out for the validation of the CFD-DEM model. The parameters used in these simulations can be seen in Table 2. Figure 5 presents the snapshots of the initial packed bed in simulations and experiments, and the heights of the initial packed beds are nearly the same between experiments and corresponding simulations. Some snapshots during the experiments and simulations of three kinds of polyhedral particles can be seen in Figures 6-8, of which the gas superficial velocity is 1.6 m/s. The simulated flow behaviors of different polyhedral particles can be seen are similar as the corresponding experiments with reasonable accuracy from visual observation. For further validation, the pressure drops changing with the superficial velocity are compared between experiments and simulations. As shown in Figure 9, the results calculated from simulations match well with the pressure drops of corresponding experiments. As for the minimum fluidization velocity, it also presents good agreement with comparative accuracy between simulations and corresponding experiments. The range of the minimum fluidization velocities can be found is between 0.8 m/s and 1.0 m/s for all simulations and experiments, which predicts well between the experiments and the corresponding simulations.      Furthermore, the time-averaged bed heights of three polyhedral particles during 8-16 s can be seen in Figure 10, and the superficial velocities are 1.6 m/s and 1.8 m/s, respectively. The time-based standard deviations of heights are denoted by the error bars. The method of Mahajan et al. [53] is employed to calculate the bed height in this study. From Figure 10 we can see that the time-averaged bed heights become higher for all three kinds of polyhedral particles as the gas superficial velocity increases. In general, the time-averaged bed heights of experiments are comparable with corresponding simulations. As a result, the above results indicate that the CFD-DEM coupling method can predict the behaviors of polyhedral particles with reasonable accuracy.  Furthermore, the time-averaged bed heights of three polyhedral particles during 8-16 s can be seen in Figure 10, and the superficial velocities are 1.6 m/s and 1.8 m/s, respectively. The time-based standard deviations of heights are denoted by the error bars. The method of Mahajan et al. [53] is employed to calculate the bed height in this study. From Figure 10 we can see that the time-averaged bed heights become higher for all three kinds of polyhedral particles as the gas superficial velocity increases. In general, the time-averaged bed heights of experiments are comparable with corresponding simulations. As a result, the above results indicate that the CFD-DEM coupling method can predict the behaviors of polyhedral particles with reasonable accuracy. Furthermore, the time-averaged bed heights of three polyhedral particles during 8-16 s can be seen in Figure 10, and the superficial velocities are 1.6 m/s and 1.8 m/s, respectively. The time-based standard deviations of heights are denoted by the error bars. The method of Mahajan et al. [53] is employed to calculate the bed height in this study. From Figure 10 we can see that the time-averaged bed heights become higher for all three kinds of polyhedral particles as the gas superficial velocity increases. In general, the timeaveraged bed heights of experiments are comparable with corresponding simulations. As a result, the above results indicate that the CFD-DEM coupling method can predict the behaviors of polyhedral particles with reasonable accuracy.

Fluidization Behaviors of Different Particles
Besides the simulations mentioned above, the simulations of the fluidization with the volume equivalent sphere are also included for the purpose of comparison. Figure 11 shows the initial packed beds of four kinds of particles, we can see that the fill height of Polyhedron 1 is the highest of the four packed beds and the height of Polyhedron 3 is the second highest. The fill heights and the voidage fractions are listed in Table 3. For Polyhedron 1, its shape is a triangular prism, which is more irregular than the other two cuboid polyhedral particles. As such, for the initial pack, the void fraction of Polyhedron 1 is higher than the others, and the fill height is the highest. As for Polyhedron 2 and Polyhedron 3, the aspect ratio of Polyhedron 3, which is defined as the length ratio of the height to width in this study, is larger than that of Polyhedron 2. As a result, the void fraction of Polyhedron 3 is larger than Polyhedron 2, which leads to a higher fill height for Polyhedron 3 compared to Polyhedron 2.

Fluidization Behaviors of Different Particles
Besides the simulations mentioned above, the simulations of the fluidization with the volume equivalent sphere are also included for the purpose of comparison. Figure 11 shows the initial packed beds of four kinds of particles, we can see that the fill height of Polyhedron 1 is the highest of the four packed beds and the height of Polyhedron 3 is the second highest. The fill heights and the voidage fractions are listed in Table 3. For Polyhedron 1, its shape is a triangular prism, which is more irregular than the other two cuboid polyhedral particles. As such, for the initial pack, the void fraction of Polyhedron 1 is higher than the others, and the fill height is the highest. As for Polyhedron 2 and Polyhedron 3, the aspect ratio of Polyhedron 3, which is defined as the length ratio of the height to width in this study, is larger than that of Polyhedron 2. As a result, the void fraction of Polyhedron 3 is larger than Polyhedron 2, which leads to a higher fill height for Polyhedron 3 compared to Polyhedron 2.

Fluidization Behaviors of Different Particles
Besides the simulations mentioned above, the simulations of the fluidization with the volume equivalent sphere are also included for the purpose of comparison. Figure 11 shows the initial packed beds of four kinds of particles, we can see that the fill height of Polyhedron 1 is the highest of the four packed beds and the height of Polyhedron 3 is the second highest. The fill heights and the voidage fractions are listed in Table 3. For Polyhedron 1, its shape is a triangular prism, which is more irregular than the other two cuboid polyhedral particles. As such, for the initial pack, the void fraction of Polyhedron 1 is higher than the others, and the fill height is the highest. As for Polyhedron 2 and Polyhedron 3, the aspect ratio of Polyhedron 3, which is defined as the length ratio of the height to width in this study, is larger than that of Polyhedron 2. As a result, the void fraction of Polyhedron 3 is larger than Polyhedron 2, which leads to a higher fill height for Polyhedron 3 compared to Polyhedron 2.   Figure 12 shows the pressure drops of four particles which change with the superficial velocity, and the error bar denotes the time-based standard deviation of the pressure drop. It is generally known that during the packed bed stage, the fluidized bed keeps stable and the pressure drop has little fluctuation. As the gas superficial velocity increases, the fluidized bed starts to expand and bubbles will occur in the incipient fluidization stage, which leads to a strong fluctuation of pressure drop. From Figure 12, the minimum fluidization velocities of all kinds of particles can be acquired, and the differences of the minimum fluidization velocity between Polyhedral 1 and Polyhedron 2 are small. The minimum fluidization velocity of Polyhedron 3 is a little larger than the velocities of Polyhedron 1 and Polyhedron 2.   Figure 12 shows the pressure drops of four particles which change with the superficial velocity, and the error bar denotes the time-based standard deviation of the pressure drop. It is generally known that during the packed bed stage, the fluidized bed keeps stable and the pressure drop has little fluctuation. As the gas superficial velocity increases, the fluidized bed starts to expand and bubbles will occur in the incipient fluidization stage, which leads to a strong fluctuation of pressure drop. From Figure 12, the minimum fluidization velocities of all kinds of particles can be acquired, and the differences of the minimum fluidization velocity between Polyhedral 1 and Polyhedron 2 are small. The minimum fluidization velocity of Polyhedron 3 is a little larger than the velocities of Polyhedron 1 and Polyhedron 2. According to Equations (15)- (17), it is clear that decreasing the bed voidage fraction and increasing the projected face area will cause the increase of drag force, and the particle motion can be promoted at last. For the particles mentioned above, the voidage fractions in the packed bed are summarized in Table 3. The orientation of non-spherical particles has influence on the projected face area. In this study, the Euler angles are employed to describe particle orientation. The angle between the z-axis and z′-axis is the nutation angle (θ), which can be seen in the inset of Figure 13a. The angle between the x-axis and the intersection line between the x-y plane and x′-y′ plane is the precession angle (ψ), which can be seen in the inset of Figure 13b.The ranges of the nutation angle and precession angle are both 0°-180° because of the symmetry of polyhedrons. Figure 13 shows the distributions of the nutation angle and precession angle of three polyhedrons in the packed bed. It can be seen that the orientations are similar on all three polyhedrons According to Equations (15)- (17), it is clear that decreasing the bed voidage fraction and increasing the projected face area will cause the increase of drag force, and the particle motion can be promoted at last. For the particles mentioned above, the voidage fractions in the packed bed are summarized in Table 3. The orientation of non-spherical particles has influence on the projected face area. In this study, the Euler angles are employed to describe particle orientation. The angle between the z-axis and z -axis is the nutation angle (θ), which can be seen in the inset of Figure 13a. The angle between the x-axis and the intersection line between the x-y plane and x -y plane is the precession angle (ψ), which can be seen in the inset of Figure 13b.The ranges of the nutation angle and precession angle are both 0 • -180 • because of the symmetry of polyhedrons. Figure 13 shows the distributions of the nutation angle and precession angle of three polyhedrons in the packed bed. It can be seen that the orientations are similar on all three polyhedrons in the packed bed, so the surface areas in Table 3 can be used to indicate the differences of the projected face areas between these polyhedrons. As shown in Table 3, because of the relatively large voidage fraction and small surface area, the minimum fluidization velocity of Polyhedron 3 is a little higher compared to Polyhedron 1 and Polyhedron 2. in the packed bed, so the surface areas in Table 3 can be used to indicate the differences of the projected face areas between these polyhedrons. As shown in Table 3, because of the relatively large voidage fraction and small surface area, the minimum fluidization velocity of Polyhedron 3 is a little higher compared to Polyhedron 1 and Polyhedron 2. As for the differences between polyhedral particles and spherical particles, particle shape can have influence on the motion of different particles and Figure 14 shows the snapshots of the behaviors of these four kinds of particles at 1.2 m/s, in which the snapshots are obtained at the same time for each case respectively. It can be seen that more bubbles appear and bed heights are higher for polyhedral particles when comparing with spherical particles, especially for Polyhedron 1. In other words, the motion of Polyhedron 1 is stronger than spherical particles at the same superficial velocity. It can be seen from Table 3 that the surface area of Polyhedron 1 is much higher than the surface area of the spherical particle, which could be the cause of the stronger motion of Polyhedron 1 compared to the spherical particles at the same superficial velocity. As for the differences between polyhedral particles and spherical particles, particle shape can have influence on the motion of different particles and Figure 14 shows the snapshots of the behaviors of these four kinds of particles at 1.2 m/s, in which the snapshots are obtained at the same time for each case respectively. It can be seen that more bubbles appear and bed heights are higher for polyhedral particles when comparing with spherical particles, especially for Polyhedron 1. In other words, the motion of Polyhedron 1 is stronger than spherical particles at the same superficial velocity. It can be seen from Table 3 that the surface area of Polyhedron 1 is much higher than the surface area of the spherical particle, which could be the cause of the stronger motion of Polyhedron 1 compared to the spherical particles at the same superficial velocity.
in the packed bed, so the surface areas in Table 3 can be used to indicate the differences of the projected face areas between these polyhedrons. As shown in Table 3, because of the relatively large voidage fraction and small surface area, the minimum fluidization velocity of Polyhedron 3 is a little higher compared to Polyhedron 1 and Polyhedron 2. As for the differences between polyhedral particles and spherical particles, particle shape can have influence on the motion of different particles and Figure 14 shows the snapshots of the behaviors of these four kinds of particles at 1.2 m/s, in which the snapshots are obtained at the same time for each case respectively. It can be seen that more bubbles appear and bed heights are higher for polyhedral particles when comparing with spherical particles, especially for Polyhedron 1. In other words, the motion of Polyhedron 1 is stronger than spherical particles at the same superficial velocity. It can be seen from Table 3 that the surface area of Polyhedron 1 is much higher than the surface area of the spherical particle, which could be the cause of the stronger motion of Polyhedron 1 compared to the spherical particles at the same superficial velocity.

Particle Mixing in the Fluidized Bed
Understanding the behaviors of particle mixing is important for relevant processes. Figure 14 presents the snapshots of particle behaviors for spheres, Polyhedron 1, Polyhe- dron 2, and Polyhedron 3, respectively, of which the superficial velocity is 1.2 m/s. The fluidized bed is divided into four equal parts along the direction of the z-axis, of which different parts are colored by black and white. It is clear that polyhedral particles have a higher mixing degree than the spherical particles at the same superficial velocity. This is because particle orientation has an influence on polyhedral particles but spheres are not influenced by orientation. As a result, the particle fluid forces change more drastically for polyhedral particles compared to spheres, which can cause higher mixing degrees [7,46,47].
Although the particles at the bottom of the bed form a cluster and have a tendency to rise, some particles at the bottom of the vessel nearly do not move up at all. This is particularly obvious for the spherical particles in Figure 14a. Figure 15 shows the snapshots of the particle behaviors at 8 s for all four kinds of particles, of which the superficial velocity is 1.2 m/s. It is clear that spherical particles have a minimum mixing degree visually. As shown in Figure 15, for spherical particles, there are more black particles still staying at the bottom of the bed. For the three kinds of polyhedral particles, particles in black stay at the bottom of the bed in various degrees, and for Polyhedron 3 this degree is a little larger than the degrees of Polyhedron 1 and Polyhedron 2, visually. More black particles can be seen in the upside of the bed for Polyhedron 1 and Polyhedron 2.

Particle Mixing in the Fluidized Bed
Understanding the behaviors of particle mixing is important for relevant processes. Figure 14 presents the snapshots of particle behaviors for spheres, Polyhedron 1, Polyhedron 2, and Polyhedron 3, respectively, of which the superficial velocity is 1.2 m/s. The fluidized bed is divided into four equal parts along the direction of the z-axis, of which different parts are colored by black and white. It is clear that polyhedral particles have a higher mixing degree than the spherical particles at the same superficial velocity. This is because particle orientation has an influence on polyhedral particles but spheres are not influenced by orientation. As a result, the particle fluid forces change more drastically for polyhedral particles compared to spheres, which can cause higher mixing degrees [7,46,47].
Although the particles at the bottom of the bed form a cluster and have a tendency to rise, some particles at the bottom of the vessel nearly do not move up at all. This is particularly obvious for the spherical particles in Figure 14a. Figure 15 shows the snapshots of the particle behaviors at 8 s for all four kinds of particles, of which the superficial velocity is 1.2 m/s. It is clear that spherical particles have a minimum mixing degree visually. As shown in Figure 15, for spherical particles, there are more black particles still staying at the bottom of the bed. For the three kinds of polyhedral particles, particles in black stay at the bottom of the bed in various degrees, and for Polyhedron 3 this degree is a little larger than the degrees of Polyhedron 1 and Polyhedron 2, visually. More black particles can be seen in the upside of the bed for Polyhedron 1 and Polyhedron 2. To further investigate the particle mixing, the fluidized bed is divided into two equal parts along the direction of the z-axis, as shown in the inset of Figure 16a. The Lacey mixing index [54] is employed to quantify the mixing process. Figure 16a shows the variation of the mixing index changing with time for all kinds of particles of which the superficial velocity is 1.2 m/s. As mentioned above, the mixing degree of spherical particles is the lowest among these four kinds of particles, which can also be observed from Figure 16a. The mixing degree of Polyhedron 3 is a little lower than the mixing To further investigate the particle mixing, the fluidized bed is divided into two equal parts along the direction of the z-axis, as shown in the inset of Figure 16a. The Lacey mixing index [54] is employed to quantify the mixing process. Figure 16a shows the variation of the mixing index changing with time for all kinds of particles of which the superficial velocity is 1.2 m/s. As mentioned above, the mixing degree of spherical particles is the lowest among these four kinds of particles, which can also be observed from Figure 16a. The mixing degree of Polyhedron 3 is a little lower than the mixing degrees of Polyhedron 1 and Polyhedron 2. The mixing degree is influenced by particle motion. As a result, the spherical particles demonstrate a lower mixing degree because of their more evenly distributed forces. Furthermore, Figure 16b shows the variation between the Lacey mixing index and the time while the superficial velocity is 1.6 m/s. At the beginning of the process, the mixing degrees of all kinds of particles are almost equal because of the strong particle motion caused by the sudden high superficial velocity. As the simulations continue, the mixing degree of spherical particles remains lower than the mixing degrees of the other three polyhedral particles. their more evenly distributed forces. Furthermore, Figure 16b shows the variation between the Lacey mixing index and the time while the superficial velocity is 1.6 m/s. At the beginning of the process, the mixing degrees of all kinds of particles are almost equal because of the strong particle motion caused by the sudden high superficial velocity. As the simulations continue, the mixing degree of spherical particles remains lower than the mixing degrees of the other three polyhedral particles.

Conclusions
In this study, a numerical CFD-DEM investigation of polyhedral particles described by the polyhedron approach in a fluidized bed is studied. The mass and the inertia tensor of a polyhedral particle can be calculated by the accumulation of the masses and the inertia tensors of its tetrahedron elements. Contact detection with the polyhedron approach is based on the combination of the overlapping volume method and the deepest point method. The tetrahedron elements of the polyhedral particle are applied to calculate the voidage fraction. Based on the comparison results between experiments and corresponding simulations, the established CFD-DEM method for the polyhedron approach can reproduce the fluidization of polyhedral particles with reasonable accuracy. The fluidization behaviors of polyhedral particles are also investigated with this CFD-DEM method. Comparing with spherical particles, the motion for polyhedral particles is stronger and the mixing degree is higher under the same fluidization gas velocity. In summary, this study provides an effective CFD-DEM coupling method for the polyhedron approach which can be applied for future investigations involving the fluidization of non-spherical particles encountered in industrial applications.

Conclusions
In this study, a numerical CFD-DEM investigation of polyhedral particles described by the polyhedron approach in a fluidized bed is studied. The mass and the inertia tensor of a polyhedral particle can be calculated by the accumulation of the masses and the inertia tensors of its tetrahedron elements. Contact detection with the polyhedron approach is based on the combination of the overlapping volume method and the deepest point method. The tetrahedron elements of the polyhedral particle are applied to calculate the voidage fraction. Based on the comparison results between experiments and corresponding simulations, the established CFD-DEM method for the polyhedron approach can reproduce the fluidization of polyhedral particles with reasonable accuracy. The fluidization behaviors of polyhedral particles are also investigated with this CFD-DEM method. Comparing with spherical particles, the motion for polyhedral particles is stronger and the mixing degree is higher under the same fluidization gas velocity. In summary, this study provides an effective CFD-DEM coupling method for the polyhedron approach which can be applied for future investigations involving the fluidization of non-spherical particles encountered in industrial applications.