DEM Calibration Approach: design of experiment

The problem of DEM models calibration is considered in the article. It is proposed to divide models input parameters into those that require iterative calibration and those that are recommended to measure directly. A new method for model calibration based on the design of the experiment for iteratively calibrated parameters is proposed. The experiment is conducted using a specially designed stand. The results are processed with technical vision algorithms. Approximating functions are obtained and the error of the implemented software and hardware complex is estimated. The prospects of the obtained results are discussed.


Introduction
The discrete element method (DEM) is the most popular numerical method for bulk materials handling simulation [1][2]23]. The corresponding software with a graphical user interface is an effective tool for optimizing the equipment of the mining and metallurgical industry. DEM is now widely used in integration with Computational Fluid Dynamics (CFD) and Finite Element Analysis (FEA), and allows one to simulate complex multiphase and multidiscipline processes that were not previously covered by numerical methods [3][4][5].
Before DEM simulation starts, it is necessary to input parameters which determine the bulk material behavior: particle shape and size distribution, density, static and dynamic friction coefficients, coefficient of restitution, Young's modulus, Poisson's ratio, and rolling resistance if modeling spherical particles. While some of the input parameters can be acquired from technology data (e.g. .Young's modulus, particle size distribution), the remaining parameters have to be calibrated [6]. The same input parameters of various contact models used in the DEM codes can affect the behavior of the bulk material differently [7].
The whole set of existing calibration methods can be divided into two groups: Bulk Calibration Approach and Direct Measuring Approach [7]. In the first case, a laboratory-scale experiment is performed with the bulk material, which is then simulated using DEM. DEM input parameters iteratively changes until results become similar to the actual experiment's outputs [8][9][10]. This approach is time consuming when using simple parameters search methodology and is dependent on application: suitable parameters set for simulation of one process can be inaccurate for another. Also in this case the physical meaning of the input parameters is lost. The second approach is based on direct measurement of input parameters values separately using unified particles of small quantity. Corresponding material tests (e.g. share test) are carried out using well-known measuring techniques [11][12][13]. This method will be accurate only if such parameters as particle size distribution or particle shape are perfectly reproduced and when using an approved and adopted DEM contact model. This approach is time consuming and needs to be done using a certified expensive set of laboratory testers. Bulk Calibration Approach is often used with Direct Measuring Approach [14]. For example, size distribution and share modulus can be measured directly, and the friction coefficients are obtained. This allows reducing the number of variables by reducing combinations of parameters that lead to the same behavior of the bulk material, and the time for calibration respectively. Using measurable macroscopic parameters (e.g. the angle of repose or the flow time) for the calibration process, let us get adequate behavior of bulk solids relatively fast. The particle shape can be simplified to a sphere, and the size distribution could be replaced by an equivalent diameter [15]. In this case, the rolling resistance coefficient has to be calibrated [16].
While it is impossible to get rid of the loss of physical sense of the calibrated parameters, the problem of having various combinations of input parameters is completely solvable. Thus, Derakhshani et al. [17] for the calibration of the coefficient of sliding friction and rolling friction used only two macroscopic parameters -the angle of repose and the time of the complete outflow of material from the funnel. Using two parameters only, the unique combination of DEM input data was found.
One of the task of this study is to classify DEM input parameters into two groups: 1) parameters that should be measured directly, 2) parameters that are easier to obtain iteratively. Further work assumes learning how to describe the bulk material with a set of unique macroscopic parameters that can be converted into a single combination of the DEM input parameters.

Introducing the idea
The first step is to consider the input parameters typical for any DEM software. Let us divide it into four categories: related to Boundary, to Particle, Paired and Common. Common parameters include acceleration of gravity (usually 9.81 m/s 2 ), as well as the type of the contact model of interaction (linear or non-linear). Boundary parameters include the density of the boundary's material, the Poisson ration and Young's modulus. Particle parameters include Poisson's ratio, Young's modulus, bulk density, particle shape, particle size distribution and rolling resistance of spheres particles. Paired parameters include paired particle-particle (SFp-p) and particle-boundary (SFp-b) static friction coefficients, particle-particle (DFp-p) and particle-boundary (DFp-b) dynamic friction coefficients and restitution coefficients (also called damping coefficient) for particle-particle (CoRp-p) and particleboundary (CoRp-b) interaction. All parameters are presented in table 1.
Each of the presented parameters has a physical meaning. But the measurement of all Paired parameters and rolling resistance is technically difficult, time consuming and, as mentioned above, direct measurement does not guarantee the accurate results. The contact model is a mathematical description of particle interaction, but since a number of parameters will be chosen iteratively, the choice of the contact model will not have a significant effect on the result [18][19]. Other parameters are measured at the macroscopic level by well-known methods [20]. Based on this, it was decided to divide the parameters into those that are better measured directly and those that will be obtained iteratively (table 2).
The number of iteratively obtained parameters is seven when modeling spheres or six when modeling non-spherical particles. When values of directly measured parameters are known, bulk material behavior will depend on seven parameters. It means that any macroscopic parameter will be a function of the seven coefficients. The solution of the inverse task (the choice of coefficients over the measured macroscopic parameters) implies the existence of a functional relationship between the parameters and the coefficients. Therefore, in order to obtain a single combination of the seven modeling coefficients, it is necessary to measure seven unique characteristics of the bulk material (to have seven functional dependencies). It is worth considering that this expression is valid for the linear influence of the coefficients on the parameters. If the influence of higher orders is detected, the results may not be so unambiguous.    The purpose of analyzing the bulk material behavior is to identify seven unique parameters. The choice of parameters for the measurement is based on empirical data. Since the variable parameters carry the physical meaning of the dynamic and static characteristics of the bulk material, it was decided to analyze similar macroscopic parameters. To solve this problem, technical vision was applied. First of all, technical vision helps to automatize the experiments, significantly reducing the time for analysis. Secondly, frame-by-frame video processing allows one to obtain dynamic material characteristics (e.g. flow time) in addition to static (e.g. the angles of repose and rupture). The algorithms used in this study were implemented using the Python programming language and OpenCV image processing library [22].
The factorial experiment (FE) was used for evaluating the correlation of coefficients and measured parameters. Each of the measured parameters is affected by seven coefficients varying from 0.1 to 0.9. Real-word variation of the coefficients does not exceed 0.5-0.6 for most materials, but the shape and size of the particles are a simplification of real data and the parameters in the simulation are more the mathematical equivalents of their physical twin. After getting coefficients correlation and variation levels, experiment planning matrix 2 7 was designed, where the variable parameters are coded (lower level -1, upper +1) and the output data at the same level averaged (it is necessary to ensure repeatability). The dimensionless coefficients of the approximating function are calculated by: Then an approximating function is constructed, which in the general case has the following form: The obtained dependence is checked for adequacy with the F-test and coefficients at factors for signification with Student's t-test. After conducting the FE for each of the seven characteristics, the resulting 7 equations are combined into a single matrix. Taking into account that all equations are linear, the matrix solution gives one single combination of coefficients for the model.

Handling the experiment
It was decided to take the following input modeling parameters to carry out the experiment (table 3). Parameters for particles corresponded to fine crushed stone, and surfaces to plexiglas.
In equations set (4), 1 α and 2 α denote the angle of rupture and the angle of repose accordingly. In addition to them, 5 other equations are combined into the set. The solution of this set gives a single combination of input parameters for modeling a material with specified parameters. The next stage was the approbation of the results. Experiment was repeated for 10 points with different coefficients. The given coefficients and coefficients obtained by solving dependences are given in table 4. As can be seen from the results of the approbation, the developed approach realized on the basis of a specialized software and hardware complex allows calibrating the coefficients for arbitrary bulk material with an accuracy of 6%.

Conclusion
The proposed approach gives fairly accurate coefficients for modelling, and it is also important that the combination is unique. It should be noted that the original stand was designed for particles with a diameter of not more than 8 mm. To calibrate bulk materials with a large diameter, it is either needed to proportionally change the dimensions of the stand, or conduct a separate research with varying particle size. In addition, the approach involves carrying out 128 numerical experiments for the investigated material, which requires considerable computing power.
The obtained empirical data can be used not only within the framework of the proposed approach but also applied in iterative frameworks or learning algorithms (e.g. machine learning).