Discrete Element Method Simulations of the Inter-Particle Contact ...

5 downloads 19477 Views 5MB Size Report
May 11, 2017 - Keywords: angle of repose; contact parameters; iron ore particles; DEM .... systems within the EDEM software package (DEM—Solutions Ltd., Edinburgh, ... ct and n are the total mass of the clump, center of mass, and number ...
materials Article

Discrete Element Method Simulations of the Inter-Particle Contact Parameters for the Mono-Sized Iron Ore Particles Tongqing Li 1,2 , Yuxing Peng 1,2, *, Zhencai Zhu 1,2 , Shengyong Zou 3 and Zixin Yin 1,2 1 2 3

*

School of Mechatronic Engineering, China University of Mining and Technology, Xuzhou 221116, China; [email protected] (T.L.); [email protected] (Z.Z.); [email protected] (Z.Y.) Jiangsu Key Laboratory of Mine Mechanical and Electrical Equipment, China University of Mining and Technology, Xuzhou 221116, China State Key Laboratory of Mining Heavy Equipment, CITIC Heavy Industries Co., Ltd., Luoyang 471000, China; [email protected] Correspondence: [email protected]; Tel.: +86-138-0520-9649

Academic Editor: Yong-Cheng Lin Received: 22 March 2017; Accepted: 10 May 2017; Published: 11 May 2017

Abstract: Aiming at predicting what happens in reality inside mills, the contact parameters of iron ore particles for discrete element method (DEM) simulations should be determined accurately. To allow the irregular shape to be accurately determined, the sphere clump method was employed in modelling the particle shape. The inter-particle contact parameters were systematically altered whilst the contact parameters between the particle and wall were arbitrarily assumed, in order to purely assess its impact on the angle of repose for the mono-sized iron ore particles. Results show that varying the restitution coefficient over the range considered does not lead to any obvious difference in the angle of repose, but the angle of repose has strong sensitivity to the rolling/static friction coefficient. The impacts of the rolling/static friction coefficient on the angle of repose are interrelated, and increasing the inter-particle rolling/static friction coefficient can evidently increase the angle of repose. However, the impact of the static friction coefficient is more profound than that of the rolling friction coefficient. Finally, a predictive equation is established and a very close agreement between the predicted and simulated angle of repose is attained. This predictive equation can enormously shorten the inter-particle contact parameters calibration time that can help in the implementation of DEM simulations. Keywords: angle of repose; contact parameters; iron ore particles; DEM simulation; mills

1. Introduction With the rapid development of computing power and advanced contact algorithm, the discrete element method (DEM) has been extensively applied as a leading tool to describe diverse issues in granular processes, including many industrial applications, such as mining, chemical, cement, and agricultural industries, particularly tumbling mills [1–12]. DEM simulations hitherto have been demonstrated their extremely desirability to predict what happens in reality, as well as the quantitatively accurate information representation inside the mills, while the accuracy of outcomes for DEM simulations depends highly on the input parameters [13–18], including parameters such as contact parameters (restitution coefficient, static friction coefficient, and rolling friction coefficient), mechanical properties (shear modulus, Poisson’s ratio, and density), and particle shape. Researchers studied that the DEM simulation of mill behavior was only a weak variation with the mechanical properties [13–16,19–21]. To dramatically accelerate the computational efficiency of

Materials 2017, 10, 520; doi:10.3390/ma10050520

www.mdpi.com/journal/materials

Materials 2017, 10, 520

2 of 14

DEM simulations, the reduction in shear modulus was commonly employed without changing the charge behavior [15,22,23]. In the last few decades the determination of contact parameters of charge were mainly classified into two approaches [14]. The first approach was to determine the contact parameters separately using the experimental setups [24–26], whereas the restrictions on particle size and shape were difficult to be measured, especially for the inter-particle contact parameters at a small scale [14,26]. At present, the DEM parameters at the mesoscale can be truly determined on the basis of molecular dynamics simulations [27,28]. The second approach was to set a series of arbitrary input parameters in DEM simulations until the results were in very close agreement to the experimental results, namely, sand-pile calibrations [29–31]. The angle of repose, as the key parameter in describing the granule flow characteristic, has been widely employed in calibrating the contact parameters of granule materials. Numerous experiments and DEM simulations proposed that the angle of repose was highly sensitive to the contact parameters [13,14,16–18]. Generally speaking, the repeated calibrations consist of two parts, namely particle-wall contact parameters and inter-particle contact parameters. Therefore, the determination of contact parameters by means of the second approach is time-consuming and sophisticated. Currently, the contact parameters between particles and the wall were commonly determined with various experimental setups in the first place, and the repeated sand-pile test was subsequently employed in calibrating the inter-particle contact parameters, in order to reduce the number of calibrations. Iron ores are the primary source for iron- and steel-making industries, which use tumbling mills for further comminution to obtain the required particle size distribution. The shapes of iron ore particles are highly irregular due to the impact-breakage, abrasion, and attrition. Currently, the exact simulation of irregular shapes hitherto is still one of the key issues to be solved in the DEM simulations [32,33]. To trustworthily predict what happens in reality, rather than qualitative prediction, the accuracy of DEM simulations requires special consideration with the particle shape, due to the influence of packing density on the particle flow properties. Arguably, the sensitivity of the angle of repose dependence on the particle shape varied from a low value, for smooth granules, to a high value, for irregular granules [13,25,34,35]. Therefore, modelling the contact model with irregular shapes is essential and worthwhile for DEM simulations. Many investigators employed the simple spherical shape in the majority of DEM simulations to achieve a reasonable simulation time, but the shortcomings of unrealistic simplifications and assumptions are also clear. Such simplified models are an unrealistic measure on both the physical properties and contact parameters, as well as a poor particle representation on the packing characteristics. Along with the fast development of X-ray and image processing techniques, the multi-element spheres method has been used successfully to represent the geometrical model of irregular iron ore particles [36–38]. To achieve an exact representation of particle shape, a large number of spheres were used to describe the information of the details, but at the expense of much greater computational effort [39]. Thus, considering the accuracy in describing the shape of iron ore particles and computational effort, the appropriate geometrical shape of irregular particles should be carefully modelled. To reduce the simulations, as well as the computational effort, the determination of the contact parameters should be studied quantitatively for determining which parameter has the greatest impact on the angle of repose. In this paper, the irregular shapes of iron ore particles were determined, and a quick and accurate sphere clump method was employed in modelling the irregular shape. The appropriate number of spheres was quantitatively selected on the basis of the multi-element sphere method. The inter-particle contact parameters were systematically changed whilst the particle-wall contact parameters were arbitrarily assumed to purely assess their impact on the angle of repose using the split cylinder method. A predictive equation was formed, which will provide the basic data for the DEM simulations of tumbling mills.

Materials 2017, 10, 520

3 of 14

2. DEM Model The exact representation of irregular particles for DEM simulations in an efficient manner is still a major limitation in the current state for industrial applications. The exact definition of the irregular shape element hitherto is extremely sophisticated. Therefore, developing the contact model with the non-spherical methods was phenomenally complicated compared with the multi-element sphere methods. To achieve a better computational efficiency for contact detection, the multi-element sphere method was employed in modelling the approximation of any desired shape. In this study, the non-linear model, namely the no-slip Hertz–Mindlin model, was employed to figure out the contact force between particle systems within the EDEM software package (DEM—Solutions Ltd., Edinburgh, UK) because of the accurate representation of the physical situation and less computational effort. In EDEM software any kind of irregular geometrical shape of particles can be generated as a clump composed of several touching or overlapping spheres. Therefore, the contact detection between the sphere clumps is sphere-based and, therefore, the discrete element algorithm of the sphere clumps is fully available for calculating the contact forces. Although the simplified Hertz–Mindlin model has been addressed extensively, it is worthwhile to review the main equations again and the corresponding inputs are essential to the model. According to the Newton’s second law of motion, the translational and rotational motions of particle j (Figure 1) can be expressed as [31]: dv j mj = m j g + ∑ (Fn,ij + Ft,ij ) (1) dt i Ij

dw j wj = −µr,ij Fij R j + ∑ R j × Ft,ij dt wj i

(2)

where mj , vj , wj , Ij are mass, velocity, angular velocity, and the moment of inertia of particle j, respectively. F 2017,and F are the normal force and total tangential force between particle i and j Materialsn,ij 10, 520 t,ij 3 of 14 given as: r 1 multi-element sphere non-spherical methods was phenomenally p withthe 4 ∗ √ ∗ 3 complicated 20  compared 2 ∗ ∗ Fn,ija = − computational R∗ δdetection, vn,ij the multi-element (3) E R δ 2 + efficiency β mfor Econtact methods. To achieve better n 3 3 sphere method was employed in modelling the approximation of any desired shape. In this study, r 1 figure out the  was employed  to the non-linear model, namely the no-slip p Hertz–Mindlin80 model, p 2 ∗ ∗ (4) Ft,ijbetween = minparticle [µs,ij Fn,ij , 8G within δt RtheδnEDEM + software β m∗package G ∗ R∗(DEM—Solutions δn vt,ij ] contact force systems Ltd., 3 Edinburgh, UK) because of the accurate representation of the physical situation and less * = of 2 of particles can be computational In EDEM software any geometrical shape where E* is equivalenteffort. Young’s modulus, 1/Ekind (1 irregular − u1 2 )/E i + (1 − u2 )/Ej , Ei , Ej , ui , uj are the generated as a clump composed of several touching or overlapping spheres. Therefore, contact Young’s Modulus and Poisson’s ratio of particles i and particle j, respectively; R* isthe equivalent radius, detection between the sphere clumps is sphere-based and, therefore, the discrete element algorithm * * 1/R = 1/Rofi the + 1/R are available the contact radius ofthe particles i and particle j. m is equivalent mass, j , Rclumps i and R sphere isj fully for calculating contact forces. 1/m* = 1/mi +Although 1/mj , mthe and m are the mass of particles i and particle j. δn and δitt isisworthwhile the normal overlap simplified Hertz–Mindlin model has been addressed extensively, i j to review the main equations again and the corresponding inputs are essential to the and tangential overlap, respectively; µs,ij is the coefficient of static friction betweenmodel. particles i and to the Newton’s second law of motion, the translational and rotational motions of particle particle j. vAccording n,ij and vt,ij are the relative normal velocity and tangential velocity, respectively.

j (Figure 1) can be expressed as [31]:

Figure The illustration of the contact forces particle i and particle j. Figure 1. The 1.illustration of the contact forcesbetween between particle i and particle j.

mj

Ij

dw j dt

dv j dt

 m j g   ( Fn,ij  Ft ,ij )

 -μr ,ij Fij R j

(1)

i

wj w

  R j  Ft ,ij i

(2)

Materials 2017, 10, 520

4 of 14

As mentioned above, the discrete element algorithm of the sphere clumps is fully available for calculating the contact forces. The properties of total mass and products of inertia of the clump were given as follows [36]: n mct = ∑k=1 m[k] (5) xict = Iii =

1 n m [k ] xi [k ] mct ∑k=1

(6)



 2 2 m[k] ( x j [k] − x j ct ) + m[k] r [k] 5 o n 2 n Iij = ∑k=1 m[k] ( x j [k] − x j ct ) ; i 6= j n

∑ k =1

2 

(7) (8)

where mct , xi ct and n are the total mass of the clump, center of mass, and number of the balls, respectively, m[k] , x[k] , and r[k] are the mass, center of mass, and radius of kth ball, respectively. DEM employs the modelling of each particle as a rigid body and determining the displacement and velocity on the basis of Newton’s second law at a time step. To make a trustworthy simulation, a time steps, in the order of a millionth of a second, is required to predict what happens in reality, which leads to greater computationally expense. Researchers proposed that the appropriate time step for simulating the dense particle motion was in range of 20–80% [40,41]. The time step, defined as the time between the particle iterations, must be less than the critical Rayleigh time step ∆TR : ∆Tstep

πR p < ∆TR = (0.163v p + 0.8766)

s

ρp Gp

(9)

where Rp is the particle radius, νp the particle Poisson’s ratio, ρp the particle density, Gp the particle shear modulus. In the current study, 30% of the critical Rayleigh time step was used due to the consideration of computational effort and accuracy. 3. Methodology 3.1. Sphere Clump Method In mineral processing, the real shape of iron ore particles is highly irregular due to the impact-breakage, abrasion, and attrition. Generally speaking, the exact definition of real particle shapes hitherto is still one of the key issues to be solved. So far there are various definitions qualitatively describing the particle shape, including parameters such as aspect ratio, sphericity, and shape factor [42]. Nevertheless, sphericity indicates how closely the particle geometry is to a perfect sphere and is employed extensively in powder technology [36]: SAes Ψ= = SArp

√ 3

36πV 2 SArp

(10)

where SArp is the real surface area of the iron ore particle, mm2 ; SAes is the surface area of the sphere determined by the same volume of the iron ore particle, mm2 ; V is the real volume of the iron ore particle, mm3 . In the current study, the raw iron ores were obtained from an iron mine in Xuzhou with and iron grade of 67.46%, which are currently used for iron- and steel-making industries in China. The larger chunks of raw iron ore were processed three times by an industrial jaw crusher. Subsequently, the processed particles were reprocessed in a laboratory-scale jaw crusher to obtain the desired product size distribution. Finally, the products were sieved carefully (10 min) by a vibrating screen to obtain the required mono-sized particles. The mono-sized iron ore particles in the range of 2–4 mesh were selected for this study.

Materials 2017, 10, 520

5 of 14

the processed particles were reprocessed in a laboratory-scale jaw crusher to obtain the desired product size distribution. Finally, the products were sieved carefully (10 min) by a vibrating screen to obtain the required mono-sized particles. The mono-sized iron ore particles in the range of 2–4 Materials 2017, 10, 520 5 of 14 mesh were selected for this study. In the event that more accurate shape representation was desired, thirty-six particles were In the thatthe more shape representation thirty-six selected selected toevent generate 3Daccurate digitized and meshed shapeswas by adesired, high-accuracy 3Dparticles scannerwere (TEXU-BLU, to generate the 3D digitized and meshed shapes by a high-accuracy 3D scanner (TEXU-BLU, 2M pixel), 2M pixel), as shown in Table 1. The 3D scanner used in this study was of 50 × 38 mm in scanning as shown Tablebreadth, 1. The 3D used in this study wasmm of 50 × 38 mmaccuracy. in scanning range of a range of ain single 2 s scanner in scanning speed, and ±0.015 in scanning Subsequently, single breadth, 2 s in scanning speed, and ±was 0.015employed mm in scanning accuracy. Subsequently, quick a quick and accurate sphere clump method in generating the 3D shapes. Theaelegant and accurate sphere clump method was employed in generating the 3D shapes. The elegant software, software, named Automatic Sphere-clump Generator (ASG), was employed to realize this method named Automatic Sphere-clump Generator (ASG), wasinemployed to realize this methodinfor generating for generating a geometrical model that can be used DEM simulations, as shown Figure 2. The amethod geometrical modelinto thattwo can steps: be used in DEM simulations, as shown in Figure The 2. The method is is divided sphere detection and sphere optimization. sphere-clump divided twoby steps: sphere candidate detection and sphere optimization. Thethe sphere-clump method starts by methodinto starts detecting spheres and then refines solution using a non-linear detecting candidate spheres and then the solution non-linear least-squares optimization. least-squares optimization. Since therefines detection process using uses aarandomized vertex sampling method Since the detection process uses a randomized vertex sampling method to find spheres, different, to find spheres, different, but equally valid, sphere-clumps can be obtained on consecutive runs. In but equally valid, sphere-clumps can be quantitatively, obtained on consecutive runs. assessing In order to characterize order to characterize the sphere-clumps two parameters and describingthe the sphere-clumps quantitatively, twoerror parameters and describing modelling error, namely modelling error, namely volume and EITassessing error, were estimated. EITthe error shows the percentage volume EIT error, error were estimated. EIT error axes. showsVolume the percentage average distribution averageerror massand distribution along the principal error shows themass percentage error error along principal axes. error shows the percentage error between the mesh volume between thethe mesh volume andVolume clump volume, expressed as: and clump volume, expressed as:

Vm  V c Volum error V −  100% m VVmc × 100% Volum error =

(11) (11)

Vm where Vm is the mesh volume, and Vc is the clump volume. where Vm is the mesh volume, and Vc is the clump volume.



Figure 2. 2. Iron Iron ore ore particle particlemodelling modellingby bythe thesphere sphereclump clumpmethod method(60 (60spheres). spheres). Figure Table1.1. The The real realphysical physicalproperties propertiesof ofthirty-six thirty-sixparticles particlesusing usingaahigh-accuracy high-accuracy3D 3Dscanner. scanner. Table No 1 2 3 4 5 6 7 8 9 10 11 12 13 14

No x (mm) y (mm) z (mm) SArp (mm2) x (mm) y (mm)10.13z (mm)9.04 SArp 306.59 (mm2 ) 1 12.00 2 12.00 17.28 3 16.33 4 17.08 5 11.99 6 8.90 7 9.06 8 7.27 9 9.87 10 7.68 11 10.87 8.66 10.06 6.90

17.28 10.13 10.96 16.33 17.088.81 11.997.66 8.9011.10 7.80 9.06 10.32 7.27 6.81 9.876.98 7.688.26 10.879.99 4.68 7.24 6.68

10.96 9.04 9.32 8.81 9.32 10.56 7.66 10.5610.56 11.1010.56 9.00 7.80 9.00 6.57 6.57 10.32 8.40 8.40 6.81 7.92 7.92 6.98 6.58 6.58 8.26 7.16 7.16 9.99 5.76 5.76 9.18 7.52 5.60

384.21 306.59 384.21 362.99 362.99 356.01 356.01 315.12 315.12 177.12 177.12 207.58 207.58 157.98 157.98 190.06 190.06 164.40 164.40 216.31 216.31 196.37 188.76 138.72

V (mm3) SAes (mm2) Ψ 2) V (mm3 ) SAes (mm0.815 371.26 249.76 450.97 284.34 371.26 249.76 0.740 450.97 284.34 0.744 417.46 270.08 417.46 270.08 0.733 396.36 260.90 396.36 260.90 0.814 386.53 256.56 386.53 256.56 0.737 140.36 130.59 140.36 130.59 184.78 156.86 0.756 184.78 156.86 132.59 125.72 0.796 132.59 125.72 172.74 149.97 172.74 149.97 0.789 138.55 129.47 138.55 129.47 0.787 198.35 164.45 198.35 164.45 0.760 179.84 154.05 165.99 146.04 124.48 120.54

Ψ 0.815 0.740 0.744 0.733 0.814 0.737 0.756 0.796 0.789 0.787 0.760 0.785 0.77 0.87

Materials 2017, 10, 520

6 of 14

Table 1. Cont. No

x (mm)

y (mm)

z (mm)

SArp (mm2 )

V (mm3 )

SAes (mm2 )

Ψ

15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

9.78 7.25 7.73 14.14 13.98 10.51 14.23 14.40 11.70 10.23 9.44 14.01 10.33 12.54 10.26 9.04 10.84 9.02 7.10 6.94 8.66 8.11

5.50 13.93 15.17 7.09 8.13 8.15 8.60 10.72 13.77 8.57 11.33 7.77 6.63 6.90 7.30 8.21 6.59 7.13 8.25 8.27 6.82 7.29

8.02 10.87 9.13 8.76 8.80 8.50 9.19 5.92 6.36 10.07 5.86 7.67 10.22 7.06 9.09 7.40 6.823 6.232 7.374 8.773 7.359 6.48

177.04 330.06 335.25 272.36 319.42 246.43 340.58 289.56 328.38 272.33 254.76 258.21 237.68 230.27 209.29 207.62 190.17 167.01 171.66 178.29 167.92 146.37

145.33 358.60 350.82 283.56 409.57 265.3 383.07 272.20 343.20 328.38 284.11 239.39 236.24 217.36 203.88 211.68 150.06 144.34 141.85 163.69 148.69 120.78

133.65 244.05 240.511 208.69 266.66 199.65 255.03 203.08 237.02 230.14 208.96 186.41 184.78 174.80 167.49 171.73 136.54 133.05 131.51 144.69 135.71 118.14

0.75 0.74 0.72 0.77 0.83 0.81 0.749 0.701 0.722 0.845 0.820 0.722 0.777 0.759 0.800 0.827 0.718 0.797 0.766 0.812 0.808 0.807

3.2. Simulation Conditions and Input Parameters Currently, there are different measurements to determine the angle of repose of particles, including measuring methods, such as the fixed-cone method, fixed-height method, rotating cylinders, and titling method [25,43]. In general, the sand-pile height was measured directly to determine the angle of repose. However, there are some limitations on determining the angle of repose accurately, due to the friction between the particle and the wall. To determine the inter-particle contact parameters purely, a series of DEM simulations in this study were studied on the basis of the swing-arm slump test, originally used by Grima and Wypych [14,44]. Shown in Figure 3a, a base of 50 mm in height and 195 in diameter was fixed horizontally on the ground. A split cylinder of 300 mm in height and 100 mm in diameter was placed on the base. A certain mass of iron ore particles (3 kg) was poured into the cone with a generation rate of 5 kg/s and then the split cylinders were pulled away at a high velocity to avoid the friction between particle and wall. Correspondingly, particles would be collapsed under the natural force of gravitation forming a pile of materials, as shown in Figure 3b, accumulated from the center of the plane and gradually expanding to the boundary. The angle of repose, as shown in Figure 3c, is defined as the angle between the inclined surface between the sand-pile and the plane, and was used to qualitatively analyze the influence of inter-particle contact parameters on particle flow characteristics. Then angle of repose can be approximately calculated and expressed as: α AOR =

45 4 [ αi ] π i∑ =1

(12)

where αi (i = 1, 2, 3, 4) is the measured angle in the orthogonal axes with a unit of radians (rad). Considering the input parameters selected in the EDEM software package (DEM—Solutions Ltd., Edinburgh, UK), the mechanical properties and contact parameters should be determined accurately. The mechanical properties of iron ore particles were obtained using the hydraulic MTS machine, and the mechanical properties of the cone were summarized from the literature. The iron ore particles, in principle, do not have exactly identical mechanical properties because of their distinct microstructure. However, it is extensively recognized that all iron ore particles in this study have the same material

Materials 2017, 10, 520

7 of 14

Materials 2017, 10, 520

7 of 14

properties and contact reductions in the assignment of EDEM input parameters. have the same materialparameters, properties due and to contact parameters, due to reductions in the assignment of Researchers presented that the shear modulus was commonly artificially reduced to achieve EDEM input parameters. Researchers presented that the shear modulus was commonly reasonable artificially simulation [22]. Hence, the shear modulus usedHence, in this the study formodulus the iron ore particles has been reduced to times achieve reasonable simulation times [22]. shear used in this study for reduced by a factor of 100 relative to the typical value, and found that the reduced value has no effect the iron ore particles has been reduced by a factor of 100 relative to the typical value, and found that on angle of repose. parameters forCurrently, the wall-particle system can be measured thethe reduced value hasCurrently, no effect the on contact the angle of repose. the contact parameters for the accurately, butsystem the inter-particles contact parameters to determine.contact To reduce the number wall-particle can be measured accurately, are butdifficult the inter-particles parameters are of DEM simulations, the contact parameters for the wall-particle system were arbitrarily assumed to difficult to determine. To reduce the number of DEM simulations, the contact parameters for the be constant, whereas the contact parameters fortointer-particles selected in a range from a low wall-particle system were arbitrarily assumed be constant, were whereas the contact parameters for value to a high value (Table 2). inter-particles were selected in a range from a low value to a high value (Table 2).

Figure 3. Forming process of the particle pile on the plane. (a) Initial state of the simulation; (b) being Figure 3. Forming process of the particle pile on the plane. (a) Initial state of the simulation; (b) being in the the particle’s particle’s accumulating; accumulating;(c) (c)sand-pile sand-pilestabilized. stabilized. in Table 2. Input parameters for EDEM simulations. Table 2. Input parameters for EDEM simulations.

Material Parameters Material ParticleParameters density (kg m−3)

Symbols Symbolsρp

Particle m−3 ) (Gpa) Particledensity shear (kg modulus ParticleParticle shear modulus (Gpa) Poisson’s ratio Particle Poisson’s ratio −3 Wall density (kg m ) Wall density (kg m−3 ) Wall shear modulus Wall shear modulus (Gpa)(Gpa) Wall Poisson’s Wall Poisson’s ratio ratio Particle-wall restitution coefficient Particle-wall restitution coefficient Particle-wall static friction coefficient Particle-wall static friction coefficient Particle-wall rolling friction coefficient Particle-wall rolling friction coefficient Particle-particle restitution coefficient Particle-particle coefficient Particle-particle static restitution friction coefficient Particle-particle rolling friction coefficient Particle-particle static friction coefficient

ρ p Gp Gp νp νp ρw ρw Gw G w νw νw epw epw µs-pw μs-pw µr-pw epp μr-pw µs-pp epp µr-pp μs-pp

Value Value 3886 3886 2.587 2.587 0.283 0.283 1200 1200 1.05 1.05 0.41 0.41 0.5 0.5 0.6 0.6 0.05 0.05 0–0.6 0–0.6 0–0.8 0–0.2 0–0.8

Particle-particle rolling friction coefficient μr-pp 0–0.2 4. Results and Discussion 4. Results and Discussion 4.1. Particle Shape Estimation 4.1. Particle Shape Estimation To determine the mechanical and flow behavior of iron ore particles, thirty-six iron ore particles are selected study thethe irregularly shape. geometrical of iron ore particles areore determined To to determine mechanical andThe flow behavior ofdimensions iron ore particles, thirty-six iron particles on the basis of the 3D digitized and meshed shapes. Considering the equivalent sphere with the same are selected to study the irregularly shape. The geometrical dimensions of iron ore particles are volume of theonreal surface areshapes. calculated in Table 1,the as shown in Figure 4. determined theparticle, basis ofthe theequivalent 3D digitized and areas meshed Considering equivalent sphere The surface area versus the real surface area is fitted byareas the method of least in squares. The withequivalent the same volume of the real particle, the equivalent surface are calculated Table 1, as fitting has a4. sphericity of 0.718surface and a correlation of 0.963area which is slightly smaller than shownvalue in Figure The equivalent area versuscoefficient the real surface is fitted by the method of is casesquares. for the average value (0.776), both are of less thanand the asphericity of acoefficient sphere. It of is demonstrated least The fitting value has abut sphericity 0.718 correlation 0.963 which is slightly smaller than is case for the average value (0.776), but both are less than the sphericity of a sphere. It is demonstrated that the iron ore particle is highly irregular comparing with the sphere so

Materials 2017, 10, 520 Materials 2017, 10, 520 Materials 2017, 10, 520

8 of 14 8 of 14 8 of 14

that theiron simplifications and assumptions of the irregular particle as that simple spherical shape is that ore particle isand highly irregular comparing with theparticle sphere as so the simplifications that the the simplifications assumptions of the irregular simple spherical shapeand is unrealistic. assumptions unrealistic. of the irregular particle as simple spherical shape is unrealistic. In this this study, the no. 13 and 3333are selected to represent the 2 ×24×mesh mono-sized particles, due In study, the theno. no.13 13and and33 are selected represent 4 mesh mono-sized particles, In this study, are selected to to represent thethe 2 × 4 mesh mono-sized particles, due to the consideration between the average sphericity and the equivalent radius. Other particle sizes due toconsideration the consideration between the average sphericity andequivalent the equivalent radius. particle to the between the average sphericity and the radius. OtherOther particle sizes are scaled on theonbasis of theof selected particle size, as shown in Table 3. Based on the on sphere clump sizes are scaled the basis the selected particle size, as shown in Table 3. Based the are scaled on the basis of the selected particle size, as shown in Table 3. Based on the sphere sphere clump methods, the 3D shapes are generated with a with cluster of spheres to represent the irregular model. clump methods, 3D shapes are generated a of cluster of spheres to represent the irregular methods, the 3D the shapes are generated with a cluster spheres to represent the irregular model. Studies show that the accuracy of the irregular particle increases with the increasing number of model. showthe that the accuracy of irregular the irregular particle increases with theincreasing increasingnumber number of of StudiesStudies show that accuracy of the particle increases with the spheres, butatat the expense of much greater computational effort. Therefore, considering the spheres, expense of much greater computational effort. Therefore, considering the accuracy spheres, but but atthethe expense of much greater computational effort. Therefore, considering the accuracy in describing the shape of iron ore particles and computational effort, different numbers of in describing the shapethe of iron ore and computational effort, different numbersnumbers of spheres accuracy in describing shape ofparticles iron ore particles and computational effort, different of spheres are selected and summarized in Figure 5. It is interesting to see that the accuracy in are selected summarized in Figure 5.inItFigure is interesting see that the describing spheres are and selected and summarized 5. It istointeresting toaccuracy see that inthe accuracythe in describing the geometrical shape of iron ore particles increases with increase in the number of geometrical shape of iron ore particles increases with increase in the number of spheres. Therefore, describing the geometrical shape of iron ore particles increases with increase in the number of spheres. Therefore, considering the accuracy in describing the shape of iron ore particles and considering the accuracy in describing the shape iron ore particles andofcomputational effort,and the spheres. Therefore, considering the accuracy inofdescribing the shape iron ore particles computational effort, the appropriate sphere clump numbers should be carefully selected. appropriate sphere clump numbers should be carefully selected. computational effort, the appropriate sphere clump numbers should be carefully selected.

Figure 4. 4. Sphericity Sphericity of of thirty-six thirty-six iron iron ore particles. particles. Figure Figure 4. Sphericity of thirty-six iron ore ore particles.

Figure 5. The geometrical model of iron ore particles with various numbers of sphere clumps. Figure 5. The geometrical model of iron ore particles with various numbers of sphere clumps. Figure 5. The geometrical model of iron ore particles with various numbers of sphere clumps. Table 3. Percent of the particle size used in the simulations. Table 3. Percent of the particle size used in the simulations. Table 3. Percent of the particle size used in the simulations. 3

Volume Intervals (mm ) Percent Volume Intervals (mm3) Percent 100–200 3 44.44% Percent Volume Intervals 100–200(mm ) 44.44% 200–300 25% 100–200 44.44% 200–300 25% 300–400 22.22% 200–300 25% 300–400 22.22% 400–500 8.33% 300–400 22.22% 400–500 8.33% 400–500

8.33%

In this model, the volume error and EIT error describing the details between the 3D model and In this model, the volume error and EIT error describing the details between the 3D model and the geometrical model are obtained, shown in Figure 6. It is also evident from Figure 6 that the In this model, the volume error and EIT error describing details between 3D model the geometrical model are obtained, shown in Figure 6. It isthe also evident from the Figure 6 that and the change in volume error and EIT error with number of sphere clumps decreases dramatically and the geometrical model areand obtained, shown in number Figure 6.of It is also evident Figuredramatically 6 that the change change in volume error EIT error with sphere clumpsfrom decreases and then remains approximately unchanged. Once the number of sphere clumps is greater than thirty, in volume error and EIT error unchanged. with numberOnce of sphere clumps decreases and then then remains approximately the number of sphere dramatically clumps is greater thanremains thirty, the maximum volume error and EIT error are 1.3% and 0.9%, and the minimum volume error and approximately unchanged. Once the number of sphere clumps is greater than thirty, the maximum the maximum volume error and EIT error are 1.3% and 0.9%, and the minimum volume error and EIT error are 0.2% and 0.1%, respectively. Hence, to achieve computational accuracy without EIT error are 0.2% and 0.1%, respectively. Hence, to achieve computational accuracy without

Materials 2017, 10, 520

9 of 14

Materials 2017, 10, 520 Materials 2017, 10, 520

9 of 14 9 of 14

volume error and EIT error are 1.3% and 0.9%, and the minimum volume error and EIT error are 0.2% increasing the amount amount of of calculation, calculation, the computational geometrical model model usedwithout in the the work work is aa cluster cluster of thirty thirty and 0.1%, respectively. Hence, to achieve accuracy increasing the amount of increasing the the geometrical used in is of spheres. calculation, the geometrical model used in the work is a cluster of thirty spheres. spheres.

Figure 6. The change in volume error and EIT error with number number of of sphere sphere clumps. Figure 6. 6. The The change change in in volume volume error error and and EIT EIT error with with Figure error number of sphere clumps. clumps.

To avoid avoid the the friction friction between between the the particle particle and and wall, wall, the the split split cylinder cylinder should should be be pulled pulled away away To quickly enough. Figure 7 presents the variation of angle of repose as a function of velocity for the quickly enough. enough. Figure Figure77presents presentsthe the variation of angle of repose a function of velocity the variation of angle of repose as a as function of velocity for thefor given given contact parameters. parameters. The The geometry geometry of of the the sand-pile sand-pile and and the the angle angle of repose repose are are only only aa weak weak given contactcontact parameters. The geometry of the sand-pile and the angle of reposeofare only a weak variation variation with the varying of the velocity if the velocity is greater than 0.5 m/s. Hence, in this study, variation with theofvarying of theifvelocity if theisvelocity greater than Hence, 0.5 m/s.in Hence, in thisa study, aa with the varying the velocity the velocity greateristhan 0.5 m/s. this study, velocity velocity of 1 m/s is carried out to conduct a series of DEM simulations. velocity m/s is out carried out to conduct DEM simulations. of 1 m/sof is 1carried to conduct a series aofseries DEMof simulations.

Figure 7. Effect of velocity on angle of speed when epp = 0.05, μs-pp = 0.15, μr-pp = 0.01. Figure7.7.Effect Effectofofvelocity velocityon onangle angleofofspeed speedwhen whenepp epp==0.05, 0.05, µμs-pp s-pp = Figure = 0.15, 0.15, μµr-pp 0.01. r-pp= =0.01.

4.2. Effect of of Inter-Particle Contact Contact Parameters 4.2. 4.2. Effect Effect of Inter-Particle Inter-Particle Contact Parameters Parameters To examine the impact of contact parameters on on angle angle of of repose, repose, varying varying inter-particle inter-particle contact contact To examine the the impact impact of of contact contact parameters parameters To examine on angle of repose, varying inter-particle contact parameters are are systematically changed. changed. The The variation variation of of the the angle angle of of repose repose as as aa function function of of the parameters parameters are systematically systematically changed. The variation of thethe angle ofofrepose asplotted a function of the the restitution coefficient is primarily studied. Figure 8 presents angle repose against the restitution coefficient is primarily studied. Figure 88 presents the angle of repose plotted against the restitution coefficient is primarily studied. Figure presents the angle of repose plotted against the restitution coefficient coefficient for for the given given rolling rolling friction friction coefficient coefficient and and static friction friction coefficient. coefficient. It It is is clear clear restitution restitution coefficient for the the givenlittle rolling friction coefficient and static static friction coefficient. It isrange clear that the angle of repose varies with the changing restitution coefficient over the that theangle angleof of repose varies little with the changing restitution coefficient over considered. the range that the repose varies little with the changing restitution coefficient over the range considered. Additionally, the angle of repose determined by the high friction coefficient is much much considered. Additionally, the angle of reposeby determined by the coefficient high friction coefficient is Additionally, the angle of repose determined the high friction is much higher than is higher than than is is case case for for the the low low friction friction coefficient. coefficient. It It is is apparent apparent that that the the angle angle of of repose repose has has strong strong higher case for theto low coefficient. It is coefficient, apparent that of repose strong sensitivity to the sensitivity thefriction rolling/static friction butthe theangle impacts of the thehas rolling friction coefficient sensitivity to the rolling/static friction coefficient, but the impacts of rolling friction coefficient rolling/static friction coefficient, but the impacts of the rolling friction coefficient and static friction and static static friction friction coefficient coefficient on on the the angle angle of of repose repose may may be be aa single single factor factor or or an an interrelated interrelated and and and coefficient on the angle of repose may be a single factor or an interrelated and interdependent effect. interdependent effect. effect. As As seen in in Figure Figure 8, 8, it it is is clear clear that that the the impacts impacts of of the the rolling rolling friction friction coefficient coefficient interdependent As seen infriction Figure coefficient 8, it is seen clear that the impacts ofa the rolling frictionThe coefficient and static friction and static are interrelated and combined effect. restitution coefficient, as aa and static friction coefficient area interrelated and aThe combined effect. The restitution coefficient, as coefficient are interrelated and combined effect. restitution coefficient, as a basic property of basic property property of of collision collision energy, energy, representing representing the the energy energy dissipation dissipation during during the the collisions, collisions, is is used used basic collision energy, the energy dissipation during the collisions, is used to determine the to determine determine the representing damping coefficient coefficient of inter-particle inter-particle and particle-wall. particle-wall. Thus, varying the damping damping to the damping of and Thus, varying the coefficient slightly slightly leads leads to to the the variation variation of of angle angle of of repose. repose. Investigators Investigators proposed proposed that that the the coefficient restitution coefficient coefficient is is not not sensitively sensitively affected affected by by the the collision collision velocity velocity so so that that the the fixed fixed value value of of the the restitution

Materials 2017, 10, 520

10 of 14

damping coefficient of inter-particle and particle-wall. Thus, varying the damping coefficient slightly 14 of angle of repose. Investigators proposed that the restitution coefficient10 isof not 10 of 14 sensitively affected by the collision velocity so that the fixed value of the restitution coefficient is restitution commonly used the Therefore, aa constant value commonly coefficient used in theis simulations. constant value of the restitution is restitution coefficient isDEM commonly used in inTherefore, the DEM DEM asimulations. simulations. Therefore, constantcoefficient value of of the the restitution coefficient is used for the rest of the DEM simulations. used for the rest of the DEM simulations. restitution coefficient is used for the rest of the DEM simulations. Materials 520 leads to2017, the 10, variation Materials 2017, 10, 520

Figure Figure 8. 8. Effect Effect of of the the restitution restitution coefficient Figure 8. Effect of the restitution coefficient on on the the angle angle of of repose. repose.

Figure Figure 99 shows shows the the angle angle of of repose repose as as aa function function of of the the static static friction friction coefficient coefficient under under the the given given restitution coefficient and rolling friction coefficient. The angle of repose is strongly dependent restitution coefficient and rolling friction coefficient. The The angle of repose is strongly dependent on on the static friction coefficient, forming the sand-pile higher, friction the static friction coefficient, forming the sand-pile higher, with with an an increase increase in in the the static static friction coefficient. coefficient. It It is is recognized recognized that that the the DEM is dependent on the the soft DEM simulation simulation is dependent on soft contact contact approach, approach, and and the deformation will be generated during the interaction collisions between particles. static the deformation will be generated during the interaction interaction collisions collisions between between particles. particles. Larger Larger static friction friction coefficients coefficients can can tolerate tolerate with with aa larger larger elastic elastic deformation deformation in in the the tangential tangential direction, direction, so so it it has has great potential to form the pile higher for a large static friction coefficient. Additionally, increasing great potential potential to toform formthe thepile pilehigher higherfor fora large a large static friction coefficient. Additionally, increasing static friction coefficient. Additionally, increasing the the rolling friction coefficient for static friction coefficient increases of the rolling friction coefficient for the the given given static friction coefficient increases the angle angle of repose, repose, and rolling friction coefficient for the given static friction coefficient increases the the angle of repose, andand the the investigation of the of friction coefficient on the of be the investigation of impact the impact impact of rolling rolling friction coefficient onangle the angle angle of repose repose will be presented presented investigation of the of rolling friction coefficient on the of repose will bewill presented in the in the following section. in the following section. following section.

Figure 9. The effect of the static friction coefficient on Figure on the the angle angle of of repose. repose. Figure 9. 9. The The effect effect of of the the static static friction friction coefficient coefficient on the angle of repose.

Obviously, Obviously, the the rolling rolling friction friction coefficient coefficient and and static static friction friction coefficient coefficient are are both both significant significant in in Obviously, the rolling friction coefficient and static friction coefficient are both significant in determining determining the the angle angle of of repose, repose, increasing increasing rolling rolling friction friction coefficient coefficient can can also also obviously obviously increase increase the the determining the angle of10 repose, increasing rolling friction coefficient can also obviously increase the the angle angle of of repose. repose. Figure Figure 10 shows shows the the angle angle of of repose repose affected affected by by rolling rolling friction friction coefficient coefficient for for the given given restitution restitution coefficient coefficient and and static static friction friction coefficient. coefficient. Increasing Increasing the the rolling rolling friction friction coefficient coefficient evidently increases the height of sand-pile, while the bottom width of sand-pile evidently increases the height of sand-pile, while the bottom width of sand-pile decreases decreases correspondingly. correspondingly. There There is is aa strong strong dependence dependence with with the the angle angle of of repose repose forming forming higher higher as as the the rolling rolling friction friction coefficient coefficient increases. increases. It It is is mainly mainly because, because, in in the the contact contact model, model, the the rolling rolling friction friction coefficient coefficient

Materials 2017, 10, 520

11 of 14

angle of repose. Figure 10 shows the angle of repose affected by rolling friction coefficient for the given restitution coefficient and static friction coefficient. Increasing the rolling friction coefficient evidently increases the height of sand-pile, while the bottom width of sand-pile decreases correspondingly. There is2017, a strong Materials 10, 520dependence with the angle of repose forming higher as the rolling friction coefficient 11 of 14 increases. It is mainly because, in the contact model, the rolling friction coefficient actually gives an impactful torque resistancetorque to control the rotational motion of individual particles. A large rolling actually gives an impactful resistance to control the rotational motion of individual particles. friction coefficient will provide a large resistance force to control the rotational motion of particles A large rolling friction coefficient will provide a large resistance force to control the rotational and, therefore, thereand, is more potential to form the potential sand-pile to of form iron ore higher than case for motion of particles therefore, there is more theparticles sand-pile of iron oreisparticles the lowthan rolling largelyfriction due to coefficient, the reduction of a large of kinetic of energy of higher is friction case forcoefficient, the low rolling largely due amount to the reduction a large particle systems. amount of kinetic energy of particle systems.

Figure on the the angle angle of of repose. repose. Figure 10. 10. The The effect effect of of the the rolling rolling friction friction coefficient coefficient on

4.3. Formulation of A Predictive Equation 4.3. Formulation of A Predictive Equation The above results show that the angle of repose is highly sensitive to the static/rolling friction The above results show that the angle of repose is highly sensitive to the static/rolling friction coefficient, while it appears to be invariant to the restitution coefficient. The simulated angle of coefficient, while it appears to be invariant to the restitution coefficient. The simulated angle of repose repose is fitted to the following relationship using a non-linear equation, giving: is fitted to the following relationship using a non-linear equation, giving:

α

 56.61  e

0.01

μ0.084μ0.25

0.01 0.084 AOR= 56.61 · e pp r  ppµ0.25 s  pp α AOR µr − pp pp s− pp

(13) (13)

It is interesting to note that varying the combination of static friction coefficient and rolling It is interesting to note that varying the combination of static friction coefficient and rolling friction friction coefficient can provide an effective prediction to determine the angle of repose of iron ore coefficient can provide an effective prediction to determine the angle of repose of iron ore particles. particles. Based on the non-linear equation, the simulated angle of repose is plotted against the Based on the non-linear equation, the simulated angle of repose is plotted against the predicted angle predicted angle of repose, as shown in Figure 11. It is apparent that a very close agreement between of repose, as shown in Figure 11. It is apparent that a very close agreement between the predicted the predicted and simulated value is attained. The maximum error is approximately 1.801 degrees, and simulated value is attained. The maximum error is approximately 1.801 degrees, which is almost which is almost within the scope of errors. within the scope of errors. To investigate further, the indices of each parameter in this equation can also provide the information to determine which inter-particle contact parameter has the greatest impact on the angle of repose. It is clear that the index of the rolling friction coefficient (0.084) and static friction coefficient (0.25) are evidently greater than is the case for the restitution coefficient (0.01). Therefore, it can also be demonstrated that the impact of the restitution coefficient can be neglected whilst the impact of the rolling friction coefficient and static friction coefficient are both significant, compared to the restitution coefficient, in determining the angle of repose, whereas the impact of the static friction coefficient is more profound than that of the rolling friction coefficient. The results are the same to that reported by researchers [45], indicating the applicability of non-linear equations. However, the geometries of iron ore particles are highly irregular, rather than the spherical balls, and the sphere clump method

Figure 11. The effect of the simulated angle of repose with the predicted angle of repose.

The above results show that the angle of repose is highly sensitive to the static/rolling friction coefficient, while it appears to be invariant to the restitution coefficient. The simulated angle of repose is fitted to the following relationship using a non-linear equation, giving: Materials 2017, 10, 520

αAOR  56.61  e pp 0.01μr0.084 μ0.25  pp s  pp

(13) 12 of 14

It is interesting to note that varying the combination of static friction coefficient and rolling is exclusively employed in determining its geometries. is essentialthe to angle recognize that particles in friction coefficient can provide an effective prediction toItdetermine of repose of iron ore granular processes are usually non-spherical, of varying humidity and particle size in nature. Therefore, particles. Based on the non-linear equation, the simulated angle of repose is plotted against the the investigations irregular particles more spheres. This predictive predicted angle ofofrepose, as shown inhave Figure 11. universalities It is apparent than that athose veryofclose agreement between equation significantly shortened theisinter-particle contact parameter time, which help the predicted and simulated value attained. The maximum error calibration is approximately 1.801 can degrees, in the implementation of EDEM simulations. which is almost within the scope of errors.

Figure angle of of repose repose with with the the predicted predicted angle angle of of repose. repose. Figure 11. 11. The The effect effect of of the the simulated simulated angle

To investigate further, the indices of each parameter in this equation can also provide the 5. Conclusions information to determine which inter-particle contact parameter has the greatest impact on the angle The results presented above contribute to a better understanding of the impact of inter-particle contact parameters of mono-sized iron ore particles, such as the restitution coefficient, rolling friction coefficient, and static friction coefficient, on the angle of repose on the basis of discrete element method (DEM). In particular, it was shown that the iron ore particle was highly irregular with a sphericity of 0.718, and then the particles were modelled by a sphere clump method. The angle of repose was found to show a weak function of the restitution coefficient, whereas it had strong sensitivity to the static friction coefficient and rolling friction coefficient with a combined effect. Therefore, the fixed value of the restitution coefficient is commonly used in the DEM simulations. Furthermore, to address which parameter has the greatest impact on the angle of repose, a predictive equation was established to determine the inter-particle contact parameters. A very close agreement between the predicted and simulated angle of repose was attained. It was clear that the indices of the rolling friction coefficient (0.084) and static friction coefficient (0.25) were evidently greater than is the case for the restitution coefficient (0.01). Notably, the static friction coefficient plays a primary role on the angle of repose, followed by a less, but still profound, role for the rolling friction coefficient. The predictive equation can be very helpful to determine the inter-particle contact parameters of iron ore particles and can also significantly shorten the calibration time. Acknowledgments: The research reported here was supported by the National Nature Science Foundation of China (grant no. 51475458), key project of National Natural Science Foundation of China (grant no. U1510205), the Program for Changjiang Scholars and Innovative Research Team in University (IRT _16R68). The authors also wish to thank the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD) and the Top-notch Academic Programs Project of Jiangsu Higher Education Institutions (TATP). Author Contributions: Yuxing Peng proposed the research topic; Shengyong Zou performed the DEM simulations; Zhencai zhu guided and supervised the entire project; Tongqing Li wrote the manuscript and analyzed the data; and Yuxing Peng and Zixin Yin made the revisions of the manuscript. Conflicts of Interest: The authors declare no conflict of interest.

Materials 2017, 10, 520

13 of 14

References 1. 2. 3. 4. 5.

6. 7. 8. 9. 10. 11. 12. 13. 14.

15. 16.

17.

18. 19. 20. 21.

22. 23.

Guo, Y.; Curtis, J.S. Discrete element method simulations for complex granular flows. Ann. Rev. Fluid Mech. 2015, 47, 21–46. [CrossRef] Bhargava, K.C.; Thompson, B.; Malmstadt, N. Discrete elements for 3D microfluidics. Proc. Natl. Acad. Sci. USA 2014, 111, 15013–15018. [CrossRef] [PubMed] Majidi, B.; Taghavi, S.M.; Fafard, M.; Ziegler, D.P.; Alamdari, H. Discrete element method modeling of the rheological properties of coke/pitch mixtures. Materials 2016, 9, 334. [CrossRef] Lee, J.; Yun, T.S.; Choi, S.U. The effect of particle size on thermal conduction in granular mixtures. Materials 2015, 8, 3975–3991. [CrossRef] Pennec, F.; Alzina, A.; Tessier-Doyen, N.; Nait-Ali, B.; Mati-Baouche, N.; De Baynast, H.; Smith, D.S. A combined finite-discrete element method for calculating the effective thermal conductivity of bio-aggregates based materials. Int. J. Heat Mass Transf. 2013, 60, 274–283. [CrossRef] Hare, C.; Ghadiri, M.; Guillard, N.; Bosworth, T.; Egan, G. Analysis of milling of dry compacted ribbons by distinct element method. Chem. Eng. Sci. 2016, 149, 204–214. [CrossRef] Chen, Y.; Munkholm, L.J.; Nyord, T.A. Discrete element model for soil–sweep interaction in three different soils. Soil Tillage Res. 2013, 126, 34–41. [CrossRef] Ghodki, B.M.; Goswami, T.K. DEM simulation of flow of black pepper seeds in cryogenic grinding system. J. Food Eng. 2017, 196, 36–51. [CrossRef] Müller, P.; Tomas, J. Simulation and calibration of granules using the discrete element method. Particuology 2014, 12, 40–43. [CrossRef] Bracey, R.J.; Weerasekara, N.S.; Powell, M.S. Performance evaluation of the novel multi-shaft mill using DEM modeling. Miner. Eng. 2016, 98, 251–260. [CrossRef] Cleary, P.W. A multiscale method for including fine particle effects in DEM models of grinding mills. Miner. Eng. 2015, 84, 88–99. [CrossRef] Cleary, P.W.; Owen, P.J. Using DEM to understand scale-up for a HICOM® mill. Miner. Eng. 2016, 92, 86–109. [CrossRef] Zhou, Y.C.; Xu, B.H.; Yu, A.B.; Zulli, P. An experimental and numerical study of the angle of repose of coarse spheres. Powder Technol. 2002, 125, 45–54. [CrossRef] Grima, A.P.; Wypych, P.W. Investigation into calibration of discrete element model parameters for scale-up and validation of particle-structure interactions under impact conditions. Powder Technol. 2011, 212, 198–209. [CrossRef] Ng, T.T. Input parameters of discrete element methods. J. Eng. Mech. 2006, 132, 723–729. [CrossRef] Yan, Z.; Wilkinson, S.K.; Stitt, E.H.; Marigo, M. Discrete element modelling (DEM) input parameters: Understanding their impact on model predictions using statistical analysis. Comput. Part. Mech. 2015, 2, 283–299. [CrossRef] Wilkinson, S.K.; Turnbull, S.A.; Yan, Z.; Stitt, E.H.; Marigo, M. A parametric evaluation of powder flowability using a Freeman rheometer through statistical and sensitivity analysis: A discrete element method (DEM) study. Comput. Chem. Eng. 2016, 97, 161–174. [CrossRef] Coetzee, C.J. Calibration of the discrete element method and the effect of particle shape. Powder Technol. 2016, 297, 50–70. [CrossRef] Cleary, P.W. Charge behaviour and power consumption in ball mills: Sensitivity to mill operating conditions, liner geometry and charge composition. Int. J. Miner. Process. 2001, 63, 79–114. [CrossRef] Franke, J.; Cleary, P.W.; Sinnott, M.D. How to account for operating condition variability when predicting liner operating life with DEM–A case study. Miner. Eng. 2015, 73, 53–68. [CrossRef] Delaney, G.W.; Cleary, P.W.; Morrison, R.D.; Cummins, S.; Loveday, B. Predicting breakage and the evolution of rock size and shape distributions in AG and SAG mills using DEM. Miner. Eng. 2013, 50, 132–139. [CrossRef] Lommen, S.; Schott, D.; Lodewijks, G. DEM speedup: Stiffness effects on behavior of bulk material. Particuology 2014, 12, 107–112. [CrossRef] Malone, K.F.; Xu, B.H. Determination of contact parameters for discrete element method simulations of granular systems. Particuology 2008, 6, 521–528. [CrossRef]

Materials 2017, 10, 520

24. 25. 26.

27. 28. 29. 30.

31. 32.

33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45.

14 of 14

Barrios, G.K.P.; de Carvalho, R.M.; Kwade, A.; Tavares, M. Contact parameter estimation for DEM simulation of iron ore pellet handling. Powder Technol. 2013, 248, 84–93. [CrossRef] Wang, W.; Zhang, J.; Yang, S.; Zhang, H.; Yang, H.; Yue, G.G. Experimental study on the angle of repose of pulverized coal. Particuology 2010, 8, 482–485. [CrossRef] Just, S.; Toschkoff, G.; Funke, A.; Djuric, D.; Scharrer, G.; Khinast, J.; Knop, K.; Kleinebudde, P. Experimental analysis of tablet properties for discrete element modeling of an active coating process. AAPS PharmSciTech 2013, 14, 402–411. [CrossRef] [PubMed] Ostanin, I.; Ballarini, R.; Potyondy, D.; Dumitrica, T. A Distinct Element Method for Large Scale Simulations of Carbon Nanotube Assemblies. J. Mech. Phys. Solids 2013, 61, 762–782. [CrossRef] Ostanin, I.; Ballarini, R.; Dumitrica, T. Distinct Element Method Modeling of Carbon Nanotube Bundles with Intertube Sliding and Dissipation. J. Appl. Mech. 2014, 81, 061004. [CrossRef] Li, Y.; Xu, Y.; Thornton, C. A comparison of discrete element simulations and experiments for ‘sandpiles’ composed of spherical particles. Powder Technol. 2005, 160, 219–228. [CrossRef] Nakashima, H.; Shioji, Y.; Kobayashi, T.; Aoki, S.; Shimizu, H.; Miyasaka, J.; Ohdoi, K. Determining the angle of repose of sand under low-gravity conditions using discrete element method. J. Terramech. 2011, 48, 17–26. [CrossRef] Chen, H.; Liu, Y.L.; Zhao, X.Q. Numerical investigation on angle of repose and force network from granular pile in variable gravitational environments. Powder Technol. 2015, 283, 607–617. [CrossRef] Bourcier, D.; Féraud, J.P.; Colson, D.; Mandrick, K.; Ode, D.; Brackx, E.; Puel, F. Influence of particle size and shape properties on cake resistance and compressibility during pressure filtration. Chem. Eng. Sci. 2016, 144, 176–187. [CrossRef] Barua, S.; Yoo, J.W.; Kolhar, P.; Mitragotri, S. Particle shape enhances specificity of antibody-displaying nanoparticles. Proc. Natl. Acad. Sci. USA 2013, 110, 3270–3275. [CrossRef] [PubMed] Ileleji, K.E.; Zhou, B. The angle of repose of bulk corn stover particles. Powder Technol. 2008, 187, 110–118. [CrossRef] Cleary, P.W. The effect of particle shape on simple shear flows. Powder Technol. 2008, 179, 144–163. [CrossRef] Majidi, B.; Azari, K.; Alamdari, H.; Ziegler, D. Simulation of vibrated bulk density of anode-grade coke particles using discrete element method. Powder Technol. 2014, 261, 154–160. [CrossRef] Majidi, B.; Melo, J.; Fafard, M.; Alamdari, H. Packing density of irregular shape particles: DEM simulations applied to anode-grade coke aggregates. Adv. Powder Technol. 2015, 26, 1256–1262. [CrossRef] Caulkin, R.; Tian, W.; Pasha, M.; Jia, X. Impact of shape representation schemes used in discrete element modelling of particle packing. Comput. Chem. Eng. 2015, 76, 160–169. [CrossRef] Ketterhagen, W.R. Modeling the motion and orientation of various pharmaceutical tablet shapes in a film coating pan using DEM. Int. J. Pharm. 2011, 409, 137–149. [CrossRef] [PubMed] O’Sullivan, C.; Bray, J.D. Selecting a suitable time step for discrete element simulations that use the central difference time integration scheme. Eng. Comput. 2004, 21, 278–303. [CrossRef] Chung, Y.C.; Ooi, J.Y. Influence of discrete element model parameters on bulk behavior of a granular solid under confined compression. Part. Sci. Technol. 2008, 26, 83–96. [CrossRef] Mora, C.F.; Kwan, A.K.H. Sphericity, shape factor, and convexity measurement of coarse aggregate for concrete using digital image processing. Cem. Concr. Res. 2000, 30, 351–358. [CrossRef] Geldart, D.; Abdullah, E.C.; Hassanpour, A.; Nwoke, L.C.; Wouters, I. Characterization of powder flowability using measurement of angle of repose. China Part. 2006, 4, 104–107. [CrossRef] Grima, A.P.; Wypych, P.W. Development and validation of calibration methods for discrete element modeling. Granul. Matter 2011, 13, 127–132. [CrossRef] Zhou, Y.C.; Xu, B.H.; Yu, A.B.; Zulli, P. Numerical investigation of the angle of repose of monosized spheres. Phys. Rev. E 2001, 64, 1–8. [CrossRef] [PubMed] © 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).