DETERMINING THE RELATIONSHIP BETWEEN THE SIMULATION DURATION BY THE DISCRETE ELEMENT METHOD AND THE COMPUTER SYSTEM TECHNICAL CHARACTERISTICS

The object of this study is the relationship between the technical characteristics of a computing system and the duration of modeling the motion of particles of granular materials by the discrete element method. The scheme of the calculation algorithm is presented; its main stages are analyzed. A 3D model of a belt feeder and a mathematical model of particle motion were developed for calculations by the discrete element method in the EDEM 2017 environment. The physical and mechanical properties of bulk material are defined; the structural and technological parameters of equipment are determined. The parameters of the algorithm and computing system are analyzed. Those parameters are defined, the change of which does not affect the accuracy of calculations but can change the volume of computing resources used. These include the number of particles of loose material, the «grid» step, and the number of processor cores. The influence of these parameters on the duration of simulation was determined using a complete factor experiment. Experimental studies have shown that for the duration of the simulation, the determining parameters are the number of particles and the number of processor cores. It was established that there is a linear relationship between the duration of the simulation and the number of particles. The regression equation is built, which makes it possible to predict the simulation time. It was also established that the software does not fully use all available computing resources; the maximum load on the processor when utilizing all available cores is 57 %. The use of RAM and disk subsystem almost did not change during simulation. The results reported here make it possible to plan the use of computing resources for research using the discrete element method and to predict the simulation time


Introduction
Software for simulating the movement of granular materials is widely used in various industries [1][2][3][4]. It makes it possible to predict the efficiency of technological equipment, determine its optimal parameters, take into account the influence of the physical and mechanical properties of particles on the nature of their movement. At the same time, mathematical methods do not make it possible to accurately model the movement of granular materials and their interaction with working bodies of the equipment. To solve these problems, various approaches were used, based on the mathematical models of motion of continuous media, rheological models, and others [4]. With the development of computing technology, the discrete element method (DEM) has become widespread [5]. It makes it possible to calculate the physical and mechanical interactions of individual particles between themselves and with equipment, that is, unlike other models, it takes into account the discrete nature of the movement of loose material. This provides an opportunity to improve the accuracy of modeling, in particular, to simulate gaps in particle flows and the appearance of vaults. This information can significantly reduce the design time of technological equipment for processing bulk materials used, in particular, in construction, light, food, mining, pharmaceutical, and other industries. The ability to simulate the interaction of solid particles (bulk medium) with working bodies of equipment significantly reduces the design time and optimization of their parameters.
One of the problems with using DEM is that improving accuracy is achieved by increasing the number of mathematical calculations and, accordingly, increasing the required computing resources. In practice, the duration of simulation by this method, even with the use of modern computer equipment, can take weeks or months [6][7][8]. As a result, various approaches have emerged to reduce the number of calculations. Thus, a common practice that reduces the total duration of calculations is to conduct preliminary simulations with a minimum number of particles. This makes it possible to relatively quickly check the parameters of the model, the nature of the interaction of particles with working bodies of the equipment, and then carry out modeling with the full number of particles. However, despite the possible options for optimizing the number of calculations, the total amount of required computing resources remains high, which makes the task of their effective use urgent.
Cloud services make it possible to significantly increase the available amount of resources. At the same time, to optimize the cost of their use, information on the actual amount of resources to be used in the calculations is needed.
These problems render relevance to the task to forecast the cost of computing resources in modeling by the discrete element method.

Literature review and problem statement
In [9], the authors reported the results of research on the phenomenon of superlinear acceleration in modeling the motion of particles of complex shape by DEM. The use of RAM, the efficiency of the processor cache, and the influence of the number of computing nodes are determined. It is shown that the use of parallel computing in the modeling by DEM can significantly reduce the time of calculations. It should be noted that the studies were conducted using super computers. This does not make it possible to conclude on the speed of simulation when using existing software for modeling by DEM on single-processor multi-core systems.
Study [10] proposes an algorithm for parallelizing the calculations of the motion of particles of complex shape based on DEM. Optimization of a number of parameters of the algorithm that affect the speed of parallel computing is performed. In particular, the cost of data transmission was mini mized, dynamic load balancing was maintained, redundant information about the contacts between the particles was removed. The proposed algorithm was checked for three-axis ellipsoidal particles using super computers. It has been established that the simulation time is a descending function of the number of computational nodes. At the same time, the results of the study do not show how much the characteristics of the proposed parallelization algorithm correlate with the characteristics of the algorithms implemented in the software products common today.
The results of the study of strategies for parallelizing calculations using DEM are reported in [11]. The authors consider two methods: OpenMP (parallel processes use a common area of RAM) and MPI (each process works with its own area of RAM). It is shown that with an increase in the number of processor cores and the number of particles, the time for data transfer between processes increases, which makes OpenMP 50-90 % faster than MPI. The results can be applied in the design of new software, but they are difficult to use to predict the simulation time using existing systems.
In [12], it is proposed to increase the speed of calculations by DEM by dividing the load between the central and graphics processor. It is shown that the performance of the proposed method is better compared to the methods of parallelization based on the central processor. It should be noted that the approach reported in the work imposes significant restrictions on the hardware of the computing system.
Paper [13] reports the results of using a hybrid DEM-PBM approach to reduce the computational costs of modeling the mixing of powder materials. The proposed solution is based on the fact that PBM (population balance method) uses significantly fewer computing resources compared to DEM. The authors have determined the optimal duration of DEM modeling, which ensures the preservation of the accuracy of the mixing dynamics. They also defined predictive accuracy of PBM models obtained on the basis of DEM modeling. Verification of the results was carried out on the basis of comparison with the full modeling by DEM, which confirms the effectiveness of the proposed method. The method can significantly reduce the number of calculations but requires monitoring the accuracy of the simulation results.
The accuracy and computational costs of the LS-DEM technique are compared with DEM in [14]. The LS-DEM technique, compared to DEM, makes it possible to take into account particles of arbitrary shape by storing the distance values to the surface on the grid for each discrete element. Boundary nodes located on the surface of the discrete element are also taken into account. It is shown that using the OpenMP parallelization approach makes it possible to significantly reduce the time spent on modeling when using LS-DEM. It should be noted that the results relate to the LS-DEM technique, which does not allow them to be extended to DEM.
In [15], the results of software development for DEM based on the method of domain decomposition are presented. The effectiveness of the code was estimated by modeling the fluidized bed for 2.5 million particles. It is shown that the effectiveness of the proposed solution reaches 81 %. At the same time, the reported results do not make it possible to compare it with the implementations of DEM used in the software distributed today.
The results of the development of parallel algorithms for DEM software are presented in [16]. For clusters of personal computers with distributed memory, the authors have developed a new algorithm. It combines communication cells to detect contacts, static domain decomposition, and MPI data transfer for processors. It is shown that when using 16 processors, an acceleration of 11 was obtained. It should be noted that the results are difficult to use for prediction of the speed of the algorithm for single-processor multi-core computing systems.
The results of parallelization of the Lagrange-Euler code for DEM/CDF (discrete element method/computational hydrodynamics) are reported in [17]. The authors showed that their proposed parallelization strategy demonstrated good scalability when using up to 32 processors. An increase in the number of processors leads to an increase in the overhead of interprocessor data transmission. The results in the study were obtained for the algorithm for determining the movement of a layer of fluidized particles, which complicates their use for processes with a different nature of the movement of loose material.
Our review of the scientific literature makes it possible to identify two main approaches to reducing the time of modeling the movement of granular materials: 1. Parallelization of calculations between different computing devices.
2. Modification of the algorithms used by DEM in order to replace part of the operations with less resource-intensive ones.
The first approach causes the problem of ensuring effective interaction between processors. It is shown, in particular in [11], that the use of shared RAM (OpenMP) simplifies the solution to this problem. Thus, the use of single-processor systems can provide greater relative efficiency in the use of computing resources, although the absolute speed of multiprocessor systems will certainly be greater.
The second approach requires checking the accuracy of the results. Such checks are carried out for each algorithm separately, which significantly complicates its application in practice.
When designing equipment for working with granular materials, it is necessary to be able to change the structural and technological parameters within a wide range. This leads to a significant increase in the number of experiments and calculations. Therefore, it is important for researchers to be able to predict computational times when using widespread software and computing devices.

The aim and objectives of the study
The purpose of this research is to determine the relationship between the duration of modeling the operation of technological equipment by the discrete element method and the technical characteristics of the computing system. This will make it possible to predict the simulation time and optimize the use of computing resources.
To accomplish the aim, the following tasks have been set: -to build a model for research using DEM; -to analyze the parameters of the mathematical model and technical characteristics of the computing system and determine those that can be changed without reducing the accuracy of DEM calculations; -to conduct experimental studies to quantify the influence of the parameters of the mathematical model and the computing system on the duration of DEM modeling.

The study materials and methods
The object of this study is the relationship between the technical characteristics of the computing system and the duration of modeling the motion of particles of granular materials by the discrete element method. The main hypothesis assumes that determining this relationship will make it possible to calculate the time of modeling by DEM and make informed decisions on the use of computing resources when making calculations using DEM.
DEM belongs to the group of numerical methods used to study large groups of particles: molecules, granules, grains of sand, pebbles, and others [5]. The basis of DEM is the assumption that the material under study consists of discrete particles that may have different properties and shapes. Modern software, in order to reduce the volume of calculations, makes it possible to simulate particles of complex shape as a set of spheres of different radii [18][19][20].
The scheme of the calculation algorithm is shown in Fig. 1. Before starting the simulation, we used Solidworks 3D modeling software (France) to build an equipment model that includes the housing, working bodies, and other elements with which the particles directly interact. This model in STL format was imported into the EDEM 2017 system (USA) [21], which makes it possible to make calculations using DEM.
The choice of EDEM 2017 among analogs (LIGGGHTS (Republic of Austria), EDEM (USA), LAMMPS (USA), YadeDEM (USA), and others) was primarily due to the ability to adjust the parameters of the computing algorithm and manage the use of computing resources, namely: -to set the method for determining the modeling step; -to determine the size of the cells of the «grid»; -to set t simulation time; -to set the time intervals at which data are stored in a permanent computer memory; -to determine the number of processor cores that will be involved in the calculations.

Fig. 1. Scheme of the algorithm for calculating by DEM
Before starting the calculation, the positions are set (x i ; y i ; z i ), initial linear (v xi ; v yi ; v zi ) and angular (ω xi ; ω yi ; ω zi ) particle velocities and physical and mechanical properties of their materials. At the third step of the algorithm, the simulation step Δt is determined. In the study, the Rayleigh criterion was used to calculate Δt [20], which is one of the most common in practice: where R, ρ, G, υ are the radius, density, shear modulus, and Poisson coefficient, respectively. In step 4, a grid is built. In this case, the modeling area is divided into cells of a given size. This operation minimizes the number of calculations since the calculations of the interaction between the particles are carried out only inside the cells. In the fifth step, the so-called «active» cells are determined, that is, those containing 2 or more particles. Between all particles inside the «active» cell, a check is carried out for possible interactions.
The sixth step of the algorithm involves the calculation of forces and moments acting on particles according to the following formulas [18]: where t is time; m i is the mass of the particle; I i is the moment of inertia of the particle; M ij is the moments of forces acting between the particles; F ij -forces acting on a particle on the side of other particles and objects; θ i -angular velocity of the particle, rad/s. Based on the obtained values, new particle positions are calculated, after which the calculation cycle is repeated starting from step 5.
As part of this study, the possible influence of the forces of electrostatic interaction between the particles was not taken into account.
We determined the amount of computing resources used during the simulation by using Windows system resource monitoring tools (USA).
During the research, all third-party applications were disabled, in particular, antivirus software, Windows update system, and others.
The calculation time was recorded in reports with the results of the simulation, which is generated by the EDEM 2017 system.
Mathematical processing of the results was carried out using a program written in the Python programming language [22], software Jupyter [23], and the NumPy library [24]. To calculate the coefficients of regression equations, the least squares method was used, to check the adequacy of the mathematical model -Fisher's criterion.

Results of investigating the duration of simulating
the operation of the belt feeder by the discrete element method

1. Model for research using the discrete element method
For the study of DEM, a belt feeder was chosen in this work. Equipment of this type is widely used in various industries for the transportation of granular materials [4]. Fig. 2 schematically shows the design and principle of operation of the belt feeder. The movement of the material is enabled by the friction forces between the particles and belt (3) stretched between the two drums. The granular material under the action of gravity from hopper (1) is poured onto the belt, which moves progressively. As a result, the particles move in the direction of movement of the belt. The thickness of the material layer is adjusted by means of shutter (2).

Fig. 2. Belt feeder design scheme
To study the simulation time of such a feeder, a mathematical model was built in the EDEM 2017 system, which included particles of loose material and two elements that interact with them: a hopper and a belt (Fig. 3).

Fig. 3. Belt feeder model built in EDEM 2017
The physical and mechanical parameters of particles, belt and hopper materials are selected based on the results of research [19,25,26] and are given in Table 1. The radius of the particles was 3 mm, the mass was 0.26·10 -3 kg, volume -1.13·10 -7 m 3 , moment of inertia -9.36·10 -10 kg·m 2 , the shape is spherical. The speed of movement of the feeder belt is 0.1 m/s. The values of the parameters of the interaction of objects with each other are given in Table 2. In this case, the parameters of the bulk material correspond to the granules of polymeric materials, the parameters of the belt -rubber, and the parameters of the hopper -steel [4,22].

2. Analysis of the parameters of mathematical mo dels and technical characteristics of the computing system
The parameters of the mathematical model and the settings of the DEM calculation algorithm include: 1) the physical and mechanical properties and shape of particles. The use of computational resources is directly affected by the shape of the particles. Particles of complex shape are modeled as a set of spheres of different radius. Increasing their number automatically increases the number of calculations in steps 5 and 6 of the algorithm (Fig. 1). At the same time, these parameters directly affect the trajectories of particle movement and should be set as close as possible to the physical properties of the object of study. Changing these parameters during preliminary studies is undesirable since it can lead to the formation of a misconception about the nature of the movement of the material; 2) the physical and mechanical properties of materials and the geometric parameters of equipment. These parameters are decisive for particle trajectories (simulation results). Reducing the number of computing resources is possible only by simplifying the design, which, at the same time, leads to a decrease in the accuracy of the simulation; 3) parameters of the interaction of objects with each other. The impact on the use of computational resources is minimal since these parameters determine the values of the coefficients in the equations used in step 6 of the algorithm (Fig. 1). At the same time, their change leads to significant changes in the simulation results; 4) modeling time. It determines the duration of calculations but does not affect the use of system resources at each individual point in time. It is set by the researcher based on the principle of operation of the equipment; 5) simulation step. It determines the number of iterations (repeated execution of steps 5-8) of the algorithm (Fig. 1), that is, it directly affects the volume of calculations. Reducing this value makes it possible to improve the accuracy of modeling, increasing -reduces accuracy and can lead to completely incorrect results. In most cases, the value is determined by one of the known criteria, for example, the Rayleigh criterion [20]; 6) the frequency of preservation of current results. The EDEM 2017 system makes it possible to set the frequency of full preservation of current calculation results. This parameter affects the load on the disk subsystem; 7) the «grid» step. It does not affect the result of calculations, only the duration of simulation. Reducing the step increases RAM usage. EDEM 2017 recommends using the «grid» pitch value be equal to 3 radii of the smallest particle; 8) the number of particles used in the simulation process. This parameter affects the intensity of interaction of particles with each other and directly affects the amount of computational resources since calculations for DEM are performed separately for each particle; 9) the number of processor cores used in the calculations. EDEM 2017 determines the available number of processor cores (in the case of support for hyperthreading technology, each physical core is considered as two). The user can determine the number of cores that the program uses and thus change the load on the system. Based on the results of the analysis, the following parameters were selected to study their impact on the duration of modeling (TM): 1. The number of particles (PN). 2. The «grid» step (CS).

Number of processor cores (CN).
The main selection criterion was the minimal impact on the accuracy of the simulation.

3. Experimental study of the influence of the param eters of the mathematical model and the computing system on the duration of DEM modeling
To determine the influence of the selected parameters, a complete factor experiment of type 2 3 was conducted. The simulation time for all experiments was 3 seconds. The parameter variation intervals are as follows: 1. PN = 2,000…3000 -selected based on the particle size and feeder so that during the simulation there is a supply of material in the bunker.
2. CS = 2R…4R (where R is the particle radius) -selected based on the recommendations of the EDEM 2017 system to the «grid» step.
3. CN = 1…8 -the minimum and maximum number of available processor cores.
The results of our study are given in Table 3. During each experiment, the following parameters were determined: 1) TM -duration of modeling, s; 2) PL -CPU usage by EDEM process, %; 3) RAM -RAM usage, %; 4) DL -load on the disk subsystem, %. The results show that the change in input parameters almost did not affect the parameters RAM and DL. This makes it possible to conclude that the use of RAM depends only on the complexity of the model. The load on the disk subsystem depends on the amount of data that needs to be stored (in all cases, this is a constant value).
For the parameters TM and PL, mathematical models are built in the following form:

TM B B CN B PN B CS B CN PN B PN CS B CN CS B CN
Models (3), (4) take into account the influence of both individual parameters on the original value and the mutual influence of the parameters.
The calculation of the coefficients A and B is carried out by the method of least squares; the result is the following values of the coefficients: . .
The coefficients A and B were checked for significance according to Student's criterion. As a result, it was established that the coefficients A2…A7 and B3...B7 are not significant. As a result, equations (5), (6) Verification of the obtained mathematical models was carried out according to the Fisher criterion, which confirmed their adequacy.
The results showed that the use of the CPU by the EDEM process depends solely on a given number of cores, and the duration of the simulation depends on the number of cores and the number of particles. Since in most cases researchers are interested in the highest possible speed of modeling, it is possible to take the CN parameter equal to the maximum value.
In order to identify the nature of the relationship between the simulation time and the number of particles, additional studies were carried out, during which the number of particles varied from 2000 to 3000 in increments of 200, CN = 8, CS = 4R. Simulation results confirm the linear nature of the relationship between the number of particles and the modeling time (Fig. 4).
Checking the model according to the Fisher criterion confirmed its adequacy.

Discussion of results of investigating the simulation time of the work of a tape feeder by the discrete element method
The algorithms used by DEM do not allow for a complete parallelization of the problem. This is due to the fact that the beginning of the calculations in step 5 (Fig. 1) should be performed only after determining the new positions of all particles. At the same time, the calculations in steps 5, 6, and 7 of the algorithm (Fig. 1) can be performed in parallel. Thus, the presence in the system of several physical computing devices (processors, processor cores, video cards) makes it possible to reduce the calculation time. In this case, the problem of uniform separation of the load on all available computing resources arises. DEM does not make it possible to predict the volume of calculations that will be necessary to determine the new positions of the particles. Since for this it is necessary to check the contacts between all particles, which in itself is one of the most resource-intensive steps [21].
The results of our study have made it possible to numerically assess the effectiveness of the use of system resources by the software product EDEM 2017. The study was conducted using a system with one processor and 4 physical cores. At the same time, the use of the processor varied from 13 % to 57 % depending on the number of cores used ( Table 3). The value of this range confirms that the calculations are distributed between the processor cores. At the same time, the maximum value (57 %) indicates incomplete use of available computing resources.
In addition, the results showed that the loads on the disk subsystem and the use of RAM practically do not change during research. This is largely due to the fact that the number of particles practically did not change during the study. That is, the amount of memory used, and the amount of data stored per disk is determined by the complexity of the 3D model of the equipment and the number of particles.
Changing the size of the «grid» within the ±1/3 limits that differ from those recommended by the developers of EDEM 2017 showed that this parameter has practically no effect on the duration of calculations.
Decisive, from the point of view of modeling time, is the number of particles. Our results showed that there is a linear relationship between these parameters, which makes it possible to predict the time required for calculations.
The phenomenon of superlinear acceleration, presented in [9], is not observed when performing calculations on a single-processor system. The speed of calculations linearly depends on the number of processor cores. The nature of the relationship between the duration of the simulation and the number of computational nodes (processor cores) corresponds to the results reported in [10]. This makes it possible to conclude that the developers of the EDEM 2017 system use methods of optimization of parallel calculations similar to the algorithms proposed by the authors of study [10]. The EDEM 2017 system parallelizes calculations by creating threads that use a common area of RAM with the main process. This mode of operation corresponds to the OpenMP method presented in [11]. The amount of RAM used by the system almost did not change during the simulation.
The method of load separation between the graphics and central processor, reported in [12], showed better results than those obtained for a system with a single processor. At the same time, researchers have the opportunity to increase performance through the use of processors with more cores, which avoids the cost of purchasing a GPU. The EDEM 2017 system implements the «classic» version of DEM calculations, which involves determining all possible interactions between particles according to the algorithm (Fig. 1). Methods and approaches to reduce the vo lume of calculations, namely, DEM-PBM [13], LS-DEM [14], and the method of domain decomposition [15], allow for modeling in significantly less time. However, the accuracy of the results may be lower compared to DEM, which requires more research.
Parallel algorithms for multiprocessor systems, which are presented in works [16,17], showed a nonlinear relationship between the number of processors for the duration of the simulation. At the same time, the effect of the introduction of additional processors into the system decreases with an increase in their number, which is a consequence of the impact of overhead costs on interprocessor data transmission. At the same time, the results of our study for a single-processor system showed a linear relationship between the duration of calculations and the number of processor cores. This makes it possible to conclude about the effectiveness of using multicore single-processor systems for DEM research.
The limitations of the study are due to the fact that the results are influenced by the type of equipment and the nature of the movement of particles. In belt feeders, most particles move in a continuous stream at an average speed not exceeding the speed of the belt (0.1 m/s). At the same time, the number of particle-to-particle interactions is much greater than the number of particle-equipment interactions. This determines the number of calculations in steps 5 and 6 of the algorithm (Fig. 1). This mode of movement is characteristic of equipment that is used to transport bulk materials. At the same time, in other types of equipment, in particular centrifugal mixers [19], particles move in rarefied flows, which significantly changes the number of particle-particle interactions. Accordingly, the number of calculations in steps 5 and 6 will also change, which will not make it possible to use the results obtained for equipment of this type without additional research.
The disadvantage of the study is a decrease in measurement accuracy due to the use of computational resources by other processes. Despite the fact that during the research all third-party applications were disabled, it is almost impossible to avoid the influence of system processes.
As part of the further development of the study, it is advisable to determine the influence of the nature of the movement of particles at the time of modeling by DEM, which will make it possible to predict the duration of simulation for equipment of various types.

Conclusions
1. A mathematical model of particle motion in a belt-type feeder has been built. The physical and mechanical proper-ties of particles and their geometric dimensions correspond to similar parameters of granules of polymeric materials, which are the basis for the manufacture of various plastic products. The properties of the materials of the feeder belt and its body correspond, respectively, to rubber and steel. The technological parameters of the feeder, in particular the speed of movement of the belt, are characteristic of light industry equipment. Thus, the proposed mathematical model for DEM research makes it possible to simulate the operation of equipment used in the manufacture of plastic products.
2. It was established that the software used for DEM research makes it possible to change the number of particles, the «grid» pitch, and the number of processor cores without affecting the accuracy of the calculation results.
3. The results of our experimental studies showed that for the duration of the simulation, determining parameters are the number of particles and the number of processor cores. It was established that there is a linear relationship between the duration of the simulation and the number of particles. The regression equation is derived, which makes it possible to predict the simulation time. It is also established that the software does not fully use all available computing resources, 57 % at best.