DEM Investigation of the Influence of Particulate Properties and Operating Conditions on the Mixing Process in Rotary Drums: Part 1—Determination of the DEM Parameters and Calibration Process

This paper’s goal was to select methods and a calibration procedure which would lead to the determination of relevant parameters of a discrete element method (DEM) and virtual material creation. Seven particulates were selected with respect to their shape (spherical and non-spherical), size and density. The first calibration experiment involved “packing test” to determine the shape accuracy and bulk density of virtual packed particulates. The series of simulations were compared with real experiments, and the size, shape and density of virtual particles were optimized. Using three apparatuses, the input parameter values were experimentally determined for a contact model that defines the behavior of particulates in DEM simulations. The research part of the paper examines the influence of factors such as particle number; pile formation method; and the method of evaluation of the angle of repose on the process of the calibration of virtual material. The most reproducible results were achieved by the “pilling” method and by the rotating drum—both evaluated by the geometric method. However, it is always advisable to make an overall visual comparison of the slope shape between the calibration simulation and the experimental curves. The bowl’s diameter to particle size ratio should be greater than 25, and the calibration experiment should contain approximately 4000 particles to ensure representative results during angle of repose calibration experiment.


Introduction
The selection of suitable methods and procedures of experimental measurements focused on calibrating DEM models is crucial for individual processes. A suitable method should produce sufficiently reliable values and reflect differences in the properties of interparticle contact and particulate materials used. Such measurement methods should be reproducible using basic laboratory tests to determine interaction coefficients. These coefficients tend to be crucial for the discrete element method (DEM) calibration of the mechanical, physical and process properties of materials. The calibration applications based on the determination of AoR (angle of repose) [1][2][3][4][5][6] are the most common. The basic input values that are set in the simulations include density, shear modulus and Poisson's ratio [7,8].
Once these values are set, it is necessary to determine interaction coefficients (static and rolling friction, restitution), which can significantly influence the static and dynamic flow, and are very often used as calibration parameters [9,10].
The DEM input values may not reflect the real measured values exactly. The most common deviations are found in the shape and size of the particles used in DEM [11][12][13][14]. The resulting static The aim of study is to show a complex approach to DEM calibration of spherical and non-spherical particles via different calibration experiment. Complete calibration methodologies with some new findings are presented in the conclusions of the paper. The procedures used are designed for solid homogenization processes but they can be used also in many others areas of industry. They are focused on flowability and interacting behavior, which are important for homogenization processes.

Materials and Methods
The DEM model used in this study uses the "soft-sphere" method originally developed by Cundall and Strack [37]. With this method, the particles in contact are able to withstand small deformations, and these deformations are used to calculate the forces acting between the particles. Tsuji et al. (1992) proposed a nonlinear contact model based on an adaptation of the original model proposed by Cundall and Strack [38]. This widely used contact model of Hertz-Mindlin is described in commercial software EDEM by Equations (1) and (2). This model is particularly suitable for the simulation of non-cohesive particulate materials.
E* is the equivalent Young's modulus of two colliding particles, R* is the equivalent contact radius, δ ij n refers to the normal displacement of particles under the influence of the normal force, m* is the equivalent mass of particles and v ij n is the normal component of relative velocity. The normal contact stiffness is then calculated as k n = 2E * √ Rδ n . Damping coefficient t is a function of the coefficient of restitution e and ranges from 0 (absolutely viscous) to 1 (absolutely elastic). The tangential force F ij t is given by the tangential displacement δ ij t , the relative tangential velocity v ij t and the tangential stiffness k t = 8G * √ Rδ t . In EDEM, the tangential force is limited by the condition defined by Coulomb's law of friction.

Material
Seven different particulate materials were used in this paper. These materials were selected to represent particles of different sizes, shapes and densities. The particulate materials are shown in Figure 1. The selected materials include basic geometric bodies, such as spheres, cubes and cylinders that reflect the morphology of real particles used in different industries. For example, the spherical shape represents a number of agricultural products, such as peas, chickpeas, nuts and seeds. Similarly, the cylindrical shape illustrates pellets utilized in the energy industry (biomass processing), while the cubic shape epitomizes, for example, sharp-edged particles of construction aggregate, ore or coal.

Materials and Methods
The DEM model used in this study uses the "soft-sphere" method originally developed by Cundall and Strack [37]. With this method, the particles in contact are able to withstand small deformations, and these deformations are used to calculate the forces acting between the particles. Tsuji et al. (1992) proposed a nonlinear contact model based on an adaptation of the original model proposed by Cundall and Strack [38]. This widely used contact model of Hertz-Mindlin is described in commercial software EDEM by Equations (1) and (2). This model is particularly suitable for the simulation of non-cohesive particulate materials. (2) E* is the equivalent Young's modulus of two colliding particles, R* is the equivalent contact radius, δij n refers to the normal displacement of particles under the influence of the normal force, m* is the equivalent mass of particles and vij n is the normal component of relative velocity. The normal contact stiffness is then calculated as = 2 * . Damping coefficient t is a function of the coefficient of restitution e and ranges from 0 (absolutely viscous) to 1 (absolutely elastic). The tangential force Fij t is given by the tangential displacement δij t , the relative tangential velocity vij t and the tangential stiffness = 8 * . In EDEM, the tangential force is limited by the condition defined by Coulomb's law of friction.

Material
Seven different particulate materials were used in this paper. These materials were selected to represent particles of different sizes, shapes and densities. The particulate materials are shown in Figure 1. The selected materials include basic geometric bodies, such as spheres, cubes and cylinders that reflect the morphology of real particles used in different industries. For example, the spherical shape represents a number of agricultural products, such as peas, chickpeas, nuts and seeds. Similarly, the cylindrical shape illustrates pellets utilized in the energy industry (biomass processing), while the cubic shape epitomizes, for example, sharp-edged particles of construction aggregate, ore or coal. Steel, aluminum alloy, glass and plexiglass (PMMA) were used as the construction (contact) material of the test equipment. The particulate material included wood (ash wood), plastic (ABS) and steel. The values of mechanical properties (density, Young's modulus, Poisson's ratio) entered into the used Hertz-Mindlin numerical model for DEM simulations are shown in Table 2. The tabulated values of density and Young's modulus served as default values and were further refined and calibrated based on the resulting values for individual materials and particles. Poisson's ratio is used in the calculation only to calculate the effective contact Young's modulus and the shear modulus. The Steel, aluminum alloy, glass and plexiglass (PMMA) were used as the construction (contact) material of the test equipment. The particulate material included wood (ash wood), plastic (ABS) and steel. The values of mechanical properties (density, Young's modulus, Poisson's ratio) entered into the used Hertz-Mindlin numerical model for DEM simulations are shown in Table 2. The tabulated values of density and Young's modulus served as default values and were further refined and calibrated based on the resulting values for individual materials and particles. Poisson's ratio is used in the calculation only to calculate the effective contact Young's modulus and the shear modulus. The influence of Poisson's ratio is negligible in terms of the behavior of bulk material. In the case where the exact value is not known, Poisson's ratio is usually defined as υ = 0.3 [36]. One of the most commonly used methods for describing complex particle geometries is to merge multiple spherical surfaces to form a more complex particle (clump). Due to the grain shape of our selected particulate materials, this paper focuses on the shape accuracy of the P6 × 15 wood pellets, particularly the C6 wood cubes, which represent a considerable challenge for defining the optimum shape. For this reason, a number of different virtual particles were created for these two particulate materials, the use of which will be investigated in terms of storage, flowability and computing time. Figure 2 shows an overview of the shapes of the created virtual particles used in the following experiments. The abbreviation system of virtual particles is as follows: particle_number of spheres which the particle is made of-variation; e.g., C6_8s-a means a C6 particle made by 8 particles clumped in variant a. influence of Poisson's ratio is negligible in terms of the behavior of bulk material. In the case where the exact value is not known, Poisson's ratio is usually defined as υ = 0.3 [36]. One of the most commonly used methods for describing complex particle geometries is to merge multiple spherical surfaces to form a more complex particle (clump). Due to the grain shape of our selected particulate materials, this paper focuses on the shape accuracy of the P6 × 15 wood pellets, particularly the C6 wood cubes, which represent a considerable challenge for defining the optimum shape. For this reason, a number of different virtual particles were created for these two particulate materials, the use of which will be investigated in terms of storage, flowability and computing time. Figure 2 shows an overview of the shapes of the created virtual particles used in the following experiments. The abbreviation system of virtual particles is as follows: particle_number of spheres which the particle is made of-variation; e.g., C6_8s-a means a C6 particle made by 8 particles clumped in variant a.   Table 3 shows the number of particles used in the experiments, and their total weight from which the weight of one particle was determined. Subsequently, the volume of each of the created virtual particles and the specific density were determined so that at a given volume of a virtual particle, its weight had a realistic value. From several variants of cubic and cylindrical particles, one was selected for each material which best represents the properties of the particulate material and is the least computationally demanding. The following chapters describe the methodology of several different calibration experiments, with the help of which the definitions of particular particulate materials were created and optimized.

Determination of Interaction Coefficients
The interaction coefficients define how particles of particulate matter interact with each other or with other materials in contact. The three basic DEM interaction coefficients used in the Hertz-Mindlin contact model include the coefficient of static friction µs, the coefficient of rolling friction µr and the coefficient of restitution e. While some values of these interaction coefficients (e.g., static friction between a particle and a contact material) can be measured quite easily using test apparatus, the experimental determination of the coefficients defining interactions between the particles themselves (and their subsequent use in DEM) is often challenging. In order to measure all three interaction parameters, a test apparatus designed with regard to accuracy and simplicity of measurements was constructed. To measure the coefficient of restitution, the device features a height-adjustable arm with a head. A vacuum source can be connected to this head, allowing one to suck a test particle onto a fine screen of the head. The tilting arm controlled by a screw facilitates a fine pitch for smooth operation serves to measure the friction parameters. The static friction measurement scheme is shown in Figure 3a. The coefficient of static friction between the particle and the contact material µ s-pw is determined from the measured angle α s-pw using the following equation: Grains of a shape approaching a sphere (i.e., grains start to roll rather than slip on the wall material) are more challenging. This phenomenon has to be eliminated by clumping several particles into one compact unit (clump). In the particle surface thusly formed (see Figure 3d), the free rotation of the grains is prevented, and only static resistance, not rolling resistance, is measured. The same methodology was used to measure the initial static friction values between particles µ s-pp , but the contact material was replaced with the surface of another particle clump.
The coefficient of rolling friction has been widely discussed, as evidenced by a number of publications dealing with this topic [16,43]. When spherical particles are used, rolling friction should be included, but this is not necessary when working with non-spherical particles [36]. Wensrich and Katterfeld (2012) discuss this issue in their publications, among others. The same methodology was used for the spherical shape particles as for the measurement of static friction. The rolling friction coefficient for the particle µr is calculated from the measured angle α r according to the following equation: The coefficient of restitution (denoted by e) is the ratio of the final to the initial relative velocity between two objects after they collide. It usually ranges from 0 to 1, where 1 refers to a perfectly elastic collision. The coefficient of restitution can be considered as the degree to which the mechanical energy of bodies is maintained when reflected from a surface or other body. In cases where the frictional forces can be neglected, and in the perfect progression of the collision and bounce, the ratio between the potential energy before the fall E P1 and the potential energy after the reflection E P2 can generally be calculated as follows: For experimental measurements of the particle-wall restitution coefficient, an apparatus for measuring the interaction coefficients was used ( Figure 3b). A vacuum source was attached to a holder mounted on an extendable arm to hold a test particle on the holder screen. After switching off the vacuum source, the particle dropped from the height h 1 vertically onto a contact material, from which it bounced to the height h 2 . The entire course of the experiment was captured using a high-speed camera. The movement was evaluated using tracking software that allows one to accurately track the particle trajectory. After reading the respective drop and rebound heights, the coefficient of restitution was determined according to Equation (5). A method published by Hlosta et al. (2018) was used to measure the coefficient of restitution of two particles. The coefficient of restitution of two particles is based on the laws of momentum conservation and energy conservation. When two particles collide, the total impulse introduced to the system (the two particles) is equal to zero. In the case of the collision of two particles that travel in one plane at different velocities before and after collision, it is assumed that the contact force between the particles is also in the same plane and the state of momentum is maintained. In this case, the coefficient of restitution is defined as the ratio of the relative velocities of the particles before and after they collide. If modelling of the flow of bulk material is required, damping behavior often becomes a factor of secondary importance. This is because many bulk materials demonstrate relatively strong damping properties. In most cases (e.g., in friction dominates processes), a value between 0.2 and 0.4 is appropriate and drop tests can be neglected in order to reduce the number of parameters to be calibrated [36].
Particle A is hanged on a hinge and released from height h A1 . At the lowest point, particle A collides with particle B. It transmits its kinetic energy to particle B, which is reflected to the height h B2 . The coefficient restitution of two particles can be determined based on the Equation (5). Another way to determine the coefficient of restitution is using the ratio of the angle of reflection α B2 to the angle of drop α A1 . The principle of this measurement is shown in Figure 3c. Here again, the coefficient of restitution is a measure of how much energy is lost in a collision. The values range is 0 < e < 1. In an ideal plastic collision (e pp = 0), the particles remain together after they collide. In a perfect elastic collision (e pp = 1), no energy is lost so that the kinetic energy of the system is maintained and the velocity of the particles after the collision is equal to the velocity preceding the collision. In this study, a two-particle collision experiment using a double pendulum is employed to measure the coefficient of restitution of two particles [44]. All measurements of interaction parameters were repeated ten times for each particulate material. From these ten values, the average value was determined including the standard deviation.
Processes 2020, 8, x FOR PEER REVIEW 7 of 27 a two-particle collision experiment using a double pendulum is employed to measure the coefficient of restitution of two particles [44]. All measurements of interaction parameters were repeated ten times for each particulate material. From these ten values, the average value was determined including the standard deviation. Figure 3. Interaction coefficient determination schemes: (a) static friction; (b) particle-wall coefficient of restitution (CoR); (c) particle-particle CoR; (d) wall and particle clump samples.

Bulk Density Measurement and Calibration
Particles packed in a bed are probably the simplest state of particulate matter, which is present both in nature and in a number of industries. The understanding of the particle deposition is also important in industrial applications due to the particle compaction (which is utilized in powder metallurgy), and the characterization of porosity, since air permeability and thermal conductivity of particulate materials can be quantitatively linked to the pore content and inter-particle contact properties. The bulk density verification parameter was the bed height of a particulate material of a known number of particles in a cylinder during the experiment and the DEM simulation. Differences in the resulting height can be resolved by shape accuracy and particle size. The weight of each individual DEM particle was always identical to the weight of actual particles according to Table 3. The second use of the packing test is the initial estimation and determination of the friction . Interaction coefficient determination schemes: (a) static friction; (b) particle-wall coefficient of restitution (CoR); (c) particle-particle CoR; (d) wall and particle clump samples.

Bulk Density Measurement and Calibration
Particles packed in a bed are probably the simplest state of particulate matter, which is present both in nature and in a number of industries. The understanding of the particle deposition is also important in industrial applications due to the particle compaction (which is utilized in powder metallurgy), and the characterization of porosity, since air permeability and thermal conductivity of particulate materials can be quantitatively linked to the pore content and inter-particle contact properties. The bulk density verification parameter was the bed height of a particulate material of a known number of particles in a cylinder during the experiment and the DEM simulation. Differences in the resulting height can be resolved by shape accuracy and particle size. The weight of each individual DEM particle was always identical to the weight of actual particles according to Table 3. The second use of the packing test is the initial estimation and determination of the friction parameters µ s and a µ r , as the values of static and rolling friction affect the height of a bed [45]. A 1000 mL cylinder, stand, hopper, length gauge, tripod camera and slide shutter were used for the packing test. Each of the seven particulate materials was gradually placed in the hopper above the measuring cylinder and closed by the slide. The slide was then removed, and the material was freely poured into the cylinder. The camera took a photograph which was later used for reading the height of the material bed using the length gauge that was part of the apparatus and was placed in a vertical plane intersecting the axis of the cylinder. This procedure was repeated ten times for each material. From these ten values, the average was determined including the standard deviation. Based on DEM, the computational cost re-calculated to seconds per particle was also recorded as a possible criterion for the subsequent selection of shape accuracy for C6 and P6 × 15 particles. Material constants shown in Table 2 and specific densities of particles in Table 3 were used as input parameters. The scheme and principle used for the packing test are shown in Figure 4. cylinder and closed by the slide. The slide was then removed, and the material was freely poured into the cylinder. The camera took a photograph which was later used for reading the height of the material bed using the length gauge that was part of the apparatus and was placed in a vertical plane intersecting the axis of the cylinder. This procedure was repeated ten times for each material. From these ten values, the average was determined including the standard deviation. Based on DEM, the computational cost re-calculated to seconds per particle was also recorded as a possible criterion for the subsequent selection of shape accuracy for C6 and P6 × 15 particles. Material constants shown in Table 2 and specific densities of particles in Table 3 were used as input parameters. The scheme and principle used for the packing test are shown in Figure 4.

Shear Modulus Value Optimization
Using Young's modulus and a shear modulus, respectively, which are lower than their actual values, is common practice in DEM. The main reason is to reduce the computational time. Chen et al.
(2017) presented a study on the relationship between Young's modulus in the range 0.01-0.0001 of the real value of Young's modulus E0 and the results of DEM simulations in which mixing of particles inside a rotary drum was investigated. The results showed that Young's modulus has a significant influence on the behavior of individual particles during collisions, but only a small effect when mixing individual layers. Convective mixing due to random particle collisions is negligible compared to diffusion and friction mechanisms. The study also showed that a value of Young's modulus reduced by up to 0.0007 times of the real value does not affect the results of the mixing simulation, and the DEM models with these values corresponded to the experiments performed. Additionally, the effect of the drum rotation speed or filling level on the relationship between the Young's modulus value and the simulated mixing results was not observed in Chen's study [46]. In the case of DEM, isotropic materials, material properties, Young's modulus E, the shear modulus G and Poisson's ratio υ are related as follows: One of the key numbers in DEM simulation is the Rayleigh time step. This is the time taken for a shear wave to propagate through a solid particle. It is therefore a theoretical maximum time step for a DEM simulation of a quasi-static particulate collection in which the coordination number (total number of contacts per particle) for each particle remains above 1. It is given by:

Shear Modulus Value Optimization
Using Young's modulus and a shear modulus, respectively, which are lower than their actual values, is common practice in DEM. The main reason is to reduce the computational time. Chen et al.
(2017) presented a study on the relationship between Young's modulus in the range 0.01-0.0001 of the real value of Young's modulus E 0 and the results of DEM simulations in which mixing of particles inside a rotary drum was investigated. The results showed that Young's modulus has a significant influence on the behavior of individual particles during collisions, but only a small effect when mixing individual layers. Convective mixing due to random particle collisions is negligible compared to diffusion and friction mechanisms. The study also showed that a value of Young's modulus reduced by up to 0.0007 times of the real value does not affect the results of the mixing simulation, and the DEM models with these values corresponded to the experiments performed. Additionally, the effect of the drum rotation speed or filling level on the relationship between the Young's modulus value and the simulated mixing results was not observed in Chen's study [46]. In the case of DEM, isotropic materials, material properties, Young's modulus E, the shear modulus G and Poisson's ratio υ are related as follows: One of the key numbers in DEM simulation is the Rayleigh time step. This is the time taken for a shear wave to propagate through a solid particle. It is therefore a theoretical maximum time step for a DEM simulation of a quasi-static particulate collection in which the coordination number (total number of contacts per particle) for each particle remains above 1. It is given by: where R is the particle's radius, ρ its density, G the shear modulus and ν the Poisson's ratio. This formula assumes that the relative velocity between contacting particles is very small. Other than for quasi-static systems, in practice some fraction of this maximum value is used, and for high coordination numbers (4 and above) a typical time step of 20% ∆t max has been shown to be appropriate. For lower coordination numbers, 40% ∆tmax is more suitable [47]. Choosing a suitable time step represents a compromise between the computational complexity/time, calculation error and simulation stability. Generally, the time step ranges between 20% and 80% of the critical time step according to Rayleigh. With regard to the duration of the calculation and its accuracy, a value of about 30% ∆t max was chosen for all simulations. In the proposed experiment, the influence of the shear modulus value on the process of discharging a hopper was investigated. A virtual cylindrical hopper with a conical discharge hopper of a 50 mm diameter was modelled. The conical hopper of a height of 100 mm was connected to the cylindrical part of a diameter of 100 mm. The experiment was performed for 15,000 pieces of D6 particles generated into the hopper. Subsequently, the flow rate of particles from the discharge hopper and the simulation calculation time were monitored for 5 s.

Static and Dynamic Angle of Repose Test
The process of forming a pile of particles is important in all industries where particulate materials are used. The angle of repose represents one of the basic properties of particulate matter, which affects transport and storage. In contrast to liquids or solids, granular materials create an angle of repose that depends on many of the material's properties. Parameters such as internal friction between particles, particle size and shape, size distribution, surface, moisture and electrostatic or cohesive forces may influence the value of the natural angle of repose. The angle of repose (static: S-AoR, and dynamic: D-AoR) is also one of the most important parameters for characterizing the flowability of particulate materials. The dynamic angle of repose represents the angle of slope at which a given flowing cohesionless material comes to stable dynamic state. Calibrations based on the angle of repose are generally performed by comparing experiments and simulations on a 1:1 scale. The interaction coefficients can be gradually optimized by systematically changing the input parameters or by simulations managed by the optimization algorithm [48,49]. At the end of the calibration, a set of DEM parameters is determined that provides the closest simulation results to those of the experiment. In particular, rolling and static friction are parameters which affect the macroscopic behavior of cohesionless bulk materials.
AoR of particulate material is influenced by many factors, such as the static and rolling frictions of particles, particle sizes and shapes, the quantity of material used during measurement and the measurement method (see Figure 5a). It is necessary to choose a particular measurement method based on the application, since each method provides slightly different AoR values. Particular materials were chosen for the experiments with regard to their grain shape so it would be possible to evaluate the advantages and disadvantages of the individual methods. The "piling" and "rotating" methods (see Figure 5a) proved to be the most universal and most credible for the determination of the AoRs of particular materials. If simple conditions are observed, these methods are able to provide high quality AoR values for a wide range of bulk materials. Angle of repose has been widely discussed in the literature for many decades and has been researched by many scientists, even outside of the process engineering field. This issue is first to be divided into two parts: (1) a study of how a pile of particulate material is formed; and (2) how the angle of repose is determined. Piling, pouring, emptying and slowly rotating in a drum leads to the determination of the static angle of repose. The dynamic angle of repose may be obtained in a rotary drum at higher rotational speeds [50]. Each of these methods in a way influences the resulting angle of repose. of whether a particular material could be retained with existing properties or if any of the parameters needed to be adjusted. Each of methods was evaluated in two ways, both graphically and geometrically (see Figure 6b).

Hopper Discharge Test
The calibration of the virtual material based on the discharge time of the hopper is important not only in terms of setting adequate flow properties, but also the rates of dispersion and diffusion of particles in the volume, which are crucial for all processes involving mixing. For all calibration experiments, it is important to find a suitable compromise between the number of particles used in the calibration test and the calculation time of single experiment. The more particles in one experiment that are calibrated, the more accurate the results which can be achieved. However, the experiments must be well-performed (simple movements of the apparatus, particles must not cause changes in the position of the apparatus, etc.) and reproducible over time. The hopper model was used to calibrate free-flowing materials. After filling the hopper with a particulate material and opening the discharge orifice, the material starts to flow onto a slide and slides/bounces to the bottom part of the model, where it forms a slope of a natural angle of repose. In addition to subjective observation and comparison of the course of the experiment with the simulation, this method makes it possible to calibrate the material using several different quantifiable aspects. For example, the particle speed at various points, particle tracking, the angle of repose defined after the experiment, the flow through the discharge orifice, the discharge time of the hopper and the mass loss Calibration using the static angle of repose is primarily used to optimize interparticle interaction parameters. For the calibration via the dynamic angle of repose, a rotating drum was used with the static friction between particles and a PMMA cylinder. A rotating drum was also used for an experimental study of homogenization. The drum diameter of 140 mm corresponds to a number of about 25 particles across the surface profile. An adequate number of particles is important for determining the dynamic angle of repose. A low number of particles may lead to distorted results. Similarly, if there are too many particles in the drum, the DEM calculation times increase significantly. As some materials travel in different ways in the drum, their trajectories were in all cases compared with the recordings of the experiment. This comparison allowed for the evaluation of whether a particular material could be retained with existing properties or if any of the parameters needed to be adjusted. Each of methods was evaluated in two ways, both graphically and geometrically (see Figure 6b).

Hopper Discharge Test
The calibration of the virtual material based on the discharge time of the hopper is important not only in terms of setting adequate flow properties, but also the rates of dispersion and diffusion of particles in the volume, which are crucial for all processes involving mixing. For all calibration experiments, it is important to find a suitable compromise between the number of particles used in the calibration test and the calculation time of single experiment. The more particles in one experiment that are calibrated, the more accurate the results which can be achieved. However, the experiments must be well-performed (simple movements of the apparatus, particles must not cause

Hopper Discharge Test
The calibration of the virtual material based on the discharge time of the hopper is important not only in terms of setting adequate flow properties, but also the rates of dispersion and diffusion of particles in the volume, which are crucial for all processes involving mixing. For all calibration experiments, it is important to find a suitable compromise between the number of particles used in the calibration test and the calculation time of single experiment. The more particles in one experiment that are calibrated, the more accurate the results which can be achieved. However, the experiments must be well-performed (simple movements of the apparatus, particles must not cause changes in the position of the apparatus, etc.) and reproducible over time. The hopper model was used to calibrate free-flowing materials. After filling the hopper with a particulate material and opening the discharge orifice, the material starts to flow onto a slide and slides/bounces to the bottom part of the model, where it forms a slope of a natural angle of repose. In addition to subjective observation and comparison of the course of the experiment with the simulation, this method makes it possible to calibrate the material using several different quantifiable aspects. For example, the particle speed at various points, particle tracking, the angle of repose defined after the experiment, the flow through the discharge orifice, the discharge time of the hopper and the mass loss characteristics of the hopper over time.
Another option would be to use PIV analysis and compare the shape of vector fields with quantifiable velocity results. The inclination of the walls of the hopper and the slide is variable and the size of the discharge orifice and can be adjusted according to the material being calibrated.
In this study, calibration was based on the discharge time of the hopper, which contained a given number of particles. Furthermore, a visual assessment of the material flow from the hopper and through the slide was performed. First, an experiment was performed for each particle type. The hopper inclination was 60 • in all cases while the discharge orifice width was set for each material so that the particles could flow freely. The whole experiment was recorded at a frame rate of 240 fps. Therefore, the exported record could to be slowed down to 12% of the real time of the experiment. This made it possible to determine the time period from opening the discharge orifice to the point when the last particle left the hopper with an accuracy of 0.004 s. DEM simulations were saved in steps of 0.01 s. Interaction parameters were then optimized and validated for the static and dynamic angles of repose.

Particles' Shapes and the Packing Test
A series of simulations was performed during the course of the packing test to initially estimate the friction parameters of the contact model. For spherical particles D6, cubic particles C6_24s and cylindrical particles P6 × 15_5s, 56 simulations of the packing test described in Section 2.2.2 were performed for each particle. The simulations involved combinations of static and rolling friction in the range of µ s-pp = 0.05, 1.5 and µ r-pp = 0.005, 1.5 . Tables 4-6 show the results of percentage-wise deviations from real experiments, in which the particle bed height for a given number of particles was 184 ± 1.3 mm for D6 particles, 202 ± 1.8 mm for C6 particles and 190 ± 1.6 mm for P6 × 15 particles (as shown in Table 3.). Thus, friction parameters can be gradually optimized by systematically changing input parameters or by simulations managed by the optimization algorithm. By gradual refinement of the friction scale values at the points of best match, more accurate estimates of the input parameters can be achieved. These frictional parameter estimates can then be combined with experimental static friction measurement to optimize the rolling friction coefficient, which is experimentally more difficult to determine. This can also solve the situation where multiple combinations provide an equal match to the experiment; e.g., µ s-pp = 0.1 and µ r-pp = 0.5; and µ s-pp = 0.2 and µ r-pp = 0.5 for D6 particles (see Table 4). The results of the calibration experiment involving volume filling are shown in Table 7. Based on Tables 4-6, initial friction parameters were used as follows: µ s-pp = 0.5 and µ r-pp = 0.1 for cubical and cylindrical particles C6 and P6 × 15; and µ s-pp = 0.2 and µ r-pp = 0.01 for spherical particles D4, D6, D12, D6 ABS and D5 steel . This table demonstrates that only virtual particles D6 ABS , D5 steel and C6 wooden cubes achieved an acceptable deviation of up to 5%, but C6 particles required much longer computational time. Therefore, it is necessary to optimize the size of all other virtual particles so that it better corresponds to the real experiment. The C6_24s particle was chosen as the virtual particle of the wooden cube for this study due to its geometric representation of the actual particle and due to the acceptable computation time per particle. This ratio was used like one of the criteria during the particle clump selection. The virtual particles D4, D6, D12 and C6_24s were subsequently enlarged by 0.2 mm in their characteristic dimension. Particles D6 ABS and D5 steel were kept in their original state. Virtual particle P6 × 15_5s was chosen for simulation of 15 × 6 mm pellets. This was also unchanged due to the very good result of the calibration experiment with a 0.5% deviation. The results of the optimized particles already showed very good agreement with the experimental data. The bulk density calibration for complex polydisperse particulate materials could be solved in a similar way which would allow for the determination of the best correlations between particle size distribution, density, friction parameters and particle bed height. The mass of the particle bed has to conserve the real mass of material in the case of particle size distribution (PSD) and/or particle size simplification by particle density setting.  Table 8 shows the coefficient of static friction of particles for different contact materials. The methodology described in Section 2.2.1 was used. In all cases, the standard deviation was based on ten repetitions of individual experiments. The measurements were influenced mainly by the quality of the surfaces of contact materials, the orientation of the particle surfaces and the speed of arm tilting. The measured values reflect the actual physical properties of the particles and could be used for the contact model, but it may be necessary to adjust them in conjunction with calibration experiments. Shear friction depends primarily on the quality of the contact surfaces, which also represents a problem in practice. The movement of particulate matter along the contact surfaces causes their degradation or even polishing. For this reason, the coefficient of static friction between particles and contact surfaces may change over time. In the case of wooden beads and their static friction, it is important whether a bead hole is placed on the wall material surface or not. The results can also be affected by the processing and finish quality of individual beads or the direction of tooling. In plastic particles placed on plexiglass, the friction can also be affected by electrostatic charge, etc. Experimentally measured values of static friction show the predicted range from 0.2 to 0.5. The same methodology was used to measure the mutual static friction between the particles, only the wall material sample was replaced with the surface of another particle. The static friction coefficients of individual particle combinations are shown in Table 9.  Table 9. Results of particle-particle static friction coefficients values. The results of rolling friction measurements shown in Table 10 are based on the same measurement methodology as that described in the previous Section 2.2.1, except that the entire particle was used for the measurement instead of the particle surface. Except for C6 particles, all other particles had a spherical shape. They had a much smaller rolling friction coefficient than the static coefficient. Therefore, the tilt of the arm on the apparatus can be used to determine the moment when a particle starts to roll and determine the rolling friction coefficient using the Equation (4). For particle C6, the default value µ r-pw = 0.5 was chosen based on the packing test due to difficult/impossible experimental determination. The particles were tested on a steel wall material, and one global particle-wall rolling friction coefficient was used for each particle for the next calibration experiments. In determining the coefficient of rolling friction, it is essential that the axis of the cylindrical particle remains perpendicular to the direction of motion throughout the experiment. The values of the coefficient of restitution (CoR) between individual particles and contact materials are shown in Table 11. The drop test methodology chosen and often discussed in the literature is very well applicable to spherical particles in which the trajectory of drop and reflection takes place in one plane. However, a problem arises with non-spherical particles, in which the reflection trajectory depends on the position of the particle upon collision with the contact material. If the material is reflected off the recorded plane, the velocity vector decomposes into individual components, thereby distorting the resulting value of the coefficient of restitution. As a result, the Equation (5) becomes inaccurate, and it is necessary to evaluate the reflection based on the speed in all directions. Furthermore, in this experiment, it is necessary to ensure sufficient stiffness of the contact material sample to prevent reflection damping caused by its deformation. The pendulum method with one hinge with the test particle and a vertically placed wall sample was used for C6 particle coefficient of restitution determination. This test used a method shown in Figure 3c, but particle B was replaced by fixed wall sample. CoR was determined using the Equation (5) in the same way as in the case of a two-particle collision. The double pendulum method described in Section 2.2.1 (see Figure 3c) was used to measure the coefficient of restitution of two particles. Table 12 shows the coefficients of restitution for particle interactions. The values of the coefficient of restitution were determined according to the heights and the angles of drop and rebound. These measurements were performed only for pairs of the same particles due to their identical weights. Other values of the coefficient of restitution for mutual combinations of different particles were determined according to the following equation [51]:

Shear Modulus Value Optimization
Following the conclusions made by Chen et al. (2017) in their study on the influence of Young's modulus on the process of mixing in a rotary cylinder, and the knowledge gained during the experiments with hopper discharging, it can be stated that the reductions of Young's modulus and the shear modulus, respectively, do not affect the courses of these processes in which the particles interact at low speeds. Figure 6 and Table 13 demonstrate that even a shear modulus change in the range of 10 4 Pa did not affect the flow properties of the material during the process of discharging. Different particles shapes achieved very similar results with slight differences of shear modulus of between 1 × 10 6 and 1 × 10 10 Pa. Bin discharge simulations exhibited very good reproducibility of results. The number of emptied particle differences was 0.5% maximally during simulation repetitions. Conversely, the calculation time increases exponentially as the value of the shear module increases. With suitable optimization, qualitatively comparable results can be achieved in much less time. In complex models with large amounts of particles, this could lead to potential time savings of several days. The optimal value of 0.001 times of the actual values was chosen (i.e., they are converted from gigapascals to megapascals).

Angle of Repose Calibration Tests
At first, a comparative measurement of four methods of determination of the angle of repose for particulate materials was performed. Each of these four methods (see Figure 5a) was evaluated in two ways, both graphically and geometrically (see Figure 5b). For each of the three selected materials (spherical D6, cubical C6 and cylindrical P6 × 15) ten measurements were made for each measurement method. A comparison of the resulting AoR values is shown in the graph in Figure 7a. The results show the differences between the methods, not only in the form of a pile but also in the method of measurement evaluation. Figure 7b shows experimental results of different pile formation methods. All materials were first tested with a "piling" method. Here the measurement is influenced mainly by the ratio of the particle size to the diameter of the bowl onto which the material is piled. Furthermore, the measurement is influenced by the height at which the material is dosed, as at too great a height the pile is disintegrated by the falling grains. For the "pouring" method, the shape of pile is informed by the speed at which the cylinder is lifted, and the friction between the cylinder wall and the bulk material. If the cylinder is lifted too fast, the material does not spill out of the cylinder immediately due to friction and inertia, but rather with a delay when the cylinder is at a certain height above the bowl, which results in a larger "spill" of the material on the surface. For the "emptying" method, care must be taken to ensure that the size of the drain hole is sufficient to avoid particle bridging, while in the "rotation" method a sufficiently low speed allowing static AoR must be observed. Ten piles were created for each material, and every pile was photographed from four views. Autodesk AutoCAD was used for manual repose angle determination. Two angles of repose were determined per picture (both sides of pile). All of values obtained were averaged, including relative standard deviation (RSD) determination.  Figure 7a,b show the results of all experiments with wood particles of different shapes. Some of the methods prove to be unsuitable for certain shapes and sizes. For example, for the "piling" method, the bowl is too small for the P6 × 15 particles and the particles are unable to form a proper cone to measure the angle of inclination. The evaluation of this method by an operator is then a highly subjective. Moreover, the "emptying" method proved to inappropriate for C6 particles, since we observed bridging above drain hole or the strong core flow, when only the particles above the drain hole were dropped, creating a slope angle of almost 90°. Figure 6c shows the tendency of the individual methods of pile formation towards higher or lower AoR values. The "pouring" method had the lowest AoR values, which were close to the dynamic AoR. Thus, the method can be used, for example, for the capacity calculations of belt All materials were first tested with a "piling" method. Here the measurement is influenced mainly by the ratio of the particle size to the diameter of the bowl onto which the material is piled. Furthermore, the measurement is influenced by the height at which the material is dosed, as at too great a height the pile is disintegrated by the falling grains. For the "pouring" method, the shape of pile is informed by the speed at which the cylinder is lifted, and the friction between the cylinder wall and the bulk material. If the cylinder is lifted too fast, the material does not spill out of the cylinder immediately due to friction and inertia, but rather with a delay when the cylinder is at a certain height above the bowl, which results in a larger "spill" of the material on the surface. For the "emptying" method, care must be taken to ensure that the size of the drain hole is sufficient to avoid particle bridging, while in the "rotation" method a sufficiently low speed allowing static AoR must be observed. Ten piles were created for each material, and every pile was photographed from four views. Autodesk AutoCAD was used for manual repose angle determination. Two angles of repose were determined per picture (both sides of pile). All of values obtained were averaged, including relative standard deviation (RSD) determination. Figure 7a,b show the results of all experiments with wood particles of different shapes. Some of the methods prove to be unsuitable for certain shapes and sizes. For example, for the "piling" method, the bowl is too small for the P6 × 15 particles and the particles are unable to form a proper cone to measure the angle of inclination. The evaluation of this method by an operator is then a highly subjective. Moreover, the "emptying" method proved to inappropriate for C6 particles, since we observed bridging above drain hole or the strong core flow, when only the particles above the drain hole were dropped, creating a slope angle of almost 90 • . Figure 7a shows the tendency of the individual methods of pile formation towards higher or lower AoR values. The "pouring" method had the lowest AoR values, which were close to the dynamic AoR. Thus, the method can be used, for example, for the capacity calculations of belt conveyors, etc. On the other hand, the "emptying" method showed the highest AoR values and a large deviation between the measurements. The rotating drum method showed very stable and consistent results. We also observed differences between the two methods of evaluation of experiments, as the geometric evaluation showed lower AoR values. Based on introductory AoR experiments, pilling and rotating methods have been selected for the DEM calibration process. These two methods exhibit the most stable results and high usability for different particle shapes. In some cases, the calibration using piling method with a small bowl and a small number of particles is more difficult than the calibration of a big pile on the free surface. Performing both ways together with a rotating drum experiment is the optimal situation. Graphical evaluation of the piling method gives closer results to the rotation drum method, so it was used for AoR evaluation in all cases. However, it is advisable to observe the overall shape of the pile; e.g., the roundness of the pile tip.
The results of calibration experiments for the determination of the static angle of repose are shown in Figure 8 and Tables 14 and 15. Materials D4, D6, D12, P6 × 15 and D6 ABS kept their experimentally measured input parameters. For material C6, the value of the rolling friction coefficient was adjusted from 0.50 to 0.15; for material P6 × 15 the value of rolling friction was changed from 0.06 to 0.03, while the value of static friction between particles was changed from 0.40 to 0.25. Similarly, for steel spheres D5 steel , the value of rolling friction was adjusted from 0.02 to 0.03, and the value of static friction between particles from 0.17 to 0. 30. shown in Figure 8 and Tables 14 and 15. Materials D4, D6, D12, P6 × 15 and D6ABS kept their experimentally measured input parameters. For material C6, the value of the rolling friction coefficient was adjusted from 0.50 to 0.15; for material P6 × 15 the value of rolling friction was changed from 0.06 to 0.03, while the value of static friction between particles was changed from 0.40 to 0.25. Similarly, for steel spheres D5steel, the value of rolling friction was adjusted from 0.02 to 0.03, and the value of static friction between particles from 0.17 to 0.30.

Particle
Piling (see Figure 7) Base 90 mm in Diameter Piling (see Figure 7) Base 150 mm in Diameter Rotating Drum (see Figure 8) Figure 9 shows a comparison of the behaviors of all particulate materials at a drum speed of 6 rpm. For material D6, the static friction coefficient between the particles and the PMMA cylinder wall was adjusted from 0.54 to 0.35. The behaviors of the other materials showed satisfactory results with the original experimentally measured and calibrated parameters. Ten measurement repetitions were made for each experiment. Based on good reproducibility of DEM simulations, the first calibration simulation was performed for the AoR determination and the second simulation for agreement of the results of every particle series. The final AoR values of calibrated materials are shown in Tables 14  and 15.

Particle
Piling (see Figure 7) Base 90 mm in Diameter Piling (see Figure 7) Base 150 mm in Diameter Rotating drum (see Figure  8)

Particle
Piling (see Figure 7) Base 90 mm in Diameter Piling (see Figure 7) Base 150 mm in Diameter Rotating Drum (see Figure 8)

Discharge Test
The discharge test was performed according to the methodology described in Section 2.2.5. The results of the calibration tests based on the discharge time of the hopper are shown in Figure 10 and summarized in Table 16. Figure 10 shows the different phases of the discharge of the hopper model. Comparisons of the experiment and the DEM model for the individual particles were captured at the same time; i.e., at the same time from the release. The discharge test represents a rather versatile calibration tool. It allows one to calibrate the flow properties given by the friction parameters, and the slip and bouncing of the particle flow from the slide, and grants the option to use various contact materials. Another way to evaluate the discharge test is to monitor the movement of particle layers of different colors or to evaluate velocity fields using PIV analysis. In this study, the discharge times of the specified numbers of particles were used for comparisons. Furthermore, the recordings of experiments and DEM simulations were visually compared. Some of the input parameters of the DEM contact model were further refined by this calibration test and back-tested.

Discussion
The simple shapes of some products can be easily achieved by the "clumping" method (e.g., agricultural commodities such as rice, cereals and pellets). The higher the number of spherical surfaces, the higher the accuracy that can be obtained when approaching real morphology. However, the computational power required to calculate the contact forces between particles with a large number of surfaces increases significantly.
The first calibration experiment involved filling the volume of a cylinder (the so-called "packing test") to determine the shape accuracy and bulk density of each virtual material. The series of simulations was compared with real experiments, and the sizes, shapes and densities of virtual particles were optimized. The computational costs of a particle made up of more than one spherical surface were recorded. It represented one of the parameters that played a role in the selection of the creation method of the resulting cubic particle. Interestingly, cubic particles C6 formed from 24 spherical surfaces showed less computational complexity than particles C6 with 16 spherical surfaces.

Discussion
The simple shapes of some products can be easily achieved by the "clumping" method (e.g., agricultural commodities such as rice, cereals and pellets). The higher the number of spherical surfaces, the higher the accuracy that can be obtained when approaching real morphology. However, the computational power required to calculate the contact forces between particles with a large number of surfaces increases significantly.
The first calibration experiment involved filling the volume of a cylinder (the so-called "packing test") to determine the shape accuracy and bulk density of each virtual material. The series of simulations was compared with real experiments, and the sizes, shapes and densities of virtual particles were optimized. The computational costs of a particle made up of more than one spherical surface were recorded. It represented one of the parameters that played a role in the selection of the creation method of the resulting cubic particle. Interestingly, cubic particles C6 formed from 24 spherical surfaces showed less computational complexity than particles C6 with 16 spherical surfaces. This was due to the larger particle diameters and the longer computation time step. However, this parameter had only minimal influence in the case of pellet particle P6 × 15.
Three test apparatuses were designed and assembled for the following series of calibration experiments. Namely, an apparatus for measuring the interaction coefficients (i.e., the static friction coefficient, rolling friction coefficient and coefficient of restitution); a rotary drum apparatus for measuring the dynamic angle of repose and flow properties; and a hopper model for dynamic flow monitoring, which allows for the calibration of virtual materials according to the flow method (mass vs. core), angle of repose, discharge time, etc. Using these three apparatuses, the input parameter values were experimentally determined for a contact model that defines the behavior of particulate materials in DEM simulations.
The influence of the shear modulus value was investigated during the calibration of the virtual material. This part of the paper showed that the use of a shear modulus value three orders of magnitude lower (from GPa to MPa) than the tabulated values for a given material does not affect the material behavior in the simulations. Shear modulus values do not affect the flows of differently shaped particles. Indeed, this simplification has a significant impact on the simulation computation time, which can be reduced by up to 35 times. Furthermore, the interaction parameters were calibrated and optimized based on the static and dynamic angle of repose or discharge time of the hopper model. The AoR of a particulate material is influenced by many factors, such as the static and rolling friction of particles, particle sizes and shapes, quantity of material used during measurement and the measurement method. This study performed comparative measurements of differently shaped materials using different methods of pile formation. The results show that it is necessary to choose a particular measurement method based on the application, since each method provides slightly different AoR values. Particular materials were chosen for the experiments with regard to their particle shape, so that it would be possible to evaluate the advantages and disadvantages of the individual methods. The "piling" and "rotating" methods proved to be the most universal and most credible for determination of AoR of particular materials. If simple conditions are observed, these methods are able to provide high quality AoR values for a wide range of bulk materials.
These practical tests helped to refine some of the interaction coefficients, which are difficult to determine experimentally. Thus, an optimal relation between the set of all interaction coefficients was found, since the mutual combination of these values defines the behavior of a particulate material. However, even if the individual coefficients are determined relatively accurately, their combination in the contact model used may not reflect the actual behavior of the particular material due to some further simplifications. All calibrated DEM parameters of investigated particulates are mentioned in Tables A1-A7 in Appendix A.

1.
The packing test can be used for initial estimation and to optimize and find a suitable combination of the input parameters for the DEM contact model. Furthermore, the packing test is a suitable tool for optimizing the shape representation of particles and bulk density.

2.
For the Hertz-Mindlin contact model, the coefficient of static friction and the coefficient of restitution can be relatively easily determined experimentally. In addition to these values, it is necessary to find a suitable rolling friction coefficient to ensure both the static and dynamic behavior of the material.

3.
By optimizing the shear modulus, the computational time of the simulation can be significantly reduced. Even a shear modulus reduction by 10 4 Pa does not affect the flow properties of the material.

4.
Calibration using a static AoR is generally the most commonly used calibration experiment. However, the results are influenced by many factors, such as the number of particles, the forming method of the pile and the method of evaluating the final value of the AoR. During selection of a suitable calibration methodology, it is necessary to base it on the simulated process and material properties. 5.
The most reproducible results of the angle of repose were achieved by the "pilling" method and in the rotating drum when evaluated by the geometric method. However, it is always advisable to make an overall visual comparison of the slope shape between the calibration simulation and the experimental curves. The base to particle size ratio should be greater than 25, and the calibration experiment should contain approximately 4000 particles to ensure representative and reproducible results.

6.
Dynamic behavior must be verified in the same way as static. Discharge test is a suitable tool for assessing flow properties. One set of input parameters must then best characterize the static and dynamic behavior of the material.

Conflicts of Interest:
The authors declare no conflict of interest.  Table A3. Calibrated coefficients of restitution values (CoR particle-particle).  Table A5. Coefficient values of contact static friction (particle-particle).