2021 Volume 62 Issue 4 Pages 551-556
Using the genetic algorithm, this paper developed an automatic method for calibrating the used four input parameters (Young’s modulus, Poisson’s ratio, coefficient of restitution, and friction coefficient) in the discrete element modeling (DEM) of vibration feeders. Two objective functions were used: (1) the difference in the results between the DEM simulation and the experiment and (2) the computation time of the DEM simulation. We conducted preliminary experiments to obtain an average experimental response of the vibration feeding and found that single-item-clogging usually arises near the feeder tray outlet. The proposed method managed to find the required set of parameters for minimizing the two objective functions by multi-objective optimization using the genetic algorithm. Overall, we confirmed that the obtained parameter values were physically reasonable in comparison with those of previous studies.
The discrete element method (DEM)1) has been commonly used in quantitatively evaluating the assembly of solid objects, such as powders. The DEM is one of the particle-based simulation methods, where the contact forces between any two rigid objects are approximated using the spring-dashpot-slider system, and the motion of all the contacting particles is explicitly calculated using time integration schemes. This is an essential difference between the DEM and particle-based methods when it comes to partial differential equations (PDEs), such as smoothed particle hydrodynamics2–4) and the moving-particle semi-implicit method.5) Usually, PDE-based simulations use constitutive models as input parameters. For example, the constitutive relation of the elastic body includes Young’s modulus and Poisson’s ratio, and that of the incompressible Newtonian fluid includes the density and viscosity. Such model parameters are determined by experiments. In contrast, using the Voigt contact model, the input parameters in conventional DEM methods are normal and tangential spring constants, normal and tangential damping coefficients, and friction coefficients. Thus, due to the variations in the particle characteristics, such as the particle shape and diameter, it is impractical to experimentally determine the parameters of all the particles in a DEM model.
With the practical use of the DEM in industrial applications, the input parameters are intuitively determined on the basis of the practitioners’ experiences. In previous academic studies, preliminary parametric studies were conducted before the main simulations. The input parameters were selected based on comparisons between the results of the parametric studies and the experiments. Katagiri et al. conducted a parametric study on the spring constant and evaluated the spring constant influence on the bulk shear6) and the one-dimensional compression7) responses. Matsushima et al. conducted parametric studies on the spring constant, dashpot damping coefficient, and friction coefficient, and they evaluated the suitability of the used parameters.8) Kikkawa et al. proposed a systematic procedure for determining the input parameters of the DEM using the results of parametric studies.9) The used parameter values in such past studies were based on the researchers’ subjective judgment. There are four parameters in the Voigt contact model in ROCKY DEM, which is the commercial DEM software that is used in this study, and they are Young’s modulus, Poisson’s ratio, the restitution coefficient, and the friction coefficient. Assuming that there are ten levels for the four parameters, the comprehensive survey would require 10,000 (104) simulations. Moreover, the parameter values can just be selected from the 10 levels of each parameter.
To improve the efficiency of the parametric study for the DEM, previous studies have applied design of experiment (DoE) methods to their parameter surveying processes.10,11) The response surface is formed by fitting the results of the DEM simulations with convex-shaped functions, such as the quadratic function,12) and the best parameter set can be determined by the mathematical minimization (or maximization) of the response surface. The main advantage of the DoE-based approach is that it can significantly reduce the number of DEM simulations in comparison with performing comprehensive parametric studies, while the main disadvantage is that the set of simulation results is supposed to be approximated by the convex response surface. The overall trend of the simulation results cannot have a convex shape, so the DoE-based approach may not find the best parameter set. Both the comprehensive parametric study and DoE-based approach are time-consuming and require manual operations of practitioners, and the preliminary parameter calibration process is a big burden for practitioners who use the DEM.
There is another problem in the parameter calibration process: the multiple input parameters contribute to the varied physical quantity during the parametric study process. For example, the restitution coefficient varies in the parametric study for evaluating the damping coefficient effect of the dashpot. However, not only the damping coefficient but also the spring constant influence the restitution coefficient.13,14) The friction coefficient is another example, and it is the only coefficient that varies in the parametric study. However, the tangential spring constant and tangential damping coefficient contribute to the friction force acting on the particle. The best parameter value cannot be determined independent of the other parameters just by the parametric study of single parameters. To avoid the burden of using the DEM by practitioners, a completely automatic method is required to determine the best parameter set for the DEM. Such technology will lead to an increase in the industrial use of the DEM. Hence, the mathematical optimization method can be applied to the calibration of the DEM input parameters. In contrast to the DoE-based approach, the optimization method can automatically find the best parameter set without presetting the levels of the parameters.
In Chapter 2, the automatic parameter calibration method and the used optimization algorithm are described in detail. The laboratory experiment and the required experimental responses in the automatic parameter calibration analysis are detailed in Chapter 3. In Chapter 4, the obtained parameters from the proposed method are discussed, and the conclusion is described in Chapter 5.
A flow chart of the automatic parameter calibration method is shown in Fig. 1, and it is explained in detail in the following sections.
Flow chart of the automatic parameter calibration method using the discrete element method and multi-objective genetic algorithm.
The governing equation of DEM is the equation of motion. The particles’ positions and velocities are explicitly updated with every increment using time integration scheme. The contact force is required to solve the equation of motion; hence the contact model and its parameters strongly influence the simulated responses of DEM. The contact forces along the normal and tangential directions are calculated using linear spring dashpot model and Coulomb limit model, respectively. The models and their correlation with the DEM parameters are detailed in the cited reference;15) the overview is described in section 2.2. The used granular material in this study was waste electrical and electronic equipment (WEEE). The WEEE has been increasing over the last decade,16) and the required technical development for WEEE recycling and its handling equipment is vital for environmental sustainability.
We used 30 WEEE samples: 10 smartphones (SPs), 10 feature phones (FPs), and 10 digital cameras (DCs). The three-dimensional (3D) shapes of the WEEE were detected using a structured light 3D scanning system (David SLS-3 stereo 3D scanner), and they were incorporated in the DEM simulation. The averages and standard deviations of volume, surface area, diameter of equivalent-volume sphere, and mass of the 30 WEEE samples are listed in Table 1. Note that the DEM simulations required a very long computation time when the original 3D shapes were directly used, and the data size of the original 3D data were reduced to 1/1000 based on the past study.17) We simulated the motions of the WEEE inside a tapered tray, as shown in Fig. 2. The widths of the rectangular region for the particle inlet and outlet were 260 mm and 100 mm, respectively, and the tray had the same geometry of an actual machine, as described in the next chapter.
The tapered tray geometry for the vibration feeder.
The 30 WEEE samples were randomly generated and deposited under gravity for 1 s at the particle inlet domain, which is shown in the left side of Fig. 2. The simulation of the gravitational deposition continued for 4 s. The tray was vibrated under a frequency of 60 Hz and an amplitude of 1.5 mm along with a 25° inclination from the horizontal plane for 60 s. The snapshots taken during the DEM simulations of the particle packing and vibration are shown in Fig. 3.
Snapshots during DEM simulation of vibration feeder: (a) initial state, (b) during particle insertion, (c) after gravitational deposition, (d) 1.4 s, (e) 2.8 s, (f) 6.2 s, (g) 10.2 s, (h) 15 s, and (i) 18 s from vibration start.
The DEM software is ROCKY DEM (version 4.1.2), which includes four contact parameters to be determined: Young’s modulus (E), Poisson’s ratio (ν), the restitution coefficient (eb), and the friction coefficient (μ).
When two spheres with same diameter (D) are in contact, the linear spring dashpot model calculates the normal force (Fn) as per following equation:15)
| \begin{equation} F_{\text{n}} = \frac{ED}{2}S_{\text{n}} + \eta \sqrt{mED} \dot{S}_{\text{n}}, \end{equation} | (1) |
| \begin{equation} F_{\text{t}} = -\mu F_{\text{n}}\widehat{t_{\text{c}}}, \end{equation} | (2) |
The E and ν are known as material constants. The E and ν of the WEEE cannot be determined based on the constituent materials because their component parts are complexly arranged inside the WEEE. The eb influences the viscous damping in the damped vibration, and its value is difficult to determine because of the complexity of the component materials in the WEEE. In this study, the four DEM parameters were automatically determined using the proposed method.
2.3 Genetic algorithmThe used parameter optimization algorithm is the genetic algorithm (GA),18) which is one of the commonly used evolutionary computation algorithms. As described below, two evaluation (objective) functions were used, so this study applied the multi-objective GA (MOGA), which was implemented in the ANSYS DesignXplorer software.
In the MOGA, the parameters to be determined are regarded as gene information. The problem in this study was to determine the values of the four parameters, i.e., the genes. The flow of computation of MOGA is as follows. First, the initial population is generated: 50 of the individuals who possess different values of the four parameters. The individuals have different parameters. The DEM simulations of the vibration feeder use the parameter set of each individual. Then, the values of the objective function are calculated using the result of the DEM simulation. By comparing the values of the objective function, the individuals carried over to the next generation are chosen from the initial population by elitist selection.19) Then, the new individuals in the next generation, i.e., the children of the initial population, are generated by two mechanisms: (1) crossover of the individuals who survived in the elitist selection and (2) mutation. A series of the described processes above is applied to the individuals of the children generation and then the individuals of the grandchildren generation are generated. The optimum parameter set can be determined by a so-called Darwinism: the concept of the survival of the fittest.
The numbers in the initial population and the number of individuals carried over to the next generation were 50 and 15, respectively. We used the one-point crossover method19) and typical values of the crossover and mutation probabilities of 0.98 and 0.02, respectively.
2.4 Objective functionWe used two objective functions: the number of WEEE items remaining in the tray after the 60-s vibration simulation (Nt) and the computation time required 65-s DEM simulation including 5-s particles deposition and 60-s vibration (T). The Nt range was between 0 and 30. When Nt was the same as that obtained from the experiments, the DEM could reproduce the experimental response of the vibration feeder. The T value was greater than 0. It is preferable that the computation time is shorter if the DEM can simulate the experimental response. Then, the problem in this study is to determine the parameter set so as to equate the obtained Nt from the experiments and minimize T.
In this study, we used the electromagnetic vibration feeder: ERIEZ 36C. The tapered tray shown in Fig. 2 was set on the vibration feeder, and the vibration amplitude and frequency were 1.5 mm and 60 Hz, respectively. The 30 WEEE items, which consisted of 10 SPs, 10 FPs, and 10 DCs, were randomly inserted to the tray.
The vibration continued for 60 s. After the 60-s vibration, we recorded the value of Nt. When all the items were discharged before 60 s, the discharging time was recorded. Ten experiments were conducted because the initial configuration and orientation of the WEEE items would influence the Nt and the discharging time.
Note that the used 30 WEEE items used in the experiment were different from the set of the used 30 items used in the DEM simulation. We have conducted preliminary experiments using different sample sets of the 30 WEEE items (consists of 10 SPs, 10 FPs, and 10 DCs); as a result, we confirmed that the sample set hardly affected Nt value when the set did not contain a large enough item compared with the tray outlet.
3.2 ResultsFigure 4 shows the taken snapshots during the experiment. The WEEE items were concentrated near the narrow space that is next to the tray outlet, so the clogging frequently increased during the experiments. Clogging due to the particles’ bridge was often observed in the powder discharging from the hoppers and silos. Bridging of some of the WEEE items was observed in the experiment. In contrast, single-item-clogging, where a WEEE item is stuck around the outlet due to its orientation angle, was frequently observed in the experiments. The single-item-clogging of some of the samples was sometimes unstuck during the vibration. Such temporal single-item-clogging can be seen in Fig. 4(b), and it was unstuck after a 56-s vibration, as shown in Fig. 4(d). Note that permanent single-item-clogging was also often observed in the experiments.
Snapshots during the vibration feeder experiment at a vibration amplitude of 1.5 mm and a frequency of 60 Hz: (a) initial state, (b) 5 s, (c) 20 s, and (d) 56 s from the vibration start.
Due to the two kinds of clogging, bridging and single-item-clogging, all the WEEE items were discharged just in two out of the ten experiments. The median, average, and standard deviations of the values of Nt were 3.5, 5.7, and 6.2, respectively. Since the Nt is the integer value, we decided to use Nt = 4 in the automatic parameter determination analysis.
The typical responses of Nt and T during the parameter optimization process are shown in Fig. 5. The horizontal axis in the figure represents the evolution of the generation in the GA. The simulated values of Nt and T decreased with the evolution of the generation. The Nt value varied because various factors other than the DEM parameters could influence the clogging. However, the median value of the Nt converged to that of the experiment after the 9th generation. Moreover, the T value decreased and became constant, which implies that the two objective functions were successfully optimized using the proposed method.
The generation evolution influence of on (a) Nt and the (b) computation time.
The obtained parameter values by the proposed method are listed in Table 2. The values of E were relatively low in the parameter optimization range, and the E value considerably influenced the computation time. The particles’ displacement and rotation were explicitly calculated every time increment. When the same overlap distance was between the two objects in the DEM, the larger contact force was calculated using the higher value of E. The repulsive velocity became high as the particle jumped out from the simulation domain. To prevent such erroneous calculations, the small-time increment was often used. Also, the number of time steps was considerably increased when the time increment became small. As a result, the computation time significantly increased.
The theoretical range of ν was −1 < ν < 0.5.20) However, most of the materials had a range of 0 < ν < 0.5. The obtained ν value was 0.297, and it is similar to the ν = 0.32 of the aluminum alloy, which is a typical material used for the housing body of WEEE.21)
The μ value depends on various factors. Among them, the constituent material and surface state (wet or dry) strongly influence the friction coefficient. In the dried condition, the static (μs) and dynamic μ (μd) were evaluated in previous experimental studies. The μs values of steel–steel, aluminum–steel, and copper–steel are 0.74,22) 0.61,23) and 0.53,23) respectively. The μd values of steel–steel, aluminum–steel, and copper–steel are 0.57, 0.47, and 0.36, respectively.24) The obtained μ value was 0.501, which is a reasonable value compared with those of the past experiments. Note that the used static and dynamic μ values were the same in this study.
The damping coefficient, the elastic stiffness such as E, and the spring constant influence the value of eb. The derived vibration force from the vibrating tray moved the particles toward the outlet. An extremely low eb caused overdamping. As a result, the particles’ movement was strongly restrained. However, the particles’ behavior was unnatural with a very high eb. The experimental eb values of the aluminum alloy and polycarbonate spheres are 0.875 and 0.884, respectively.21) When the particle shape was non-spherical, its rebound direction was unpredictable. In such case, the eb value was influenced by the inclined angle from the horizontal plane.25) We could not find the experimental evaluation of eb for the WEEE in the literature. The WEEE includes void spaces inside the housing body. The derived energy from the collision was dissipated by the housing deformation, so we were convinced that the eb of the WEEE should be low.
This study developed an automatic parameter calibration method using the genetic algorithm (GA). We used the two objective functions that represent the reproducibility of the experimental response and computational efficiency. The conventional comprehensive parametric study problem was replaced by the problem of minimizing the two objective functions. The proposed method was applied to determine the input parameters for the DEM simulation of the vibration feeder. Considering the physical meaning of the parameters, we confirmed that the obtained parameter values were reasonable.
As described in Chapter 1, the 104 simulations required a comprehensive survey for 10 levels of the 4 parameters. In contrast, the parameter set, which could reproduce the experimental response, was obtained by 185 simulations during the automatic parameter calibration method, so high efficiency is an advantage of the proposed method. This paper described the application of the proposed method for the vibration feeding simulation. The WEEE recycling includes various processes such as dismantling, crushing, and physical separation. As the typical experimental response has been evaluated, the proposed method was found to have no limitation of the applied problem; application flexibility is another advantage of the proposed method.
The values of the used GA parameters were small compared with the typical values in the recent applications that used the GA. However, the influence of the GA parameters on the results of the automatic parameter calibration should be evaluated.
The authors believe that there are various combinations of the DEM parameters for reproducing the experimental response, which means that parameter calibration is a multimodal optimization problem. We aim to evaluate the reproducibility of the proposed method in future studies.
This paper is based on results obtained from a project commissioned by the New Energy and Industrial Technology Development Organization (NEDO).
