Voronoi analysis of the packings of non-spherical particles

Chemical Engineering Science 153 (2016) 330–343

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

Voronoi analysis of the packings of non-spherical particles Kejun Dong a,n, Chuncheng Wang b, Aibing Yu b a b

Institute for Infrastructure Engineering, Western Sydney University, Penrith, NSW 2751, Australia Laboratory for Simulation and Modeling of Particulate Systems, Department of Chemical Engineering, Monash University, Clayton, VIC 3800, Australia

H I G H L I G H T S

G R A P H I C A L

A B S T R A C T


A general method to perform Voronoi analysis of non-sphere packings proposed.
Voronoi analysis conducted on the packings of ellipsoids and cylinders.
Effects of particle shape and friction on Voronoi cell properties studied.
Universal relationships between microstructure properties and porosity found.

art ic l e i nf o

a b s t r a c t

Article history: Received 1 April 2016 Received in revised form 28 June 2016 Accepted 10 July 2016 Available online 12 July 2016

We present a structural analysis of the packings of identical non-spherical particles based on Voronoi cells. The packings are generated by discrete element method (DEM) simulations. The particles include axisymmetric ellipsoidal particles from oblates to prolates and cylindrical particles from disks to rods. The Voronoi cells are constructed by virtue of space discretization and surface reconstruction, which is shown to be universal for different shapes. The effects of particle aspect ratio and sliding friction coefficient on the properties of Voronoi cells, including the reduced volume, reduced surface area and sphericity, are quantified. The reduced volume and surface area are found to observe log-normal distributions, while their mean values and standard deviations have different dependencies on particle shape and friction. By analyzing the correlations and using inherent relationships between different Voronoi cell properties, we establish a group of universal equations to predict these distributions according to particle sphericity and overall packing fraction. Such findings can not only improve our understanding on the packings of non-spherical particles but also provide a basis for evaluating the transport properties and advancing the statistical mechanics theory for granular matter composed of non-spherical particles. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Particle packing Particle shape Structural analysis Voronoi analysis Discrete element method

1. Introduction Particle packing is widely encountered in various industrial processes. In many of these processes, the structure of the packing dominates the process performance as it determines the related transport properties (Dullien, 1992 ), such as the permeability of a filter cake (Dong et al., 2009b; Moreno-Atanasio et al., 2010), the heat conductivity of a packed bed heat exchanger (Cheng and Yu, 2013),

n

Corresponding author. E-mail address: [email protected] (K. Dong).

http://dx.doi.org/10.1016/j.ces.2016.07.013 0009-2509/& 2016 Elsevier Ltd. All rights reserved.

and the electro-conductivity of a dusk cake in an electrostatic precipitator (Guo et al., 2013; Yang et al., 2013). However, to understand and model the packing structure is still a scientific challenge. The difficulties lie in the following aspects: first, a packing is normally disordered with no simple way to describe (Bideau and Hansen, 1993; Tian et al., 2011); secondly, it includes both the particle network and pore network; thirdly, it is affected by a number of parameters including the properties of particles, the environment, and the packing method. Previous studies focused on the modeling of porosity (or packing fraction) (German, 1989; Zou and Yu, 1996a) and its relationship with the transport properties (Dullien, 1992 ), which can only be empirical or semi-empirical as porosity is just a macroscopic

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

parameter. In the last two decades, with the developments of experimental techniques such as X-ray tomography (Aste, 2006; Moreno-Atanasio et al., 2010; Xia et al., 2014) and numerical methods such as discrete element method (DEM) (Zhu et al., 2007, 2008), packing structures have been characterized at the particle-scale in much more details. Such characterizations have been applied to sphere packings first, including the packings of monosized spheres and multisized spheres (Yi et al., 2011, 2012, 2015), cohesionless and cohesive spheres (Yang et al., 2002, 2013), and in vacuum or liquids (Dong et al., 2009a, 2006). The characterization methods have also been developed from the two-body to many-body scales (Anikeenko and Medvedev, 2007; Tian et al., 2014). These studies have led to a significant progress in the modeling of the structures of sphere packings (Song et al., 2008). In recent years, the packings of non-spherical particles have attracted increasing attention in different disciplines. Controversial effects of particle shape on the packing structures have been shown even for the very regular non-spherical particles, such as the ellipsoids (Donev et al., 2004) and regular polyhedrons (Torquato and Jiao, 2009), which are related to applications of not only granular particles but also nano particles (Glotzer and Engel, 2011). Theoretically this is because a non-spherical particle has additional rotational degrees of freedom than a sphere, which enriches the structures of the packings. Similar to sphere packings, computer simulations have also been used to obtain the packing structures of non-spherical particles. However to model nonspherical particles in DEM is much more complicated than spherical particles. In the literatures, there are various methods used. Some methods are for specific shapes such as the geometric potential method for ellipsoidal (Zhou et al., 2011) and superquadric particles (Delaney and Cleary, 2010), and the common-plane method (Vorobiev, 2012) or GJK method (Wachs et al., 2012) for polyhedral particles, while other methods use composite (Favier et al., 1999) or cluster particles (Guo et al., 2015) to represent different shapes. There are also methods that combine DEM with the finite element method (FEM) for non-spherical shapes (Latham and Munjiza, 2004; Xiang et al., 2009). One can refer to several review papers on this topic (Džiugys and Peters, 2001; Lu et al., 2015). Recently, we have developed a novel method based on orientation discretization, which can simulate non-spherical particles in a simple, uniform and fast way in DEM (Dong et al., 2015). By these methods, the packings of particles of different shapes have been simulated, including ellipsoids/super-ellipsoids (Delaney and Cleary, 2010; Donev et al., 2004; Zhou et al., 2011), cylinders/sphero-cylinders (Kodam et al., 2010), polyhedrons (Li et al., 2013; Torquato and Jiao, 2009), and so on. The resulting packings have been characterized by different parameters, including the radial distribution function (RDF), coordination numbers (CN), nematic and orientational orders (Li et al., 2013; Zhao et al., 2012), interfacial volume fraction (Xu et al., 2015, 2016), and other structural parameters. In particular, the coordination number has been found to be able to predict by the isostatic condition for ideal frictionless particles (Donev et al., 2004), whereas for granular particles the packings can be hyperstatic due to the effect of friction (Schaller et al., 2015). These analyses have improved our understanding on the packing structures of non-spherical particles. However, compared to those of sphere packings, the structural characterizations of the non-sphere packings are far less comprehensive. Notably one important characterization, Voronoi analysis, is still limitedly conducted on non-sphere packings. Voronoi analysis can effectively characterize both the particle network and pore network of a packing, which hence is very helpful to the modeling of transport properties, particularly for those related to the fluid-packed-bed interactions (Cheng and Yu, 2013; Rong et al., 2013). Voronoi analysis has also been used in building the statistical mechanics theory for particle packing

331

(Baule et al., 2013; Song et al., 2008) since it can precisely characterize the local space around each particle. Such kind of analysis has been widely conducted for sphere packings (Yang et al., 2002; Yi et al., 2015), but just a few for non-spherical particles (Baule et al., 2013; Luchnikov et al., 1999; Schaller et al., 2015; Xia et al., 2014). This is due to not only the lacking of the data of non-sphere packings but also the difficulty in constructing the Voronoi cells for non-sphere packings, which are more complicated than sphere packings. Nevertheless, recent experimental studies have shown some interesting findings in this aspect. For example, by using X-ray tomography, Schaller et al. (2015) reported a Voronoi analysis of the packings of frictional oblates, and demonstrated the local packing fraction distribution can be correlated with overall packing fraction, while Xia et al. (2014) proposed that the particle asphericity can be treated as polydispersity effect for the packings of oblates. These studies were focused on the oblates, and the analyses were mainly on Voronoi cell volumes related to the modeling of packing fraction, but other Voronoi cell properties were less discussed, which will also be needed in the modeling of transport properties of packings. In this paper, we perform the Voronoi analysis on the packings of 3D non-spherical particles, including both ellipsoidal particles (oblates and prolates) and cylindrical particles (disks and rods). The packings are generated using a recent developed method for discrete modeling of non-spherical particles (Dong et al., 2015). The Voronoi cells for the non-sphere packings are obtained by using space discretization combined with mesh reconstruction, which will be described in Section 2. In Section 3, the effects of particle shape and friction on the packing structures will be quantified in terms of the properties of the Voronoi cells, including the reduced volume, reduced surface area and sphericity. The results reveal the complicated effects of particle shape and friction on the packing structures, whereas we also discover some universal relationships between the microstructure properties of Voronoi cells and macroscopic parameter, packing fraction. The findings can improve our understanding on the packings of non-spherical particles and provide structure models for evaluating transport properties of these packing and advancing statistical mechanics theory for non-sphere packings.

2. Method description 2.1. DEM model As aforementioned, both experimental and numerical techniques have been employed in obtaining the packing structures. However, the numerical methods are comparably more cost-effective, more controllable and can provide richer information. There are different numerical methods to simulate the packing of particles. A large number of early simulations were based on geometrical models which consider ideal hard particles with simplified interactions, and the packings were normally generated by kinetic methods such as Monte Carlo method. Later the DEM-based methods which consider more realistic granular particles have been adopted (Delaney et al., 2011; Dong et al., 2015; Kodam et al., 2010; Zhou et al., 2011). DEM has shown its advantages in more realistically representing the dynamic packing process and generating more comparable packing structures to experiments (Yi et al., 2011), owing to its usage of first principles with minimal assumptions (Zhu et al., 2007, 2008). In DEM, as originally proposed by Cundall and Strack (Cundall and Strack, 1979), two types of motions, namely, translational and rational motions, are considered for each particle, which are governed by Newton's second law of motion, given by:

332

mi

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

d vi = dt

∑ ( Fijn + Ftij ) + mi g

(1)

j

and,

Ii

d ωi = dt

∑ ( Tij + Tr, ij)

(2)

j

where mi, vi and ωi are the mass and the translational and angular velocities of particle i, respectively. Fijn and Fijt are the normal and tangential components of the contact force exerted on particles i by particle j respectively. g is the gravitational acceleration. Ii is the moment of inertia of particle i, which is a tensor but can be simplified by using body-fixed coordinate system. Tij and Tr,ij are the torques on particle i from particle j resulting from the total contact force (the sum of the normal and tangential forces) and rolling friction respectively. The most important aspect in applying these equations is to calculate the interaction between two particles, including the normal and tangential forces as well as the torques. This has been well established by the analytical method for spherical particles (Johnson, 1985; Langston et al., 1995), but is still a non-trivial problem for non-spherical particles (Dong et al., 2015; Lu et al., 2015). Recently we developed a novel method to tackle this issue in a general scheme. The method is based on orientation discretization and using the database rather than a full analytical way to obtain the overlap information, which hence is called as ODDS (orientation discretization database solution). The equations for calculating these forces and torques are listed in Table 1. The details of this method and its validation can be found in our previous work (Dong et al., 2015). Here we apply this method to generate the packings of identical axisymmetric ellipsoidal particles and cylindrical particles. The considered particle shapes are listed in Table 2. Note that the aspect ratio α, as defined in Table 2, is used as the shape parameter for both types of particles. For ellipsoidal particles, oblates are with α o1.0 and prolates α 41.0, while α ¼ 1.0 corresponds to spheres. Table 1 List of equations for contact force calculation. Force or torque

Equation

Normal elastic force, F cn ij

2 Y 3 1 −σ˜ 2

Normal damping force, F dn ij

−γn (

^ R¯ ξn3/2n

4Y 1 − σ˜ 2

(

)

^ n ^ mij R¯ ξn )1/2 vij⋅n

Tangential elastic force, F ct ij

⎡ −μ F cn ij ⎢ ⎣ 1 − 1 − min ( ξt , ξt , max )/ξt , max

Tangential damp-

⌢⌢ −γt (6μmij F cn (1 − min ( ξt , ξt , max )t )1/2/ξt , max )1/2 (vt , ij⋅ t ) t ij

ing force, F dt ij

(

Total normal force, Fnij

F cn + F dn ij ij

Total tangential force, Ftij

F ct + F dt ij ij

Torque by contact forces, Tij Torque by rolling friction, T r,ij

rij × (Fnij + Ftij )

)3/2 ⎤⎦⎥⌢t

−μr Ri |F cn |⌢ ωi ij

where Ri and Rj are the radii of particles i and j respectively, the radius is half of the particle size D; rij is the vector pointing from the center of particle i to the contact point with particle j; Y is the Young's modulus; σ˜ is the Poisson's ratio; γn and γt are the normal and tangential damping coefficients respectively; μ and μr are the sliding and rolling friction coefficients respectively; ξn is the overlap magnitude given as ξn = Ri + Rj − Rij for spheres, and is determined by contact overlap volume (Vo) for non-spherical particles ξn = Vo/(πR ) ; ξt is the total tangential displaceand ment, in each time step it is added with vt , ij⋅Δt , ⌢ ξt , max = μs ⎡⎣ ( 2 − σ˜ )/2 (1 − σ˜ ) ⎤⎦ ξn ; ⌢ n and t are the normal vector of the contact plane and the unit vector along the tangential direction respectively; and ⌢ ⌢ ω = ω / |ω | . vij = vj − vi + ωj × ri − ωj × r ji , vt , ij = (vij × n ) × n , ⌢ i i i

For cylindrical particles, disks are with α o1.0 and rods α 4 ¼ 1.0. The size of a particle D is defined as the diameter of its circumscribed sphere. A packing is formed by continuously dropping a batch of 100 particles at the top layer in a simulated zone of 10D  50D  10D in dimensions. Periodic boundary conditions are applied both in the X-axis and Z-axis directions and a plane boundary is placed at Y¼0. Gravity is along the – Y direction. After the feeding is stopped the particles will gradually become stable due to the energy dissipation between particles and between particles and walls. The maximal translational velocity and angular velocity of the particles both reach to nearly zero after the packing is stable. The final packings will be about 15D high, while the particle number will be different for particles of different shapes as their volumes are different. Such a packing protocol is similar to the normal pour packing method which has been widely used in the packing simulations (An et al., 2016; Dong et al., 2009a; Zhou et al., 2011). The packing process and a few final packings are demonstrated in Fig. 1. The simulation parameters are given in Table 3. In our previous studies we have also evidenced that the packing protocol and the numerical method can generate consistent and uniform packing structures (An et al., 2016; Wang et al., 2015; Yi et al., 2015). Compared to the literatures (An et al., 2016; Delaney and Cleary, 2010; Donev et al., 2004; Schaller et al., 2015; Torquato and Jiao, 2009; Xia et al., 2014), the system size and simulation settings including the periodic boundary conditions adopted in this work should be able to minimize the scale effect. For example, the particle number used in the simulations is similar or higher than those used in the experimental studies (Donev et al., 2004; Schaller et al., 2015; Xia et al., 2014). The obtained packing fractions have been shown to agree with previous numerical and experimental results as detailed in our previous study (Dong et al., 2015), while the local packing fraction distribution is also in accordance with the experimental data which will be shown in Section 3. These consistencies ensure the reliability of the simulation data. Also to avoid the boundary effect (Zou and Yu, 1996b) only the particles in the center region of the final packing will be used in the analysis, while those in the top and bottom layers of 2D thick are discarded, as it can be seen from Fig. 2, the packing fraction will have significant changes if these parts are included. Similar treatments have been applied in analyzing the packing structures obtained by numerical simulations or experiments (An et al., 2016; Dong et al., 2015; Schaller et al., 2015). 2.2. Voronoi cell construction method Originally Voronoi analysis is applied to volumeless points. Given a number of seed points in a space, for each seed point the region consisting of all points closer to it than to any other seed points is defined as its Voronoi cell. For the assembly of particles, the Voronoi cell of a particle contains the points with distance to its surface no greater than those to other particles. Such cells can be obtained by the intersecting of the bisect planes between a particle and its neighbors for the packings of monosized spheres, which are always polyhedrons. For the packings of multisized spheres, radical tessellations rather than original Voronoi cells need to be used to assure the obtained cells are still polyhedrons. Voronoi or radical tessellations for sphere packings can now be generated by open-source programs, e.g., see recent work done by Voroþ þ (Rycroft, 2009; Yi et al., 2012). For non-spherical particles, however, the Voronoi cells could have curved surfaces (Luchnikov et al., 1999; Medvedev et al., 2006). As the distance of a point to a non-spherical particle surface is often much more complicated than that of a spherical particle, there is no simple method to obtain the Voronoi cells for general non-spherical particles. Several advanced algorithms have been developed,

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

333

Table 2 List of 3D non-spherical particles simulated.

Shape

In body-fixed coordinate system

Aspect Ratio α=b/a

Ellipsoid b

(Oblate/Prolate)

a

b a

(a > b) Cylinder (Disk/Rod)

(a < b) α=h/d

d d

h h

(d>h)

(d<=h)

Fig. 1. Simulated packings of non-spherical particles: (a) a snapshot in the packing process of cylinders with α ¼1.0; (b) final packing state for disks with α ¼0.5; and, (c) final packing state for prolates with α ¼ 2.0.

including the Voronoi channel method based on tracking the imaginary empty sphere of variable size inside a packing (Luchnikov et al., 1999), and the set Voronoi diagram method based on a triangulation of the particles’ bounding surfaces (Schaller et al.,

2013). A more simple and universal numerical method has been proposed in the X-ray tomography (Al-Raoush and Alshibli, 2006), which is more computationally demanding than the above algorithms but relatively simpler to be generalized by computer

334

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

programs. Such a numerical method is adopted here and briefly described in below. The method is based on space discretization. For narrative convenience, we describe it in 2D in Fig. 3. The whole packing space is firstly discretized into finite pixels. The pixels inside the particles (with its center inside) are labeled as ‘0’ and other pixels as ‘1’, as Table 3 List of parameters used in the DEM simulations. Parameter

Value

Young's modulus, Y Poisson's ratio, s Normal damping coefficient , γn Tangential damping coefficient , γt Sliding friction coefficient, μs Rolling friction coefficient, μr Particle size (circumscribed sphere diameter) , D Particle density , ρp Time step, Δt Particle number, N Aspect ratio, α

1.0  107 Pa 0.29 0.3 0.3 0.3 0.005 1 mm 2.5  103 kg/m3 1.0  10  7 s 2000–12,000 0.3–2.5

Fig. 2. Packing fraction estimated by using different heights of the packings formed by different particles: Δ, prolates, α ¼1.33, μ ¼0.3; ■, oblates, α ¼0.3, μ¼0.3;  , disks, α¼ 0.3, μ¼ 0.3.

shown in Fig. 3(a). Then, the boundary pixels are identified by checking the label of a pixel and its surrounding pixels. If a pixel is labeled “0” and has one or more neighbor pixels labeled “ 1”, it will be marked as the boundary pixel. After that, we ‘burn’ those void pixels adjacent to the marked boundary pixels of all particles simultaneously, which yields a layer of pixels with 1 pixel distance from the surfaces of the particles (see pixels labeled with ‘1’ in Fig. 3 (b)). These “burned” pixels are then marked as the “new” boundary pixels and the last step will be re-iterated, which gives the second, the third, …, the nth layers of pixels consisting the space with nth pixels distance from the surfaces of the particles, as labeled in Fig. 3 (b). The iteration stops when there are no void pixels to be burnt. The common-boundary pixels are identified when the two or more ‘fires’ spreading from different particles meet, which are equal-distant to the surfaces of at least two particles, as shown in red color in Fig. 3 (b). The obtained common-boundary pixels constructs the boundaries of the Voronoi cells, and the pixels ‘burnt’ from a specific particle are closer to that particle than to any other ones, which construct the inside space of the Voronoi cell. Such a method has been successfully applied to the analysis of the packings of granular polygons (Wang et al., 2005). When applied to 3D, the finite pixels will be changed to voxels. This method is general to any convex shapes and the obtained Voronoi cells conform to the genuine definition, while the errors can be controlled by using enough fine pixels or voxels. Originally it has only been used to calculate the Voronoi cell volumes (AlRaoush and Alshibli, 2006), which can be obtained simply by counting the total number of voxels inside a Voronoi cell. However, as the Voronoi cells obtained in this way will have very rough boundaries as can be seen from Fig. 3(c), the surface areas cannot be correctly calculated. To overcome this problem, we use the Poisson surface reconstruction (Kazhdan et al., 2006) method to rebuild the surface from the face center points of the commonboundary voxels. The obtained surface of the Voronoi cell will then be smooth and the surface area can be obtained. Fig. 4 demonstrates an example of the procedures. The surface reconstruction and area calculation are conducted using the open-source program MeshLab (http://meshlab.sourceforge.net) for each particle, and then the Voronoi cells for the whole packing are obtained. In this work, we use the cubic voxels of size 1/128D in all three dimensions. From Fig. 3(c) we can see such a resolution can give enough details and the reconstructed surfaces are smooth. To quantitatively validate the method, we have compared the Voronoi cell obtained by this method and that by Voroþ þ for each

Fig. 3. Voronoi cell construction by space discretization: (a) Particle identification, pixels of particles are labeled ‘0’ and those of voids ‘  1’; (b) “Burning” voids from the boundary of particles layer by layer, the positive integers represent the number of the layers; (c) Obtained Voronoi cells for a packing of pentagons, where particles, pores, and boundaries of Voronoi cells are colored in blue, green and red, respectively.

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

335

Fig. 4. An example of Voronoi cell construction in a packing of oblates with α ¼0.75: (a) the points are the face centers of the boundary voxels of a Voronoi cell obtained from the space discretization method; (b) surface reconstructed by the points obtained in step (a); (c) the obtained Voronoi cell (yellow) enclosed by other neighbor Voronoi cells (gray) in the packing. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. The effect of the voxel size on the Voronoi cell properties: (a), Voronoi cell volume for each particle in the packing of oblates (α ¼ 0.3, μ¼ 0.3), values obtained by using voxel size 1/32D and 1/128D versus those by using voxel size 1/200D; (b), Average errors of the Voronoi cell volume and surface area as a function of the ratio of particle size to voxel size, for the packings of oblates (α ¼0.3, μ¼ 0.3) and disks (α¼ 0.3, μ¼ 0.3). The errors are calculated based on the values obtained with voxel size 1/200D.

particle in the packings of spherical particles. These packings are included in our simulations of the packings of ellipsoidal particles with α ¼ 1.0. The two methods give very close Voronoi cell volume and surface area for each particle. The differences are within 71.5% for the cell volume and 73% for the surface area. We also test the effect of voxel size on the results of two kinds of particles with the lowest aspect ratio, i.e., disk of α ¼0.3 and oblate of α ¼0.3. In Fig. 5(a), we compare the Voronoi cell volume obtained by using different voxel sizes for each particle in the packings of the oblates. It can be seen that using 1/32D voxels the obtained cell volumes deviate in a small degree from the values obtained by using 1/200D voxels, whereas using 1/128D voxels the values are almost the same as the latter. Fig. 5(b) further shows that if using the values obtained by 1/200D voxels as the base values, the errors for Voronoi cell volume and area decrease rapidly with the decrease of the voxel size, and the errors are generally smaller than 0.2% if 1/128D voxels are used. The validity of the method will be further confirmed in other aspects in our following analysis. We note that the accuracy of such method may be affected by the discretization resolution and the shape of the particle, while for the packings studied here the method gives enough accurate results.

3. Results and discussion 3.1. Voronoi cell morphology Fig. 6 demonstrates the Voronoi cells obtained in the packings of particles of different shapes. It can be seen that a Voronoi cell

generally encloses a particle with a similar flatness/slimness to that of the particle, i.e., when the aspect ratio of the particle deviates from 1.0, the Voronoi cell becomes flatter or slimmer. Such a Voronoi cell looks similar to the so-called interfacial volume surrounded by each particle (Xu et al., 2015, 2016), but here the cells of different particles cannot overlap, instead the two contacting Voronoi cells share a common surface seamlessly. Therefore the shape of the Voronoi cell is far more irregular and complicated than the enclosed particle. Particularly the Voronoi cell surfaces are made of various curved faces, which are very different from the polygon faces in the Voronoi or radical tessellations obtained in the sphere packings. Some properties studied in the sphere packings, such as those face related properties (Yang et al., 2002; Yi et al., 2012, 2015), may not be applicable. Therefore, in this work, our analysis is focused on the Voronoi cell properties as a starting point. 3.2. Voronoi cell volume We first discuss the reduced volume of the Voronoi cell, denoted as V, which is the Voronoi cell volume divided by the particle volume. Because the Voronoi cells completely divide the whole packing space, the average reduced Voronoi cell volume V can be linked to the overall packing fraction ρ for the packings of identical particles, given by V =1/ρ . Here we calculate ρ by counting the particles whose centers are inside the center region of the packing, and compare its reciprocal with V in Fig. 7. It can be seen that the two values are very close for each packing, which demonstrates the accuracy of our Voronoi cell construction method. Fig. 7 also demonstrates the effect of shape on V , which

336

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

Fig. 6. Examples of Voronoi cells obtained in the packings of non-spherical particles.

Fig. 7. Average reduced Voronoi cell volume as a function of particle aspect ratio for the packings of (a) ellipsoidal particles and (b) cylindrical particles with different sliding friction coefficients: Δ, μ¼ 0.1; ■, μ ¼0.3;  , μ¼ 0.6; ○, μ¼ 0.9. Lines are the reciprocals of the overall packing fractions.

are consistent with the previous findings on the changes of ρ . Specifically, for ellipsoidal particles, V shows minimal values at α E0.6 for oblates and at α E1.7 for prolates respectively, which are similar to the aspect ratios for the maximal packing fractions reported in the literatures (Man et al., 2005; Zhou et al., 2011). For cylindrical particles, V increases with either the increase or decrease of α from 1.0, which agrees well with the experimental observations (Zou and Yu, 1996a). Actually the dependence of ρ on α has been reported in our previous study (Dong et al., 2015) with μ ¼0.1 and μ ¼0.3. The new results here are the packings obtained

with μ ¼0.6 and μ ¼0.9, which show similar trends to those of μ ¼0.1 and μ ¼0.3 but just with higher V . The four series of data with different μ also show that generally the increase of μ leads to the increase of V , which further confirms that a packing is always more loose with a higher sliding friction coefficient. We then discuss the distribution of V. As can be seen from Fig. 6, Voronoi cells in a packing have very irregular morphologies, which represents the disordered structure of the packings. However, the statistical distributions of V for all the packings show a similar pattern as demonstrated in Fig. 8. It can be seen that V

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

337

Fig. 8. Distributions of the reduced Voronoi cell volumes for the packings of (a) ellipsoidal particles and (b) cylindrical particles with different aspect ratios and with μ¼0.3. Symbols:  , α¼ 0.3; ○, α¼ 1.0; ■, α ¼ 2.5. Lines are the best fittings of log-normal distributions.

always presents a single-peak distribution which is slightly rightskewed. Based on the feature, we find that the distribution of V can be described by the log-normal distribution function, given by:

f (V ) =

1 exp [ − (ln (V − 1) − ν )2 /2τ 2] τ 2π (V − 1)

(3)

where f (V ) is the probability density, v and τ are the geometrical parameters for the distribution, which are fitted for each packing. The fitting method is Generalized Reduced Gradient (GRG2) algorithm in Microsoft Excel Solver. Note (V 1) is the reduced volume of the void in a Voronoi cell, which is also called as free volume or excluded volume (Wang et al., 2015; Xia et al., 2014), and is no smaller than zero. Fig. 8 shows that for each packing, the fitted log-normal distribution can always describe the distribution of V well. The distribution of V or 1/V is an important topic in the statistical mechanics theory for particle packing. According to Edwards et al. (Edwards and Oakeshott, 1989; Song et al., 2008), the temperature in thermodynamics can be replaced by ρ for a granular matter, and the partition function is directly related to the distribution of local packing fraction which can be regarded as 1/V. There are also other kinds of distributions proposed for V or 1/V. For example, the k-gamma distribution is used for the packings of monosized spheres (Aste and Di Matteo, 2008) to evaluate the entropy change. For non-spherical particles, recently by X-ray tomography, Schaller et al. (2015) found that in a packing of oblates, the local packing fraction obtained as 1/V can be described by the Gaussian distribution. They also demonstrated that such distributions are similar for the packings with similar packing fractions regardless of the particle shape. We compare our results with their results in Fig. 9. Here the samples include not only the packings of oblates, but also of prolates and cylinders. It can be seen that the distributions in this work are close to the average distribution obtained by Schaller et al. (2015) at a similar packing fraction, while certain deviations can also be found which probably result from the differences in ρ as well as the particle shape. Although the Gaussian distribution can be used to model 1/V for certain packings, we find that in our studies, the log-normal distribution generally fits the distribution of V better. And this is also consistent with the previous findings in the packings of uniform cohesive spheres (Yang et al., 2002) and polygons (Wang et al., 2015), as well as with another experimental study on the packings of ellipsoidal particles (Xia et al., 2014). The form of the

Fig. 9. Distributions of local packing fractions (the reciprocal of the reduced Voronoi cell volume) for packings with ρ E0.625. Symbols: Δ, prolate, α ¼1.15, μ¼ 0.9; ■, oblate, α ¼0.3, μ¼0.6;  , rod, α ¼2.5, μ¼0.3; ○, sphere, μ¼ 0.3. Line is the averaged results obtained from X-ray tomography of the packings of oblates from Schaller et al. (2015).

distribution of V may be dependent on the packing system and requires further studies. For a log-normal distribution described by Eq. (3), its mean 2 value and standard deviation (denoted as s) are given by e ν + τ /2 þ1 and (e τ −1) e2ν + τ respectively. We use the fitted v and τ to calculate these two parameters for each packing. The calculated mean value of the distribution versus the average value (arithmetic mean) is demonstrated in Fig. 10. The two values are in agreement with less than 73% errors, which demonstrates the validity of the log-normal distribution of V. However, the mean value of the distribution is generally higher than the average value. This is because the ideal log-normal distribution has an infinite long tail with V-1, whereas actually V should have an upper bound due to the stability limit. Hence V of the distribution will be higher the average V. Such an overestimation has also been observed in the packings of 2D polygons (Wang et al., 2015). This indicates that we probably can use truncated distributions with the known bounds of V to more accurately model the distribution of V, which needs further studies. The effect of the aspect ratio on the standard deviation of the distribution of V is shown in Fig. 11. Generally for ellipsoidal 2

2

338

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

particles, s(V) decreases when α decreases or increases from 1.0, reaches two valleys at α E 0.6 and α E1.7 respectively, and finally increases with the further decrease of α for oblates or further increase of α for prolates. For cylindrical particles, the standard deviation has only one minimal value at α ¼1.0, while either the increase or decrease of α from 1.0 leads to the increase of s(V). The increase of the sliding friction coefficient always increases s(V). These changes are comparable to the changes of the average V shown in Fig. 8, indicating that the standard deviation of V may have a consistent change with V and hence packing fraction. We thus plot s(V) versus ρ in Fig. 12. Interestingly, it is shown that s(V) can be linked to ρ in one universal correlation for both cylindrical and ellipsoidal particles under different sliding friction coefficients. Quantitatively, the relationship can be fitted by an exponential function as shown in Fig. 12, given by:

σ (V ) = 6. 523 exp ( − 4. 516ρ)

(4)

For a long time in engineering applications, packing fraction (or its equivalent: porosity) has been used as a single parameter to model many properties of a packing (Cheng and Yu, 2013; Coelho

Fig. 10. The mean value of the log-normal distribution of V versus the value averaged over particles, including both ellipsoidal particles (○) and cylindrical particles (  ).

et al., 1997; Dullien, 1992 ; Rong et al., 2013), as such a macroscopic parameter is far more easy to measure than the microscopic structure. But this treatment implicitly assumes that the structure of a packing should be able to be correlated with the packing fraction. In early studies, such possible correlations have been tested for the packings of monosized coarse spherical particles in terms of CN and referred as “quasi-universality” (Jullien et al., 1996). Recently, such an assumption has been comprehensively examined in the packings of monosized spheres in a wide range of packing fraction, in which it has been confirmed that the statistical distributions of many structural parameters, including RDF, CN, and Voronoi cell properties, have almost a one-to-one correlation with ρ (An et al., 2016). Non-spherical particles make the packings more complicated, and such quasi-university correlations may not be fully held. Actually it is found that some of the parameters can be related to packing fraction for different shaped particles while others may not. For example, it is shown that the relationship between CN and ρ is different for prolates and oblates (Zhou et al., 2011) and for different polygons (Wang et al., 2015). Interestingly here we find that the relationship between the statistical distribution of the reduced Voronoi cell volume and packing fraction has such “quasi-universality”, for both ellipsoidal particles and cylindrical particles and/or under different sliding friction

Fig. 12. Standard deviation of the fitted log-normal distribution of V as a function of packing fraction, for all the packings studied: ○, ellipsoidal particles;  , cylindrical particles.

Fig. 11. Standard deviation of the distribution of the reduced Voronoi cell volume as a function of particle aspect ratio for the packings of (a) ellipsoidal particles and (b) cylindrical particles with different sliding friction coefficients: Δ, μ¼0.1; ■, μ¼ 0.3;  , μ ¼0.6; ○, μ¼ 0.9.

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

339

coefficients. Such a uniform correlation probably results from the distribution form of the reduced Voronoi cell volume. In the following analysis, we will also try to explore such kind of correlations for other Voronoi cell properties. Note in this work we focus on finding the general relationships between the local and global packing properties rather than directly fitting the dependencies of Voronoi cell properties on the particle shape and sliding friction coefficient. As the equations obtained by such fittings could be system dependent. On the contrary, as shown in the previous studies in the sphere packings, finding the general relationships between the local and global packing properties can lead to more general theoretical models to predict the effects of different variables on the packing structures (Dong et al., 2006, 2012), while these relationships are rather limited for non-spherical particles in the literatures.

obtained by fitting, we can also calculate the mean value and standard deviation. Fig. 15 shows S calculated from the distribution versus that averaged over particles. Similar to V , we can see that the two values are close but the former is generally higher than the latter. This should also be attributed to the long tail of the log-normal distribution of S, which cannot be reached by the real Voronoi cells. As pointed in the above discussion, this suggests that we can use the truncated distribution for more accurate modeling. The calculated standard deviation of the distribution of S is demonstrated in Fig. 16. Generally its changes are similar to those of S if compared to Fig. 13. Once we plot the mean value and standard deviation of S together in Fig. 17, we can clearly see that s(S) has a simple linear relationship with S for all the packings studied in this work, given by:

3.3. Voronoi cell surface area

σ ( S )=0. 5817 S −0. 5619

Fig. 13 plots the average reduced Voronoi cell surface area, denoted as S, as a function of particle aspect ratio. Note S is reduced by the surface area of the enclosed particle. It can be seen that for oblates and prolates, S generally decreases when the particle deviates from sphere, either becomes flatter for oblates (α decreases from 1.0) or slimmer for prolates (α increases from 1.0). Different from that of V , there is no critical change for oblates at α E0.6. For prolates at α E 1.7, where V reaches a minimum, S also has a minimal value yet such critical change is not as distinct as that of V , as S just slightly increase with the further rise of α. For cylindrical particles, S simply increases with the increase of α, and there is no critical change at α ¼1.0. In short, the change of S is different from V , showing that S may not be directly linked to ρ. However, the effect of friction is consistent: with the decrease of μ, S always decreases. This is probably because that the Voronoi cells are generally smaller when the packing is denser with a smaller μ. We also investigate the distribution of S. From Fig. 14 it can be seen that similar to V, S also presents a one-peak distribution with a slightly positive skewness for each packing. Thus we find that its distribution can also be described by the log-normal distribution function, given by:

f (S ) =

1 2 ⎡ ⎤ exp ⎣ − ln ( S − 1) − v /2τ 2⎦ τ 2π (S − 1)

(

)

(5)

Note similar to V, S should also be no smaller than 1 so that (S 1) no smaller than 0. With the parameters for each packing

(6)

On the other hand, neither S nor s(S) can be simply correlated with ρ, which is different from the packings of uniform cohesive spheres (Yang et al., 2002). 3.4. Voronoi cell sphericity With the volumes and surface areas of the Voronoi cells known, we can evaluate their morphologies by sphericity ψV , given by 1

ψV =

2

π 3 (6VV ) 3 . SV

Note here VV and SV are the volume and surface area of

a Voronoi cell respectively but not the reduced values. Yet we have: 1

ψV =

2

2

2

π 3 (6VP ) 3 V 3 V3 =ψP SP S S

(7)

where VP , SP and ψP are the volume, surface area and sphericity of the particle respectively. Fig. 18 plots the average sphericity ψV of the Voronoi cells for the packings of ellipsoidal and cylindrical particles with different sliding frictions coefficients. Interestingly it can be seen that μ has very little effect on ψV . For both ellipsoidal and cylindrical par-

ticles, ψV always decreases with α increasing or decreasing from 1.0, demonstrating that the sphericity of the Voronoi cell is related to the sphericity of the particle. Fig. 19 shows the averaged Voronoi cell sphericity versus the particle sphericity for all the packings studied in this work. Note because μ has little effect on ψV , here

ψV

is averaged over packings with different

μ for each particle

Fig. 13. Average reduced Voronoi cell surface area as a function of particle aspect ratio for the packings of (a) ellipsoidal particles and (b) cylindrical particles with different sliding friction coefficients: Δ, μ¼0.1; ■, μ¼ 0.3;  , μ¼0.6; ○, μ ¼0.9.

340

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

Fig. 14. Distribution of the reduced Voronoi cell surface area for the packings of (a) ellipsoidal particles and (b) cylindrical particles with different aspect ratios at μ¼ 0.3. Symbols:  , α¼ 0.3; ○, α¼ 1.0; ■, α ¼2.5. Lines are the best fittings of log-normal distributions.

shape. From Fig. 19, the mean Voronoi cell sphericity in a packing can be correlated with the particle sphericity for both ellipsoidal and cylindrical particles, given by:

ψV =0. 4962 ln ( ψP )+0. 9103

(8)

We also obtain the distributions of ψV for the studied packings and find they are generally very narrow, indicating that the Voronoi cells in a packing have similar sphericities. That can also explain why the Voronoi cell volume and surface area have similar log-normal distribution forms as they can be linked with ψV . With ψV similar for each particle in a packing, the average sphericity can be linked with the average reduced Voronoi cell volume and surface area, given by: 2

Fig. 15. The mean reduced Voronoi cell surface area calculated from the log-normal distribution versus that averaged over particles, for all the packings studied: ○, ellipsoidal particles;  , cylindrical particles.

ψV =ψP

V 3 S

(9)

Fig. 16. Standard deviation of the distribution of the reduced Voronoi cell surface area as a function of particle aspect ratio for the packings of (a) ellipsoidal particles and (b) cylindrical particles with different sliding friction coefficients: Δ, μ¼0.1; ■, μ¼ 0.3;  , μ ¼0.6; ○, μ¼ 0.9.

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

Combining Eqs. (8) and (9), S can be obtained as:

ψ S= P ψV

V

2 3=

2 ψP ρ− 3 0. 4962 ln ( ψP )+0. 9103

(10)

341

particles, which deserves further studies. In particular, these equations can be linked to the existing fundamental theories for sphere packings, such as the statistical mechanics theory (Aste and Di Matteo, 2008; Song et al., 2008) and the granocentric packing

Fig. 20 demonstrates that S predicted by Eq. (10) has a good agreement with the values measured in the simulated packings. By predicting S , we can also obtain the standard deviation of S from Eq. (6). 3.5. Implication of the equations We note that by establishing the Eqs. (4), (6), (8) and (10), we actually can predict the Voronoi cell properties of a packing of identical ellipsoidal or cylindrical particles based on the particle sphericity and overall packing fraction (or porosity). Such universal correlations are based on the distribution forms of the Voronoi cell volume and surface area, and the link between the sphericity of the Voronoi cell and the enclosed particle, which is probably due to the similar morphologies for the Voronoi cell and the interfacial volume (Xu et al., 2015, 2016). The inherent relationships between different Voronoi cell properties are also implied in the equations. These equations indicate that there probably exist some general laws for the local structures of the packings of non-spherical

Fig. 19. Average Voronoi cell sphericity as a function of particle sphericity: ○, ellipsoidal particles;  , cylindrical particles.

Fig. 17. Standard deviation of the fitted log-normal distribution of the reduced Voronoi cell surface area as a function of its average value, for all the packings studied: ○, ellipsoidal particles;  , cylindrical particles.

Fig. 20. Predicted average reduced Voronoi cell surface area versus the value obtained in the packings, including both ellipsoidal particles (○) and cylindrical particles (  ).

Fig. 18. Average Voronoi cell sphericity as a function of particle aspect ratio for the packings of (a) ellipsoidal particles and (b) cylindrical particles with different sliding friction coefficients: Δ, μ¼ 0.1; ■, μ¼0.3;  , μ¼ 0.6; ○, μ¼ 0.9.

342

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

model (Clusel et al., 2009), and help the extension and generalization of these theories for non-spherical particles, such as helping construct the partition function of local packing fraction in the statistical mechanics theory. On the other hand, these equations can also be used for evaluating the transport properties of packings based on the Voronoi cells (Cheng and Yu, 2013; Dullien, 1992 ). For example, in modeling the heat conductivity of a packing, the Voronoi cell properties are required to calculate the conduction and radiation between the surfaces of two particles (Cheng and Yu, 2013). Another example is that in determining the distribution of fluid-particle interactions in the fluid-particle multiphase flow systems (Rong et al., 2013, 2015), the local packing fraction can be used. These potential applications of the current findings will be pursued in our future studies.

4. Conclusion We have conducted the Voronoi analysis on the packings of identical axisymmetric ellipsoidal particles and cylindrical particles with different aspect ratios (α ¼0.3–2.5) and sliding friction coefficients (μ ¼ 0.1–0.9). The Voronoi cells are constructed by the combined space discretization and mesh reconstruction method, which has shown to be universal for the packings of non-spherical particles. The effects of α and μ on the Voronoi cells have been quantified in terms of the reduced volume (V), reduced surface area (S) and sphericity ( ψV ), including their average values and distributions. The changes of the average properties can be summarized below:


The effect of α on V are just opposite to that of packing fraction. For oblates, V decreases with the decrease of α from






1.0 until about 0.6, and then increases; and for prolates, V decreases with the increase of α from 1.0 until about 1.7, and then increases. For cylindrical particles, V increases with α increasing or decreasing from 1.0. For oblates and prolates, S mainly decreases with α deviating from 1.0. While for cylindrical particles, S increases monotonously with the increase of α. For both ellipsoidal and cylindrical particles, ψV decreases with α deviating from 1.0, and it can be correlated to the particle sphericity by a general equation. The increase of μ always results in a loose packing and consequently larger V and S , whereas μ has little effect on ψV .

It is found that both V and S can be described by the log-normal distributions, although the extremely high V and S in the long tails of the distributions do not exist. In addition, the mean values and standard deviations of the distributions of these properties can be universally correlated with the particle sphericity ψP and overall packing fraction ρ. Specifically, for the packings of cylindrical and ellipsoidal particles with different frictions studied in this work, the distributions of the Voronoi cell properties can be predicted by a group of equations given in below:

⎧ V =1/ρ , σ ( V )=6. 523 exp ( −4. 516ρ); ⎪ ψP 2 ⎪ ρ− 3 , σ ( S )=0. 5817 S −0. 5619; ⎨ S= ( ψ )+ 0. 4962 ln 0. 9103 P ⎪ ⎪ ψV =0. 4962 ln ( ψP )+0. 9103,σ ( ψV )≈0. ⎩ These results have improved our understanding on the packing structures of non-spherical particles. In particular, the quantitative relationships provide a useful basis for modeling the transport properties of non-sphere packings, and are also helpful to the theoretical modeling of granular matter by the statistical

mechanics theory. Such relationships are expected to be further analyzed and improved by considering the bounds of the properties and probably extended to the particles of other shapes and to the mixtures of non-spherical particles in future studies.

Acknowledgment The authors are grateful for the financial support from Australian Research Council (DE120100960).

References Al-Raoush, R., Alshibli, K.A., 2006. Distribution of local void ratio in porous media systems from 3D X-ray microtomography images. Phys. A: Stat. Mech. Appl. 361, 441–456. An, X.Z., Dong, K.J., Yang, R.Y., Zou, R.P., Wang, C.C., Yu, A.B., 2016. Quasi-universality in the packing of uniform spheres under gravity. Granul. Matter 18, 1–7. Anikeenko, A.V., Medvedev, N.N., 2007. Polytetrahedral nature of the dense disordered packings of hard spheres. Phys. Rev. Lett. 98, 235504. Aste, T., 2006. Volume fluctuations and geometrical constraints in granular packs. Phys. Rev. Lett. 96, 018002. Aste, T., Di Matteo, T., 2008. Emergence of Gamma distributions in granular materials and packing models. Phys. Rev. E 77, 021309. Baule, A., Mari, R., Bo, L., Portal, L., Makse, H.A., 2013. Mean-field theory of random close packings of axisymmetric particles. Nat. Commun., 4. Bideau, D., Hansen, A., 1993. Disorder and Granular Media. Elsevier Science, North Holland, Netherlands. Cheng, G.J., Yu, A.B., 2013. Particle scale evaluation of the effective thermal conductivity from the structure of a packed bed: radiation heat transfer. Ind. Eng. Chem. Res. 52, 12202–12211. Clusel, M., Corwin, E.I., Siemens, A.O.N., Brujic, J., 2009. A /'granocentric/' model for random packing of jammed emulsions. Nature 460, 611–615. Coelho, D., Thovert, J.F., Adler, P.M., 1997. Geometrical and transport properties of random packings of spheres and aspherical particles. Phys. Rev. E 55, 1959–1978. Cundall, P.A., Strack, O.D., 1979. A discrete numerical model for granular assemblies. Geotechnique 29, 47–65. Delaney, G.W., Cleary, P.W., 2010. The packing properties of super ellipsoids. Eur. Phys. Lett. 89, 34002. Delaney, G.W., Hilton, J.E., Cleary, P.W., 2011. Defining random loose packing for nonspherical grains. Phys. Rev. E 83, 051305. Donev, A., Cisse, I., Sachs, D., Variano, E.A., Stillinger, F.H., Connelly, R., Torquato, S., Chaikin, P.M., 2004. Improving the density of jammed disordered packings using ellipsoids. Science 303, 990–993. Dong, K.J., Wang, C.C., Yu, A.B., 2015. A novel method based on orientation discretization for discrete element modeling of non-spherical particles. Chem. Eng. Sci. 126, 500–516. Dong, K.J., Yang, R.Y., Zou, R.P., Yu, A.B., 2006. Role of interparticle forces in the formation of random loose packing. Phys. Rev. Lett. 96, 145505. Dong, K.J., Yang, R.Y., Zou, R.P., An, X.Z., Yu, A.B., 2009a. Critical states and phase diagram in the packing of uniform spheres. Eur. Phys. Lett. 86, 6. Dong, K.J., Zou, R.P., Yang, R.Y., Yu, A.B., Roach, G., 2009b. DEM simulation of cake formation in sedimentation and filtration. Miner. Eng. 22, 921–930. Dong, K.J., Yang, R.Y., Zou, R.P., Yu, A.B., 2012. Settling of particles in liquids: Effects of material properties. AIChE Journal 58, 1409–1421. Dullien, F.A., 1992. Porous Media: Fluid Transport and Pore Structure, 2nd ed. Academic press, San Diego, California, USA. Džiugys, A., Peters, B., 2001. An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers. Granul. Matter 3, 231–266. Edwards, S.F., Oakeshott, R.B.S., 1989. Theory of Powders. Phys. A – Stat. Mech. Appl. 157, 1080–1090. Favier, J., Abbaspour-Fard, M., Kremmer, M., Raji, A., 1999. Shape representation of axi-symmetrical, non-spherical particles in discrete element simulation using multi-element model particles. Eng. Comput. 16, 467–480. German, R.M., 1989. Particle Packing Characteristics. Metal Powder Industries Federations, Princeton, New Jersey. Glotzer, S.C., Engel, M., 2011. Materials science: Complex order in soft matter. Nature 471, 309–310. Guo, B.-Y., Yang, S.-Y., Xing, M., Dong, K.-J., Yu, A.-B., Guo, J., 2013. Toward the development of an integrated multiscale model for electrostatic precipitation. Ind. Eng. Chem. Res. 52, 11282–11293. Guo, Y., Wassgren, C., Hancock, B., Ketterhagen, W., Curtis, J., 2015. Computational study of granular shear flows of dry flexible fibres using the discrete element method. J. Fluid Mech. 775, 24–52. Johnson, K., 1985. Contact Mechanics. Cambridge University Press, Cambridge. Jullien, R., Jund, P., Caprion, D., Quitmann, D., 1996. Computer investigation of longrange correlations and local order in random packings of spheres. Phys. Rev. E 54, 6035–6041.

K. Dong et al. / Chemical Engineering Science 153 (2016) 330–343

Kazhdan, M., Bolitho, M., Hoppe, H., 2006. Poisson surface reconstruction, In: Proceedings of the fourth Eurographics symposium on Geometry processing. Kodam, M., Bharadwaj, R., Curtis, J., Hancock, B., Wassgren, C., 2010. Cylindrical object contact detection for use in discrete element method simulations, Part II – experimental validation. Chem. Eng. Sci. 65, 5863–5871. Langston, P.A., Tüzün, U., Heyes, D.M., 1995. Discrete element simulation of granular flow in 2D and 3D hoppers: dependence of discharge rate and wall stress on particle interactions. Chem. Eng. Sci. 50, 967–987. Latham, J., Munjiza, A., 2004. The modelling of particle systems with real shapes. Philos. Trans. – R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 362, 1953–1972. Li, S., Lu, P., Jin, W., Meng, L., 2013. Quasi-random packing of tetrahedra. Soft Matter 9, 9298–9302. Lu, G., Third, J.R., Müller, C.R., 2015. Discrete element models for non-spherical particle systems: From theoretical developments to applications. Chem. Eng. Sci. 127, 425–465. Luchnikov, V.A., Medvedev, N.N., Oger, L., Troadec, J.P., 1999. Voronoi–Delaunay analysis of voids in systems of nonspherical particles. Phys. Rev. E 59, 7205–7212. Man, W., Donev, A., Stillinger, F.H., Sullivan, M.T., Russel, W.B., Heeger, D., Inati, S., Torquato, S., Chaikin, P.M., 2005. Experiments on random packings of ellipsoids. Phys. Rev. Lett. 94, 198001. Medvedev, N.N., Voloshin, V.P., Luchnikov, V.A., Gavrilova, M.L., 2006. An algorithm for three-dimensional Voronoi S-network. J. Comput. Chem. 27, 1676–1692. Moreno-Atanasio, R., Williams, R.A., Jia, X., 2010. Combining X-ray microtomography with computer simulation for analysis of granular and porous materials. Particuology 8, 81–99. Rong, L.W., Dong, K.J., Yu, A.B., 2013. Lattice–Boltzmann simulation of fluid flow through packed beds of uniform spheres: effect of porosity. Chem. Eng. Sci. 99, 44–58. Rong, L.W., Zhou, Z.Y., Yu, A.B., 2015. Lattice–Boltzmann simulation of fluid flow through packed beds of uniform ellipsoids. Powder Technol. 285, 146–156. Rycroft, C., 2009. Voroþ þ: A Three-dimensional Voronoi Cell Library in Cþ þ . Chaos 19. Schaller, F.M., Neudecker, M., Saadatfar, M., Delaney, G.W., Schröder-Turk, G.E., Schröter, M., 2015. Local origin of global contact numbers in frictional ellipsoid packings. Phys. Rev. Lett. 114, 158001. Schaller, F.M., Kapfer, S.C., Evans, M.E., Hoffmann, M.J., Aste, T., Saadatfar, M., Mecke, K., Delaney, G.W., Schröder-Turk, G.E., 2013. Set Voronoi diagrams of 3D assemblies of aspherical particles. Philos. Mag. 93, 3993–4017. Song, C., Wang, P., Makse, H.A., 2008. A phase diagram for jammed matter. Nature 453, 629–632. Tian, Z.A., Dong, K.J., Yu, A.B., 2014. Structural evolution in the packing of uniform spheres. Phys. Rev. E 89, 032202. Tian, Z.A., Liu, R.S., Dong, K.J., Yu, A.B., 2011. A new method for analyzing the local structures of disordered systems. Eur. Phys. Lett. 96, 6. Torquato, S., Jiao, Y., 2009. Dense packings of the Platonic and Archimedean solids. Nature 460, 876–879.

343

Vorobiev, O., 2012. Simple common plane contact algorithm. Int. J. Numer. Methods Eng. 90, 243–268. Wachs, A., Girolami, L., Vinay, G., Ferrer, G., 2012. Grains3D, a flexible DEM approach for particles of arbitrary convex shape – Part I: numerical model and validations. Powder Technol. 224, 374–389. Wang, C., Dong, K., Yu, A., 2015. Structural characterization of the packings of granular regular polygons. Phys. Rev. E 92, 062203. Xia, C., Zhu, K., Cao, Y., Sun, H., Kou, B., Wang, Y., 2014. X-ray tomography study of the random packing structure of ellipsoids. Soft Matter 10, 990–996. Xiang, J., Munjiza, A., Latham, J.P., Guises, R., 2009. On the validation of DEM and FEM/DEM models in 2D and 3D. Eng. Comput. 26, 673–687. Xu, W., Wang, H., Niu, Y., Bai, J., 2016. Insight into interfacial effect on effective physical properties of fibrous materials. I. The volume fraction of soft interfaces around anisotropic fibers. J. Chem. Phys. 144, 014703. Xu, W., Duan, Q., Ma, H., Chen, W., Chen, H., 2015. Interfacial effect on physical properties of composite media: interfacial volume fraction with non-spherical hard-core-soft-shell-structured particles. Scientific Reports 5, 16003. Yang, R.Y., Zou, R.P., Yu, A.B., 2002. Voronoi tessellation of the packing of fine uniform spheres. Phys. Rev. E 65, 041302. Yang, S., Dong, K., Zou, R., Yu, A., Guo, J., 2013. Packing of fine particles in an electrical field. Granul. Matter, 1–10. Yi, L.Y., Dong, K.J., Zou, R.P., Yu, A.B., 2011. Coordination number of the packing of ternary mixtures of spheres: DEM simulations versus measurements. Ind. Eng. Chem. Res. 50, 8773–8785. Yi, L.Y., Dong, K.J., Zou, R.P., Yu, A.B., 2012. Radical tessellation of the packing of ternary mixtures of spheres. Powder Technol. 224, 129–137. Yi, L.Y., Dong, K.J., Zou, R.P., Yu, A.B., 2015. Radical tessellation of the packing of spheres with a log-normal size distribution. Phys. Rev. E 92, 032201. Zhao, J., Li, S., Zou, R., Yu, A., 2012. Dense random packings of spherocylinders. Soft Matter 8, 1003–1009. Zhou, Z.Y., Zou, R.P., Pinson, D., Yu, A.B., 2011. Dynamic Simulation of the packing of ellipsoidal particles. Ind. Eng. Chem. Res. 50, 9787–9798. Zhu, H.P., Zhou, Z.Y., Yang, R.Y., Yu, A.B., 2007. Discrete particle simulation of particulate systems: theoretical developments. Chem. Eng. Sci. 62, 3378–3396. Zhu, H.P., Zhou, Z.Y., Yang, R.Y., Yu, A.B., 2008. Discrete particle simulation of particulate systems: a review of major applications and findings. Chem. Eng. Sci. 63, 5728–5770. Zou, R.P., Yu, A.B., 1996a. Evaluation of the packing characteristics of mono-sized non-spherical particles. Powder Technol. 88, 71–79. Zou, R.P., Yu, A.B., 1996b. Wall effect on the packing of cylindrical particles. Chem. Eng. Sci. 51, 1177–1180.