Periodic Boundary Conditions for Discrete Element Method Simulation

Article pubs.acs.org/IECR

Periodic Boundary Conditions for Discrete Element Method Simulation of Particle Flow in Cylindrical Vessels Wenjing Yang,† Zongyan Zhou,† David Pinson,‡ and Aibing Yu*,† †

Laboratory for Simulation and Modelling of Particulate Systems, School of Materials Science and Engineering, The University of New South Wales, Sydney, NSW 2052, Australia ‡ Steelmaking Technology and Planning, BlueScope Steel, P.O. Box 202, Port Kembla, NSW 2505, Australia ABSTRACT: The application of the discrete element method (DEM) for industrial process simulation is often limited to twodimensional slot models for the purpose of reducing the number of particles and hence the computational effort. In this work, a three-dimensional (3D) sector model is proposed with circumferential periodic boundary conditions (PBC) to represent a full 3D model for cylindrical geometries. Two typical processes, i.e. hopper flow and piling, are used to test the sector model. This reveals that the simulated flow features in the 3D sector model with suitable PBC are consistent with those in the full 3D model qualitatively and quantitatively. The proposed sector model can be generally used to investigate granular flow in a full 3D cylindrical geometry in the DEM simulation.

1. INTRODUCTION

conditions, e.g. geometries with particles having small radial movements, or those sector domains with a large sector angle. In slot or sector models, front and rear walls have a significant effect on granular flow as demonstrated in blast furnace modeling.9,10 In order to eliminate or reduce the wall effect on solid flow, the boundary has to be treated properly. There are several options to set the boundary, such as “rigid boundary”,13,14 “frictionless boundary”,15 and “periodic boundary”.9,11,16−19 To date, periodic boundary conditions (PBC) have been recognized as the most useful treatment to eliminate the wall effect and, at the same time, can maintain a continuous internal flow involving particle−particle contacts throughout the computational domain. PBC for slot models have been well established, and the attempts to apply PBC to sector models have been made but not fully established yet. Generally speaking, to apply PBC to sector models, here referred to as “circumferential PBC”, the small region where two periodic boundaries converged has to be treated with caution; otherwise, particles may flow out of the computational domain. Due to the limited space adjacent to the center, the circumferential PBC cannot fully function at the center. Some attempts have been made to solve this issue. For example, the center line can be set as a rigid wall,13 but the contact force from the wall improperly affects the particle horizontal movement and subsequently influences the whole flow structure. As demonstrated by Cui et al.,11 the friction coefficient is difficult to determine. Another method is to remove the central region so that the problem in computation is also resolved. However, granular flow features in the center are missing, which may be the most concerning issue in some devices.6,12 On the other hand, Cui et al.11 assumed that once the particle center was located at the sector center, it could not move horizontally. Such a treatment may not affect their force test,

Particulate systems are quite common in industry, and their behavior can be complex due to the interactions between individual particles and their interactions with surrounding gas, liquid, and walls. To describe the flow of granular materials, the discrete element method (DEM) originally developed by Cundall and Strack1 is a commonly applied approach. The method considers a finite number of discrete particles interacting with each other and the translational and rotational motions of every particle in the system on the basis of Newton’s second law of motion. DEM simulations can generate rich dynamic and microscopic information, such as the trajectories of and transient forces acted on individual particles, which is important to fundamental understanding but difficult to obtain by physical experimentation. With the rapid development in computer technology and the improvements in simulation techniques, DEM has been increasingly used in the study of particulate systems as reviewed by Zhu et al.2,3 Practical particulate systems often involve millions of particles, which could cause a difficulty in computation with the current computational capability. To overcome this difficulty, various methods have been proposed to reduce the number of particles and hence the computing cost. For example, small scale models are used to cut down the particle number.4−6 In particular, the socalled two-dimensional (2D) slot model has been extensively used in the DEM simulation, not only to represent cuboidal geometries,7,8 but also to replace cylindrical geometries, due to its simple mathematical treatments.6,9 These studies have demonstrated that slot models can provide useful information for better understanding of granular flow, especially in cuboidal vessels. However, for cylindrical geometries, slot models may neglect some important flow features and hence cannot exactly describe particle flow behavior in 3D cylindrical vessels.6,10 Theoretically speaking, sector models are preferred to represent 3D cylindrical vessels with proper boundary conditions.11,12 Such a sector model may also be useful for other certain © 2014 American Chemical Society

Received: Revised: Accepted: Published: 8245

December 8, 2013 March 29, 2014 April 22, 2014 April 23, 2014 dx.doi.org/10.1021/ie404158e | Ind. Eng. Chem. Res. 2014, 53, 8245−8256

Industrial & Engineering Chemistry Research

Article kc

as particles are almost motionless during the compression. However, severe effects may be generated for the cases with particles having horizontal movements, leading to unreasonable results. Therefore, the purpose of this paper is to present a new method to establish better circumferential PBC for sector models. The treatments for the PBC are described first. Then its applicability is verified by two case studies through direct comparisons with the corresponding full 3D models. The effect of the computational domain, represented by the angle of a sector model, is also examined.

mi dvi/dt =

∑ (fc,ij + fd, ij) + mi g (1)

j=1

and kc

Ii dωi /dt =

∑ Mij (2)

j=1

where vi and ωi are the translational and angular velocities of the particle, respectively, and kc is the number of particles in interaction with the particle. The forces involved are the gravitational force mig, and interparticle forces which include elastic force fc,ij and viscous damping force fd,ij. The torques acting on particle i by particle j include the following: Mt,ij, generated by the tangential force, and Mr,ij, commonly known as the rolling friction torque. The equations to determine these forces and torques in the present study have been well established in the literature and particularly consistent with those used by Zhou et al.20 2.2. Periodic Boundary Condition. Granular flow behavior near boundaries is affected by the settings of the boundary conditions. Three types of boundary conditions are commonly used in DEM simulation, including walls, inflow/outflow, and PBC. In particular, PBC are significant in reducing the number of

2. NUMERICAL TREATMENTS 2.1. DEM Model. According to the DEM, the governing equations for the translational and rotational motion of particle i with radius Ri, mass mi, and moment of inertia Ii can be written as

Figure 1. Illustration of the rectangle periodic boundary conditions.

Figure 4. Geometry of the cylindrical hopper used in the DEM simulation, where θ is the sector angle.

Figure 2. (a) Geometry and the boundary location of a sector model and (b) illustration of the circumferential periodic boundary conditions.

Table 1. Parameters used in the DEM Simulation for Hopper Flow variables Young’s modulus Poisson ratio sliding friction rolling friction time step particle density particle diameter models particle number running timea

value 1.0 × 107 Pa 0.3 particle−particle 0.6; particle−wall 0.5 0.02d 2.26 × 10−5 s 3800 kg/m3 0.005 m 3D θ = 90° θ = 60° θ = 45° 120 000 36 000 24 000 18 000 90 h 10 h 10 h 5h

θ = 30° 12 000 5h

a

All the cases presented in this work are calculated on the National Computational Infrastructure (NCI), a Sun Constellation Cluster with 1492 nodes, each containing 2 quad core Nehalem processors.

Figure 3. Illustration of the proposed treatments near the central region at the sector center. 8246

dx.doi.org/10.1021/ie404158e | Ind. Eng. Chem. Res. 2014, 53, 8245−8256


Article

Figure 7. Spatial distribution of coordination number at t = 0.7 s (only left side shown).

Figure 8. Probability density distribution of coordination number at t = 0.7 s.

Figure 5. Snapshots of discharging process in a cylindrical hopper: (a) full 3D, (b) 90° sector, (c) 60° sector, (d) 45° sector, and (e) 30° sector.

Figure 9. Correlation between the arc length and CN.

particles in the DEM simulation. Parallel PBC for 2D slot models have been successfully applied.18,19 The principle involved is simple, as indicated in Figure 1. When a particle moves out of one periodic boundary, it can be treated as moving into the other periodic boundary. Here, particles 1 and 1′ (or particles 3 and 3′) have the same properties except for the location difference in the y direction. Specifically, the theoretical treatments mainly include two aspects: (i) Force treatment. When the particle surface, e.g. particle 1 in Figure 1, exceeds PB1, the interactions with particles near PB2, e.g. particle 2, should be taken into account. The force calculation is based on the interaction Fn′ between

Figure 6. Snapshots of discharging process in a cylindrical hopper on the horizontal planes at the height of 0.1 m. 8247



Article

where α is the rotation angle. α = −θ for clockwise rotation, and α = θ for counterclockwise rotation, where θ is the angle of sector model as shown in Figure 2b. (2) Velocity Treatment. As the particle position is relocated, its dynamic property which may affect the calculation in the next time step has to be rotated, such as velocity. The velocity maintains its magnitude to rotate with a rotation angle α given below:

particles 1′ and 2. The force Fn acting on particle 1 has the same magnitude and direction of Fn′. (ii) Position and velocity treatment. When a particle center exceeds PB2 (e.g., particle 3′ at (x′3, y′3, z′3)), it will move to a new position (e.g., particle 3 at (x3, y3, z3)) with the identical dynamic parameters (e.g., velocities). The position relocation can be represented by x3 = x′3, y3 = y′3 + L, z3 = z′3, where L is the thickness of the model in the y direction. For a sector model shown in Figure 2a, the circumferential PBC are applied on two side walls of the sector as shown in Figure 2b. The principle is similar to the rectangular PBC, but the treatment is based on the rotation instead of translation. Therefore, circumferential PBC work through the following three aspects: position, velocity, and force treatments. (1) Position Treatment. A particle center moving out of one periodic boundary, e.g. particle 3′ at (x′, y′, z′), can be treated as moving into the other periodic boundary, e.g. particle 3. The position is relocated through the following equations: x = x′ cos α + y′ sin α

(3)

y = −x′ sin α + y′ cos α

(4)

z = z′

(5)

Vx = Vx ′ cos α + Vy ′ sin α

(6)

Vy = −Vx ′ sin α + Vy ′ cos α

(7)

Vz = Vz ′

(8)

where Vx, Vy, and Vz are the three components of particle velocity in the x, y, and z directions, respectively. (3) Force Treatment. The force acting on a particle is updated at each time step. Similar to the case in the rectangular PBC, the force has to be modified. Specifically, if a particle edge exceeds one periodic boundary while the center still inside, e.g. particle 1, its edge is assumed to enter the other periodic boundary and the assumed position can be calculated according to eqs 3−5, e.g. particle 1′. Thus, the interactions with particles near the other periodic boundaries, e.g. particle 2, have to be calculated with a rotation angle α. The calculated force F′n act on particle 1 can be transformed by the following equations: Fx = Fx ′ cos α + Fy ′ sin α

(9)

Fy = −Fx ′ sin α + Fy ′ cos α

(10)

Fz = Fz ′

(11)

where Fx, Fy, and Fz are the three components of particle normal contact force in x, y, and z directions, respectively. The above treatments are valid for particles not located in the central region of a sector model. That is, for particles in a small central region with a radial distance R (shown in Figure 3 where R is equal to the radius of particle k), the circumferential PBC cannot fully function. Here, the contacts on particle k can be divided into two groups. Group 1 are contacts that make particle k move toward to the sector center (e.g., green particles in the figure), and group 2 are contacts that make particle k move outward the periphery (e.g., blue particles). In a sector model,

Figure 10. Corrected spatial distribution of CN in the hopper flow (only left side shown).

Figure 11. Spatial distribution of the magnitudes of particle vertical velocities in the hopper flow at t = 0.7 s. 8248



Article

Figure 13. Normal contact force network among particles in a hopper flow at t = 0.7 s.

with unrealistic flow patterns observed in the central region. How to overcome this problem by developing suitable radial PBC will be a subject of research in the future. An alternative method is therefore devised here, done by setting a radial velocity rather than radial force to this particle. The aim is to prevent particles moving out of the computational domain without introducing an obvious interference. Hence, the method based on velocity is rather straightforward, focused on the radial velocity of a particle when its surface reaches the center. The transformation of the velocities can be expressed as Vr 0 = Vx 0 cos β + Vy0 sin β

(12)

Vt 0 = −Vx 0 sin β + Vy0 cos β

(13)

where Vx0 and Vy0 are particle velocities in the Cartesian coordinate system, Vr0 is the radial velocity, Vt0 is the tangential velocity, and β is the angle between particle center and x-axis, calculated by β = tan−1(y/x). If Vr0 < 0, it means that the particle moves to the sector center and may pass through. Thus, a velocity at Vr0 = 0 will be assigned to this particle, stopping its further movement in the radial direction. If the particle surface leaves the sector center, such a condition will be removed. Mathematically, this criterion is written as Vr 0 < 0: Vr = 0, Vt = Vt 0 (14)

Figure 12. Probability density distributions of vertical velocities (a), radial velocities (b), and tangential velocities (c).

group 2 particles do not exist, which are however important to prevent particle k from flowing out of the sector center. Therefore, there are two types of PBC for a sector model: circumferential PBC and radial PBC. In principle, the PBC based on the force and velocity, as used for the circumferential PBC, should also be applied to the radial PBC. Thus, once particle k tends to move radially out of the sector center, suitable force and velocity should effect radially but with opposite direction to make particle move back to the sector domain. However, because the particle does not have its “mirror particle” radially, the force magnitude and direction cannot be generated automatically. They have to be assigned artificially based on the contacts of the particle with other particles from the sector domain. This will result in varying results. For example, in this study, an attempt to apply the radial force and velocity on particle k from the contacts with group 2 particles is made, but the results are not satisfactory,

Vr 0 ≥ 0:

Vr = Vr 0,Vt = Vt 0

(15)

The modified radial and tangential velocities can be mapped back to the x−y plane: Vx = Vr cos β − Vt sin β

(16)

Vy = Vr sin β + Vt cos β

(17)

In the next section, the proposed theoretical treatments will be tested by two cases that are widely examined in the literature: flow in a cylindrical hopper and formation of a conical sandpile. This is done by comparing the simulated results between sector models with different sector angles and full 3D models. The parameters examined include flow pattern, velocity, coordination number (CN), and forces. 8249



Article

Table 2. Parameters Used in the Sandpile Simulation variables particle diameter particle number and running timea rolling friction (particle− particle) rolling friction (particle−wall) sliding friction (particle− particle) Sliding Friction (particle-wall) particle density Poisson ratio Young’s modulus time step

value 0.002 m 3D 90° 24 000 6000 20 h 5h 0.0001 m

60° 4000 5h

45° 3000 3h

30° 2000 3h

0.0001 m 0.6 0.5 2500 kg/m−3 0.3 1.0 × 107 Pa 1.43 × 10−5 s

a

All the cases presented in this work are calculated on the National Computational Infrastructure (NCI), a Sun Constellation Cluster with 1492 nodes, each containing 2 quad core Nehalem processors.

its design and control and also important to the development of granular dynamic theories.21−23 Consequently, DEM-based model has been used by many researchers, as reviewed by Zhu et al.3 A cylindrical hopper with a flat bottom is used in this work and shown in Figure 4. It is similar to that used by Zhu and Yu.24 The hopper is 20d in diameter (d is particle diameter), and 50d in height, with a circular orifice of diameter 8d at the central bottom. Different sector angles (θ = 30°, 45°, 60° and 90°) are used to compare with the full 3D results. Other parameters used in the simulation are listed in Table 1. Note that the number of particles is reduced in a sector model compared to the full 3D model, and consequently the simulation time decreases with the decrease of sector angle. In a run of simulation for a full 3D or sector model, particles with small random initial velocities are first allowed to settle into the hopper under gravity, giving a packing with height around 90d. Then the orifice is opened, permitting discharge of particles from the hopper until the residual particles have negligible velocities. For visualization, the vertical section with a thickness of 4d through the center of the full 3D model or a 3D-sector model is chosen. Figure 5 shows the discharging process of hopper flow achieved by the full 3D and sector models. For better visualization of flow patterns, different colors are set to divide the bed into different layers in the vertical direction. The mirror image of 3D-sector model is shown for comparison with the full 3D model. The full 3D model shows three distinct flow zones: plug flow zone in the upper part, V-shape converging flow zone in the central bottom, and stagnant zone in the corners. The same phenomenon can be observed in each sector model, demonstrating sector models can produce consistent results with the full 3D model. Figure 6 displays the horizontal plane of each model at the height of 0.1 m (shown as dash lines in Figure 5), which further confirms the consistence between the full 3D and sector models. To be quantitative, microscopic properties such as CN, velocity, and force structures are compared. These properties are not uniform, having spatial and hence statistical distributions. In this study, such microscopic results are collected at time 0.7 s after the orifice opening. CN is the number of particles in contact with a given particle, which can be used to characterize the structure of an assembly of particles. As shown in Figure 7, for all cases considered, CN mainly varies between 5 and 7. At the lower part near the orifice,

Figure 14. Probability density distribution of normal contact forces at t = 0.7 s in (a) radial, (b) tangential, and (c) vertical directions.

Figure 15. Geometry of cylindrical container for sandpile formation.

3. RESULTS AND DISCUSSION 3.1. Hopper Flow. Hoppers are widely applied in engineering practice. Comprehensive understanding of the dynamic behavior of granular flow in a hopper is essential for optimizing 8250



Article

Figure 16. Sandpile forming process in (a) full 3D, (b) 90°, (c) 60°, (d) 45°, and (e) 30° sector models.

Figure 17. Spatial distribution of CN at t = 0.14 s in the full 3D and sector models.

However, differences still exist. In the region near the central line, the CN in a sector model is lower than the full 3D model. This region becomes larger with the decrease of the sector angle. This is because some particles near the center have been removed in a sector model. Figure 8 shows the probability density distribution of CN. The similar distributions confirm the consistency between the full 3D model and sector models. Notably, the deviation is not obvious because the amount of particles in the center is too small to be reflected in the probability density distribution. The deviation of CN in the center is caused by the geometry of a sector. It can be fixed by relating this deviation to the geometry. The thickness of the sector in the radial direction can be illustrated by the length of the arc rθ, where r is the distance between particle center and the sector center with its unit being d. Figure 9 shows the correlations between the arc length rθ and CN for different sector models, demonstrating that when the arc length is smaller than 4d, the CN deviates from the average value CN0, the CN far from the center region. In this deviation region, there exists an approximately linear increase of CN with the arc length

Figure 18. Probability density distribution of CN at t = 0.14 s.

the value is much lower due to the converging flow region. Each sector model can provide a distribution consistent with the full 3D model, not only qualitatively but also quantitatively. 8251



Article

Figure 19. Spatial distribution of vertical velocities (m/s) at t = 0.14 s in the full 3D and sector models.

plug flow in the upper part. In the stagnant region, especially the part close to the orifice, some particles move upward due to the pressure from the downward flow stream. The amount of particles with positive vertical velocities is small, suggested by the low probability density in the positive axis. The radial velocity mainly falls around −0.005 m/s, indicating that most particles move extremely slowly in the radial direction toward to the hopper center. The distribution curves of large velocities are a bit higher in the negative axis, consistent with the fact that a small amount of particles above the orifice move fast to the outlet. The tangential velocity symmetrically distributes on the axis, with a peak value at around 0 m/s, implying the tangential movement is insignificant and equally balanced positively and negatively. The force structure, which can represent the particle connectivity, is an essential local characteristic of granular matter to depict the force transmission and other dynamic behavior of particles. In this work, only the normal contact forces between particles are considered as all other contact forces are related to it. A thin vertical slice (4d thick) is chosen for visualization. In such a sector, the particles in the center are not shown due to their contacts are not fully calculated for the reason discussed above. The thickness of each line connecting two particle centers is proportional to the magnitude of the force. As shown in Figure 13, all the force networks indicate three main features of hopper flow: large forces exist at the walls around the bottom corner, smaller forces in the upper part, and the smallest force in the region adjacent to the orifice where particles have fewer contacts. The probability density distribution of the total forces acting on a particle is also analyzed and shown in Figure 14. Here, the i total force is defined as Fi = ∑kj=1 Fn,ij. As done for velocities, the normal contact force is investigated by its three components: radial, tangential, and vertical. In the radial direction, the consistence between sector models and full 3D model keeps well. For each model, the predominant value of the radial force is very small, indicating small interactions between particles in the radial direction. The larger forces distribute a little more in the negative axis, suggesting the associated particles move toward the orifice in the center. In the tangential direction, the curves are almost symmetric on the vertical axis, implying the particles move randomly and insignificantly. In the vertical direction, the curves cover a large area in the negative axis with a peak value at −0.005 N,

(the solid line in Figure 9). Thus, an equation is proposed here to describe the relationship between CN and r: CN(r ) = krθ + b

(18)

where k is the slope of the line and b is a constant. According to Figure 9, the two constants can be estimated as k = (CN0/3) − 1 and b = 4 − (CN0/3). Correcting the deviated CN back to the average value can make the CN spatial distribution more realistic. Therefore, a correction factor c can be defined as

c(r ) = CN0 /CN(r )

(19)

That is, the correction factor c is a function of the arc length rθ where the particle locates. Therefore, the coordination number should be c·CN where CN is calculated based on the DEM simulation in a sector model. Note that eqs 18 and 19 are valid only when the arc length less than about 4.0d. When the arc length is more than 4.0d, the correction factor is equal to 1. Figure 10 shows the modified CN spatial distribution, confirming the method proposed here can reasonably, though not fully, eliminate the deviation of CN in the sector center. Particle velocity is another important parameter to validate the sector model proposed. Here, three components of a velocity are considered: radial (Vr), tangential (Vt), and vertical (Vz), calculated according to eqs 5 and 6. In the hopper flow, the vertical velocity is dominant compared with the radial and tangential velocities. Thus, only the spatial distribution of this velocity is shown in Figure 11. In the lower part, particles descend at a relatively high velocity at the orifice and the particles around the bottom corner of hopper have a very low velocity. Note that the region with high velocities around the orifice in the sector model of 30° is larger than the full 3D model and other sector models. This is caused by the solid flow converged at a high position above the orifice as flow patterns shown in Figure 5. In the upper part, the velocity fields of all models are consistent. The vertical velocity becomes smaller and almost steady. The probability density distributions of the radial, tangential and vertical velocities of particles are shown in Figure 12. The results further confirm the consistency of different models. In each model, the vertical velocity principally distributes in the negative axis with a peak value at 0.1 m/s, corresponding to the 8252



Article

which is the natural result due to the discharge of particles at the bottom. The results above suggest that, as long as the sector angle is greater than a certain value, sector models are valid. For the systems considered in the present work, the smallest angle lies in between 30 and 45°. However, this smallest angle is affected by the ratio of hopper diameter to particle size. With the increase of the ratio (e.g., obtained with a larger hopper diameter or smaller particle diameter), the effect of sector central regions becomes smaller, hence the smallest angle will decrease. How to generally determine the smallest angle for different cylinder diameter/ particle diameter ratios for different flow systems may need further studies, which may require a large number of computations to establish some relationships for general application. Interestingly, the deviation observed in the CN spatial distribution and force network in sectors almost have no effect on the dynamics of particles, even for the particles near the center. This demonstrates that with the boundary conditions used in this work, particles near the sector model can move in both vertical and horizontal directions as those in the full 3D model. 3.2. Sandpile Formation. The angle of repose is one important parameter in characterizing the behavior of granular materials. Many fundamental phenomena are related to it, e.g. avalanching,25 stratification,26 and segregation.27 DEM can explicitly take into account not only the geometrical factors but also the forces involved in the formation of a sandpile. It becomes a realistic method to examine the formation of sandpiles, as briefly reviewed by Zhu et al.2 The formation of a sandpile therefore offers another important operation to examine the applicability of the proposed circumferential PBC. In this work, a cylindrical container with a fixed plate bottom and an annulus outlet is used for simulations. The geometry is shown in Figure 15. Similar to the hopper study shown in section 3.1, four 3D-sector models (30°, 45°, 60°, and 90°) are used, in addition to the full 3D model. A simulation is conducted using the procedures, as used by Zhou et al.28,29 The particles randomly generated are first dumped into the container with the outlet closed. After a gravitational settling process, a stable packing is formed. Then the discharging process commences when the outlet opens instantaneously and ends with stationary particles remaining at the bottom plate forming a stable sandpile. The parameters and particle properties used in the simulation are listed in Table 2. The particle number and running time required in each model demonstrates the efficiency of the sector models compared with the full 3D model, a sector model may only take

Figure 20. Probability density distribution of (a) vertical velocities; (b) radial velocities; and (c) tangential velocities when t = 0.14 s.

Figure 21. Normal contact network in a pilling process for different sector models at t = 0.14 s. 8253



Article

10%−25% of the computing resource required by the full 3D model. The same treatment as used for hopper flow is adopted here for result analysis. A vertical slice with thickness of 4d is used to show the formation process of a sandpile. The mirror image is generated in each sector model to provide the full view of flow pattern. Five stages of pilling process at different times and one horizontal view at t = 0.36 s (at 0.01 m as shown by the dash lines) are shown in Figure 16 (discharge starts at t = 0 s). Clearly, the stable sandpiles formed are almost identical, and the patterns during the sandpile forming process demonstrate high similarity between the sector and full 3D models. The microscopic properties including CN, velocity, and normal contact force are analyzed at t = 0.14 s. As illustrated in Figure 17, the CN in each sector model matches that in the full 3D model well, with a dominant value of 5. Note that the same method to eliminate the deviation of CN for the hopper flow is adopted. In the discharging region, the value decreases as expected. This consistency is further verified by the probability density distribution shown in Figure 18. The only difference observed is that the peak value is small in the 30° sector. As demonstrated in the hopper flow, the CN near the sector center is low due to the lack of neighbors in a small sector. Generally, such difference is small in the CN probability density distribution, because the number of particles in the center is small. However, in the case of sandpile formation, the total particle number in the 30° sector is relatively small. Thus, the percentage of particles in the central region can be reflected by the probability density as shown in Figure 18. Figure 19 shows the spatial vertical velocity distributions, indicating the results are quite comparable. High velocities locate at the region above the outlet, and the low velocity zone has a similar profile to the sandpile formed in the following stage. In the area above the stagnant zone, particles move along with the profile of stagnant zone to the orifice with increasing the vertical velocity. The probability density distributions of particle velocities are also examined as done for hopper flow. As shown in Figure 20a, the highest peak for the vertical velocity distribution located at 0 m/s indicates nearly 25% particles descend with an extremely low velocity, which contributes to the stagnant area. Another small peak appears at −0.45 m/s, corresponding to the area near the outlet, where particles move quickly out of the orifice. For the radial velocity (Figure 20b), each sector model produces results similar to the full 3D model. The radial velocity predominantly distributes in the positive axis, which means that most particles move toward to the cylinder periphery. The peak value locates at around 0.005 m/s, corresponding to the stagnant area. For the tangential velocity distribution (Figure 20c), the consistency between different models is obvious. The distribution is symmetric at the vertical axis, implying particles move equally in the positive and negative tangential directions. A thin vertical slice with a thickness of 4d is selected to show the normal contact force structure (Figure 21). It can be observed that the magnitude of the contact force gradually propagates from the middle to bottom vertically and from outside to inside horizontally. The normal contact force nearly vanishes around the outlet, where particles are discharged mainly by gravity and less contacts. The probability density distribution of the normal contact forces is also examined and shown in Figure 22. It shows that the force probability density of each sector model is consistent with the full 3D model for all the three components of a contact force.

Figure 22. Probability density distribution of normal contact force of pilling process: (a) radial force, (b) tangential force, and (c) vertical force at t = 0.14 s.

Thus, the sector models can be used to describe the pilling process, providing consistent results with full 3D model; not only the flow pattern but also the resulted microdynamics quantitatively match the full 3D model.

4. CONCLUSIONS A sector DEM model is proposed to represent the full 3D cylindrical model geometry. This is achieved by setting circumstantial PBC with a velocity-relevant treatment in the center. The treatment based on particle velocity rather than contact force can minimize the interference of boundary settings on the solid flow. To test the treatment, four sector models with different sector angles (90°, 60°, 45°, and 30°) are used to directly compare with the results from the corresponding full 3D model. Two typical processes, mainly hopper flow and sandpile formation, are 8254



■

considered. The following conclusions can be drawn based on the results obtained. In the hopper flow, three main distinct flow zones of plug flow zone, V-shape converging flow zone, and stagnant zone in the corners can be identified in both sector models and full 3D model. This consistency is further confirmed by the microscopic properties such as CN, velocity, and normal force network by both spatial and statistical distributions. The vertical velocity is higher in the orifice region where particles move fast, and lower in the corner where the stagnant particles are located. The small CN and the low contact forces between particles can be observed in the region above the orifice in both processes, as a result of the fast solid flow. In the region of stagnant particles, e.g. the corners in hopper, the CN and the particles contact forces are relatively higher. It should be noted that for the small sector angles, some deviations exist in the central region for CN and normal force network, which is caused by the limited number of neighbors due to the narrow space in the center region. A general correction is proposed to reduce the deviation of CN in the center. In the pilling formation, the process of the formation and the angle of repose in the final state in the sector models are almost the same as in the full 3D model. Similar to the hopper flow, the consistence is further confirmed by the microdynamic analysis. The high velocity region locates near the orifice, correspondingly the CN and the contact forces are smaller. The low velocity region has a similar profile to the final sandpile, where the CN and contact forces are higher. Clearly, the proposed sector models can produce flow patterns and other flow properties consistent with the corresponding full 3D cylindrical model, demonstrating the applicability of the proposed PBC treatment. However, comparing with the full 3D cylindrical model, the number of particles used in the system and hence the computational cost can be significantly reduced. As demonstrated in this work, a substantial computing time is saved by around 80%−90% for the cases considered in the present work. It should be noted that the PBC treatment proposed in this work is not perfect yet. For example, it has a problem when applied to particles in the center. The problem may be overcome by developing suitable radial PBC. Moreover, there are other problems, although minor. For example, although the results are satisfactory under the conditions considered in this work, it is useful to determine a priori the smallest sector angle for different vessel-to-particle diameter ratios. Therefore, further studies are necessary to develop better PBC for sector models for general application.

■

Article

REFERENCES

(1) Cundall, P. A.; Strack, O. D. L. Discrete numerical-model for granular assemblies. Geotechnique 1979, 29, 47. (2) Zhu, H. P.; Zhou, Z. Y.; Yang, R. Y.; Yu, A. B. Discrete particle simulation of particulate systems: Theoretical developments. Chem. Eng. Sci. 2007, 62, 3378. (3) Zhu, H. P.; Zhou, Z. Y.; Yang, R. Y.; Yu, A. B. Discrete particle simulation of particulate systems: A review of major applications and findings. Chem. Eng. Sci. 2008, 63, 5728. (4) Kuang, S. B.; Zou, R. P.; Pan, R. H.; Yu, A. B. Gas-solid flow and energy dissipation in inclined pneumatic conveying. Ind. Eng. Chem. Res. 2012, 51, 14289. (5) Zhu, H. P.; Yu, A. B. Steady-state granular flow in a 3D cylindrical hopper with flat bottom: macroscopic analysis. Granular Matter 2005, 7, 97. (6) Bai, H.; Theuerkauf, J.; Gillis, P. A.; Witt, P. M. A coupled DEM and CFD simulation of flow field and pressure drop in fixed bed reactor with randomly packed catalyst particles. Ind. Eng. Chem. Res. 2009, 48, 4060. (7) Anand, A.; Curtis, J. S.; Wassgren, C. R.; Hancock, B. C.; Ketterhagen, W. R. Predicting discharge dynamics from a rectangular hopper using the discrete element method (DEM). Chem. Eng. Sci. 2008, 63, 5821. (8) Elperin, T.; Golshtein, E. Comparison of different models for tangential forces using the particle dynamics method. Phys. A: Stat. Mech. Appl. 1997, 242, 332. (9) Zhou, Z. Y.; Zhu, H. P.; Yu, A. B.; Wright, B.; Pinson, D.; Zulli, P. Discrete particle simulation of solid flow in a model blast furnace. ISIJ Int. 2005, 45, 1828. (10) Wright, B.; Zulli, P.; Zhou, Z. Y.; Yu, A. B. Gas-solid flow in an ironmaking blast furnace – I: Physical modelling. Powder Technol. 2011, 208, 86. (11) Cui, L.; O’Sullivan, C.; O’Neill, S. An analysis of the triaxial apparatus using a mixed boundary three-dimensional discrete element model. Geotechnique 2007, 57, 831. (12) Zhou, Z. Y.; Yu, A. B.; Choi, S. K. Numerical simulation of the liquid-induced erosion in a weakly bonded sand assembly. Powder Technol. 2011, 211, 237. (13) Mochen, N. W. B. Numerical modeling in micromechanics via particle methods. Proceedings of the 1st International PFC Symposium, Gelsenkirchen, Germany, Nov 6−8, 2002. (14) O’Sullivan, C.; Bray, J. D.; Li, S. F. A new approach for calculating strain for particulate media. Int. J. Numer. Anal. Methods Geomech. 2003, 27, 859. (15) Kaneko, Y.; Shiojima, T.; Horio, M. DEM simulation of fluidized beds for gas-phase olefin polymerization. Chem. Eng. Sci. 1999, 54, 5809. (16) Thornton, C. Numerical simulations of deviatoric shear deformation of granular media. Geotechnique 2000, 50, 43. (17) Grest, G. S.; Dünweg, B.; Kremer, K. Vectorized link cell Fortran code for molecular dynamics simulations for a large number of particles. Comput. Phys. Commun. 1989, 55, 269. (18) Rapaport, D. C. The event scheduling problem in molecular dynamic simulation. J. Comput. Phys. 1980, 34, 184. (19) Dziugys, A.; Peters, B. An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers. Granular Matter 2001, 3, 231. (20) Zhou, Z. Y.; Zhu, H. P.; Yu, A. B.; Wright, B.; Zulli, P. Discrete particle simulation of gas-solid flow in a blast furnace. Comput. Chem. Eng. 2008, 32, 1760. (21) Potapov, A. V.; Campbell, C. S. Computer simulation of hopper flow. Phys. Fluids 1996, 8, 2884. (22) Flekkøy, E. G.; Jørgen Måløy, K. Continuum description of granular flows: Simulation and experiment. Phys. Rev. E 1998, 57, 6962. (23) Zhang, J. Y.; Rudolph, V. Effect of shear friction on solid flow through an orifice. Ind. Eng. Chem. Res. 1991, 30, 1977. (24) Zhu, H.; Yu, A. Micromechanic modeling and analysis of unsteady-state granular flow in a cylindrical hopper. J. Eng. Math. 2005, 52, 307. (25) Lee, J. Avalanches in (1 + 1)-dimensional piles - a moleculardynamics study. J. Phys. I 1993, 3, 2017.

AUTHOR INFORMATION

Corresponding Author

*Tel: 61 2 9385 4429. Fax: 61 2 9385 6565. E-mail address: a. [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work has been financially supported by the Australian Research Council. W.Y. is grateful to the University of New South Wales for providing the Tuition Fee Scholarship and the China Scholarship Council for providing a stipend. All the cases presented in this work are calculated on the National Computational Infrastructure, which is funded by the Australian Government. 8255



Article

(26) Makse, H. A. Stratification instability in granular flows. Phys. Rev. E 1997, 56, 7008. (27) Dury, C. M.; Ristow, G. H.; Moss, J. L.; Nakagawa, M. Boundary effects on the angle of repose in rotating cylinders. Phys. Rev. E 1998, 57, 4491. (28) Zhou, Y. C.; Xu, B. H.; Yu, A. B.; Zulli, P. Numerical study of the angle of repose of granular materials. Phys. Rev. E 2001, 64, 021301. (29) Zhou, Y. C.; Xu, B. H.; Yu, A. B.; Zulli, P. A numerical and experimental study of the angle of repose of granular particles. Powder Technol. 2002, 125, 45.

8256


Periodic Boundary Conditions for Discrete Element Method Simulation

Recommend Documents