Calibration and Verification of Dynamic Particle Flow Parameters by the Back-Propagation Neural Network Based on the Genetic Algorithm: Recycled Polyurethane Powder

The discrete element method (DEM) is commonly used to study various powders in motion during transportation, screening, mixing, etc.; this requires several microscopic parameters to characterize the complex mechanical behavior of the particles. Herein, a new discrete element parameter calibration method is proposed to calibrate the ultrafine agglomerated powder (recycled polyurethane powder). Optimal Latin hypercube sampling and virtual simulation experiments were conducted using the commercial DEM software; the microscopic variables included the static friction coefficient between the particles, collision recovery coefficient, Johnson–Kendall–Roberts surface energy, static friction coefficient between the particles and wall, and collision recovery coefficient. A predictive model based on genetic-algorithm-optimized feedforward neural network (back propagation) was developed to calibrate the microscopic DEM simulation parameters. The cycle search algorithm and mean-shift cluster analysis were used to confirm the input parameters’ range by comparing the mean value of the dynamic angle of repose measured via the batch accumulation test. These parameters were verified by the baffle lifting method and the rotating drum method. This calibration method, once successfully developed, will be suitable for use in a variety of fine viscous powder dynamic flow conditions.


Introduction
Currently, the main methods that are used for recovering and recycling thermosetting plastic waste include the mechanical and physical method, the chemical method, and the energy recovery method. The mechanical and physical method [1,2] can achieve complete recovery with high recovery efficiency and low pollution. The recycling method is based on the principle that materials accumulate under various mechanical forces, causing the accumulation of mechanical energy, destroying the crosslinked molecular structural network, generating a recycled powder with a low crosslinking degree, and realizing the recovery and recycling of the thermosetting plastic [3].
In industrial production, as shown in Figure 1, in the process of mixing and regenerating the material after mechanical force pulverization, the active regenerated material is usually added, mixed, and regenerated; then, a reagent such as a crosslinking agent is added and solidified to form a thermoplastic finished material. However, during the production process, the powder is often in calibration simulation, it is necessary to use the calibration test under complex dynamic conditions to characterize the actual mechanical behavior of the powder. Recently, the discrete element method (DEM) has become a popular tool for describing the movement of particles, including powders, seeds, and soil, especially in the mining, chemical, pharmaceutical, and environmental industries, with the improvement of the computer capabilities [4,5]. DEM has been demonstrated to be a better tool than other simulation methods at this stage because it can analyze and optimize the process parameters associated with the preparation of the particle material and the operation of equipment used in the industry. In DEM simulations, several physical properties of the materials involved in collision and the contacts that occur between them are readily available, including the parameters that describe the intrinsic properties of a single contact material such as its density, size, Poisson's ratio, and Young's modulus. However, the contact characteristics will depend upon the specific model that is used while describing the interaction between particles and between a particle and the device geometry. The commonly used contact models are the Hertz-Mindlin, Johnson-Kendall-Roberts (JKR), and Bonding models; however, these models require the calibration of a large number of contacts. Further, the parameters can be used to accurately characterize the actual working conditions [6,7]. In the DEM parameter calibration test, the application of direct measurement calibration or indirect virtual calibration is the primary method for obtaining the aforementioned parameters. The characteristic parameters that are easy to measure can be directly measured experimentally, and the measurement results can be considered to be the input value for the DEM parameters. However, because of the anisotropy of some particulate matter, the direct measurement results are observed to vary considerably. The difference between the established particle model and the actual model is large, and the parameter measurement value cannot be directly applied to the DEM numerical simulation. Therefore, numerous researchers have proposed the usage of virtual simulation experiments to calibrate the contact parameters required in discrete element simulations and have conducted extensive research [8,9]; however, these calibrations mostly use the "testing method", which exhibits high computing costs. To solve this problem, researchers have attempted to use predictive models or apply optimization algorithms. Consider Wu [10] as an example, the triaxial shear tester and stress-strain curve were used to determine the macroscopic mechanical properties of the soil. Then, the genetic neural network was used to derive the DEM parameters. Cheng [11] proposed a sequential quasi-Monte Carlo filter as a new calibration method for the DEM models of granular materials. Rackl and Hanley [12] used Latin hypercube sampling and the kriging method to describe and demonstrate a method based on the static angle method. Zhou et al. [13] calibrated the irregularly expanded graphite particles, combining optimal Latin hypercube sampling and virtual simulation experiments to obtain the DEM parameters that are to be calibrated using an adaptive simulated annealing optimization algorithm; further, the membership function was selected. The output Recently, the discrete element method (DEM) has become a popular tool for describing the movement of particles, including powders, seeds, and soil, especially in the mining, chemical, pharmaceutical, and environmental industries, with the improvement of the computer capabilities [4,5]. DEM has been demonstrated to be a better tool than other simulation methods at this stage because it can analyze and optimize the process parameters associated with the preparation of the particle material and the operation of equipment used in the industry. In DEM simulations, several physical properties of the materials involved in collision and the contacts that occur between them are readily available, including the parameters that describe the intrinsic properties of a single contact material such as its density, size, Poisson's ratio, and Young's modulus. However, the contact characteristics will depend upon the specific model that is used while describing the interaction between particles and between a particle and the device geometry. The commonly used contact models are the Hertz-Mindlin, Johnson-Kendall-Roberts (JKR), and Bonding models; however, these models require the calibration of a large number of contacts. Further, the parameters can be used to accurately characterize the actual working conditions [6,7]. In the DEM parameter calibration test, the application of direct measurement calibration or indirect virtual calibration is the primary method for obtaining the aforementioned parameters. The characteristic parameters that are easy to measure can be directly measured experimentally, and the measurement results can be considered to be the input value for the DEM parameters. However, because of the anisotropy of some particulate matter, the direct measurement results are observed to vary considerably. The difference between the established particle model and the actual model is large, and the parameter measurement value cannot be directly applied to the DEM numerical simulation. Therefore, numerous researchers have proposed the usage of virtual simulation experiments to calibrate the contact parameters required in discrete element simulations and have conducted extensive research [8,9]; however, these calibrations mostly use the "testing method", which exhibits high computing costs. To solve this problem, researchers have attempted to use predictive models or apply optimization algorithms. Consider Wu [10] as an example, the triaxial shear tester and stress-strain curve were used to determine the macroscopic mechanical properties of the soil. Then, the genetic neural network was used to derive the DEM parameters. Cheng [11] proposed a sequential quasi-Monte Carlo filter as a new calibration method for the DEM models of granular materials. Rackl and Hanley [12] used Latin hypercube sampling and the kriging method to describe and demonstrate a method based on the static angle method. Zhou et al. [13] calibrated the irregularly expanded graphite particles, combining optimal Latin hypercube sampling and virtual simulation experiments to obtain the DEM parameters that are to be calibrated using an adaptive simulated annealing optimization algorithm; further, the membership function was selected. The output values were used to obtain several sets of optimal contact parameter sets. Do et al. [14] developed a parameter calibration framework for discrete element models based on a multi-objective genetic algorithm (GA). They used quartz sand particles to optimize the model accuracy and simulation time and demonstrated that the contact parameters of the materials can be accurately calibrated within the optimal simulation time. Simultaneously, coarse granulation is the most common modeling approach for the calibration of the discrete element parameters of special particles (including wet particles and ultrafine polymer powders), as demonstrated by Alizadeh et al. [15] for polyethylene glycol 400 (PEG 400). In a study investigating the coating treatment of the EP particles, the concept of cohesive number was used in the DEM numerical simulation to measure the shear modulus of the material or to change the particle size to ensure that the surface energy could be rapidly scaled and that the calibration time could be substantially reduced. Nasato et al. [16] used DEM to numerically simulate the wet granulation of continuous granulators and investigated the possibility of seed-particle formation during a continuous process to reduce the number of trials and processes. Further, a control method was designed to define dimensionless number (cohesion number) on the basis of the effect of particle cohesion and gravitational potential energy to investigate the effects of factors, such as the drum rotation speed, particle surface energy, and particle size ratio, on the performance of the seed granulator. Among them, reducing the shear modulus of the particles and improving the calculation ability are also necessary content for coarse grain modeling. For example, J. Haervig et al. [17] gave an analytical derivation criterion for reducing computation time by reducing particle stiffness. The validity of the criterion was verified by comparing the experimental data with the particle simulation of reducing the particle stiffness. Yile Gu et al. [18] used an improved cohesion model for fluidized CFD-DEM simulations. For viscous particles, the predicted flow pattern depends on the value of the particle spring stiffness used in the simulation. An improved cohesive model was proposed, which was verified by the simulation of two-particle collision, which provides satisfactory results for the fluidization of particles.
Therefore, based on the aforementioned summary and given that the recycled polyurethane (PU) powder investigated in this study also exhibits strong adhesion, we selected the JKR model exhibiting the agglomeration bonding effect as the characterization method to select the most sensitive discrete meta-contact parameters. Simultaneously, we used the sample space obtained after Latin hypercube sampling (LHS) sampling to simulate a series of parameter values on the basis of the determined range of the contact parameters via the commercial DEM software. Finally, the back-propagation (BP) neural network is trained to obtain an accurate prediction model. The calibration process is divided into two phases. The first phase is based on the development of a BP neural network prediction model optimized by the GA (Figure 2a), whereas the second phase is based on a neural network model developed for specific particulate materials. The application of the calibration experiment resulted in the selection of the DEM parameters (Figure 2b).
Training a BP neural network requires a series of DEM simulations in the selected range of simulation parameters. Herein, the BP neural network was trained in conjunction with the improved FBS-104 powder meter test. In engineering applications [19], the correlation coefficient (R 2 ) is usually greater than 90% before the neural network can be used as the predictive model; thus, the correlation coefficient can be considered to be an evaluation criterion. Figure 2b denotes the application of the BP neural network after training. The same powder measuring device is used to test the PU sample to avoid introducing an error into the prediction model. Based on the trained BP neural network database, the output will be the angle of repose, which is input into the loop search algorithm, and the mean-shift cluster analysis method is used to solve the feasible domain of the DEM parameters. Then, the candidate DEM parameter set are determined from the database and verified using the baffle lift and rotary drum methods. This calibration method, once successfully developed, will be suitable for use in a variety of fine viscous powder dynamic flow conditions. Training a BP neural network requires a series of DEM simulations in the selected range of simulation parameters. Herein, the BP neural network was trained in conjunction with the improved FBS-104 powder meter test. In engineering applications [19], the correlation coefficient (R 2 ) is usually greater than 90% before the neural network can be used as the predictive model; thus, the correlation coefficient can be considered to be an evaluation criterion. Figure 2b denotes the application of the BP neural network after training. The same powder measuring device is used to test the PU sample to avoid introducing an error into the prediction model. Based on the trained BP neural network database, the output will be the angle of repose, which is input into the loop search algorithm, and the mean-shift cluster analysis method is used to solve the feasible domain of the DEM parameters. Then, the candidate DEM parameter set are determined from the database and verified using the baffle lift and rotary drum methods. This calibration method, once successfully developed, will be suitable for use in a variety of fine viscous powder dynamic flow conditions.

Funnel Flow Meter Test Device
The test material was obtained by pulverizing the thermosetting PU closed rigid foam plastic in a homemade thermosetting plastic pulverization and reclaiming test machine; powders with particle sizes of 80 (0.18 mm) and 120 mesh (0.125 mm) and a small amount of 200-mesh (0.0750 mm) powder were obtained. During the calibration process, the PU powder of the calibration object is prone to agglomeration; therefore, agglomeration occurs in the rotating drum, leading to the formation of a dynamic accumulation angle and the occurrence of the "avalanche" phenomenon, which is difficult to overcome [20]. Therefore, the use of a funnel flow meter (modified FBS-104 powder meter) can improve the cumulative mechanical properties caused by the angle. The device comprised a 120-mm funnel, 94-mm diameter circular loading platform, brush, and Plexiglas bracket ( Figure 3). However, the powder cannot be normally formed to reduce the "pour" phenomenon of the powder particles during the free-fall movement of the PU powder. The funnel outlet was designed to have a small 7° angle of inclination to reduce the material drop rate. During the whole test process, the material slides down along the funnel wall with the brush at 60 rpm. Finally, the angle of repose of the powder particles can be measured when the powder to be discharged accumulates on the receiving table to form a stable material pile.

Funnel Flow Meter Test Device
The test material was obtained by pulverizing the thermosetting PU closed rigid foam plastic in a homemade thermosetting plastic pulverization and reclaiming test machine; powders with particle sizes of 80 (0.18 mm) and 120 mesh (0.125 mm) and a small amount of 200-mesh (0.0750 mm) powder were obtained. During the calibration process, the PU powder of the calibration object is prone to agglomeration; therefore, agglomeration occurs in the rotating drum, leading to the formation of a dynamic accumulation angle and the occurrence of the "avalanche" phenomenon, which is difficult to overcome [20]. Therefore, the use of a funnel flow meter (modified FBS-104 powder meter) can improve the cumulative mechanical properties caused by the angle. The device comprised a 120-mm funnel, 94-mm diameter circular loading platform, brush, and Plexiglas bracket ( Figure 3). However, the powder cannot be normally formed to reduce the "pour" phenomenon of the powder particles during the free-fall movement of the PU powder. The funnel outlet was designed to have a small 7 • angle of inclination to reduce the material drop rate. During the whole test process, the material slides down along the funnel wall with the brush at 60 rpm. Finally, the angle of repose of the powder particles can be measured when the powder to be discharged accumulates on the receiving table to form a stable material pile.

Angle-of-Repose Measurement
The angle of repose is used to describe the application of material particles in specific scenarios and the corresponding mechanical behavior, including the liquid bridge force, friction force, van der Waals force, and other forces that often exist between particles. Therefore, the material particles are studied by considering the angle of repose as a macroscopic response value to ensure that the appropriate microscopic DEM parameters of the material can be indirectly identified based on these response values. The angle of repose is most commonly measured with a material in its unconstrained state, as depicted in Figure 4a; the material particles are stacked and formed on the horizontal surface without collapse [21], and the slope of the formed material is observed, which can be converted into the material accumulation angle. Materials 2019, 12, x FOR PEER REVIEW 5 of 24

Angle-of-Repose Measurement
The angle of repose is used to describe the application of material particles in specific scenarios and the corresponding mechanical behavior, including the liquid bridge force, friction force, van der Waals force, and other forces that often exist between particles. Therefore, the material particles are studied by considering the angle of repose as a macroscopic response value to ensure that the appropriate microscopic DEM parameters of the material can be indirectly identified based on these response values. The angle of repose is most commonly measured with a material in its unconstrained state, as depicted in Figure 4a; the material particles are stacked and formed on the horizontal surface without collapse [21], and the slope of the formed material is observed, which can be converted into the material accumulation angle. Figure 4b-d denote that the PU powder pile exhibits good symmetry. Alizadeh et al. [15] proposed that the high material is generally 10%-90% to avoid large fluctuations in the boundary contour of the material particles during stacking. For fitting the final boundary contour, the contour is read using the computer image processing technology and then processed by gradation and binarization; the boundary is subsequently extracted, and the boundary point is linearly fitted using the least-squares method to obtain the left and right boundaries of the powder particle stack. The equation of the contour is linearly fit, and the average value measured using the left and right angles is the angle of repose. The specific process is presented in Figure 4a-d.

Angle-of-Repose Measurement
The angle of repose is used to describe the application of material particles in specific scenarios and the corresponding mechanical behavior, including the liquid bridge force, friction force, van der Waals force, and other forces that often exist between particles. Therefore, the material particles are studied by considering the angle of repose as a macroscopic response value to ensure that the appropriate microscopic DEM parameters of the material can be indirectly identified based on these response values. The angle of repose is most commonly measured with a material in its unconstrained state, as depicted in Figure 4a; the material particles are stacked and formed on the horizontal surface without collapse [21], and the slope of the formed material is observed, which can be converted into the material accumulation angle. Figure 4b-d denote that the PU powder pile exhibits good symmetry. Alizadeh et al. [15] proposed that the high material is generally 10%-90% to avoid large fluctuations in the boundary contour of the material particles during stacking. For fitting the final boundary contour, the contour is read using the computer image processing technology and then processed by gradation and binarization; the boundary is subsequently extracted, and the boundary point is linearly fitted using the least-squares method to obtain the left and right boundaries of the powder particle stack. The equation of the contour is linearly fit, and the average value measured using the left and right angles is the angle of repose. The specific process is presented in Figure 4a-d.

Coarse-Grain Modeling of the Particle Discrete Element Based on the JKR Contact Model
In this calibration test, we mainly studied a thermosetting PU powder with a particle size of 120 mesh (0.125 mm). A scanning electron microscopy image of the powder is depicted in Figure 5, denoting that the powder particles are mostly not spherical. Therefore, the mechanical interlocking of the particles may substantially contribute to the volumetric shear strength of the particulate material and may affect its fluidity. Härtl and Ooi [22] denoted that geometric interlocking of the   [15] proposed that the high material is generally 10%-90% to avoid large fluctuations in the boundary contour of the material particles during stacking. For fitting the final boundary contour, the contour is read using the computer image processing technology and then processed by gradation and binarization; the boundary is subsequently extracted, and the boundary point is linearly fitted using the least-squares method to obtain the left and right boundaries of the powder particle stack. The equation of the contour is linearly fit, and the average value measured using the left and right angles is the angle of repose. The specific process is presented in Figure 4a-d.

Coarse-Grain Modeling of the Particle Discrete Element Based on the JKR Contact Model
In this calibration test, we mainly studied a thermosetting PU powder with a particle size of 120 mesh (0.125 mm). A scanning electron microscopy image of the powder is depicted in Figure 5, denoting that the powder particles are mostly not spherical. Therefore, the mechanical interlocking of the particles may substantially contribute to the volumetric shear strength of the particulate material and may affect its fluidity. Härtl and Ooi [22] denoted that geometric interlocking of the nonspherical particles strongly affects the in vivo friction of the particulate materials and should be introduced in the DEM simulations to capture the volumetric shear strength of the material. Simultaneously, Bharadwaj et al. [23] proposed that duplicating the exact shape of the particles is unnecessary for calibrating the powder discrete element parameters; they observed that the use of double-sphere paired particles with an aspect ratio of greater than 1 can provide a satisfactory quantitative prediction and that the aspect ratio could be increased to 2 without resulting in any significant difference. Given that the PU powder used in this study exhibits a relatively complicated mechanical behavior and avoids large errors associated with the establishment of the DEM model, it is further verified in Section 2.3.3. Thus, we can conclude that the ball pairing with an aspect ratio of 1.3 was appropriate. This model is a suitable replacement for the PU powder ( Figure 6). In this calibration test, we mainly studied a thermosetting PU powder with a particle size of 120 mesh (0.125 mm). A scanning electron microscopy image of the powder is depicted in Figure 5, denoting that the powder particles are mostly not spherical. Therefore, the mechanical interlocking of the particles may substantially contribute to the volumetric shear strength of the particulate material and may affect its fluidity. Härtl and Ooi [22] denoted that geometric interlocking of the nonspherical particles strongly affects the in vivo friction of the particulate materials and should be introduced in the DEM simulations to capture the volumetric shear strength of the material. Simultaneously, Bharadwaj et al. [23] proposed that duplicating the exact shape of the particles is unnecessary for calibrating the powder discrete element parameters; they observed that the use of double-sphere paired particles with an aspect ratio of greater than 1 can provide a satisfactory quantitative prediction and that the aspect ratio could be increased to 2 without resulting in any significant difference. Given that the PU powder used in this study exhibits a relatively complicated mechanical behavior and avoids large errors associated with the establishment of the DEM model, it is further verified in Section 2.3.3. Thus, we can conclude that the ball pairing with an aspect ratio of 1.3 was appropriate. This model is a suitable replacement for the PU powder ( Figure 6).  The contact model is the core of the particle DEM. The JKR contact model [24][25][26] is an extension of the Hertz contact theory. The model is based on elastic solid contact and focuses on the "overlap" of contact between particles. The particle-matching ball mainly has a tangential elastic force, and the normal dissipative force and the tangential dissipative force interact. Two spherical particles are assumed to agglomerate only on the contact surface (Figure 7). When the external load is P0 and the elastic displacement is δ, the particle contact surface radius α0 can be obtained using the The contact model is the core of the particle DEM. The JKR contact model [24][25][26] is an extension of the Hertz contact theory. The model is based on elastic solid contact and focuses on the "overlap" of contact between particles. The particle-matching ball mainly has a tangential elastic force, and the normal dissipative force and the tangential dissipative force interact. Two spherical particles are assumed to agglomerate only on the contact surface ( Figure 7). When the external load is P 0 and the elastic displacement is δ, the particle contact surface radius α 0 can be obtained using the Hertz contact theory if the particles are not agglomerated. However, if the particles agglomerate, the contact surface radius α 1 is greater than α 0 although the external load is still P 0 . The contact model is the core of the particle DEM. The JKR contact model [24][25][26] is an extension of the Hertz contact theory. The model is based on elastic solid contact and focuses on the "overlap" of contact between particles. The particle-matching ball mainly has a tangential elastic force, and the normal dissipative force and the tangential dissipative force interact. Two spherical particles are assumed to agglomerate only on the contact surface ( Figure 7). When the external load is P0 and the elastic displacement is δ, the particle contact surface radius α0 can be obtained using the Hertz contact theory if the particles are not agglomerated. However, if the particles agglomerate, the contact surface radius α1 is greater than α0 although the external load is still P0. In industrial DEM modeling, a method of amplifying the particle size and reducing the shear modulus is adopted because the number of particles may be in the order of several million, leading to a long calculation time; however, a scientific and reasonable scale basis is required. Thus, the accurate surface energy value, particle size, and mechanical properties should be coordinated in the JKR contact model. To achieve this, the concept of cohesion number [27] is proposed. The cohesion number is expressed as the ratio of cohesive work between particles to the gravitational potential energy of the particles.

Work of cohesion
Cohesion number Gravitational potential e = nergy (1) This can be mathematically expressed as follows: In industrial DEM modeling, a method of amplifying the particle size and reducing the shear modulus is adopted because the number of particles may be in the order of several million, leading to a long calculation time; however, a scientific and reasonable scale basis is required. Thus, the accurate surface energy value, particle size, and mechanical properties should be coordinated in the JKR contact model. To achieve this, the concept of cohesion number [27] is proposed. The cohesion number is expressed as the ratio of cohesive work between particles to the gravitational potential energy of the particles.

Cohesion number =
Work of cohesion Gravitational potential energy (1) This can be mathematically expressed as follows: where r denotes the density of the particles, g denotes the acceleration of gravity, γ denotes the average surface energy of the particles, E* denotes the equivalent Young's modulus of elasticity, and R* denotes the equivalent particle size from among various particles. The parameter R* is calculated according to Equation (3).
where R 1 and R 2 are the radii of the contacting spheres, E 1 and E 2 are their Young's moduli, and m 1 and m 2 are their Poisson's ratios. During the simulation, to avoid large overlaps in the process of contact between particles, which would increase the calculation cost, the cohesion number (Coh) is usually used to scale the particle discrete element model; further, the particle size is increased according to a certain zoom level. In addition, the equivalent surface energy is calculated by increasing the particle size and reducing the shear modulus. In addition, Alizadeh et al. [15] used this method. For example, the average surface energy of the PU particles is 0.0214 J/m 2 [28]; when this value is substituted into Equation (2), the cohesion number is observed to be 0.000622. When the average surface energy value is maintained constant, the particle size is magnified 10 times. However, when the modulus is reduced by two orders of magnitude, the equivalent surface energy used in the simulation is 4.0 J/m 2 . Simultaneously, to accurately calibrate the calculated parameter values, we refer to the surface energy values of similar materials in the GEMM database provided by the EDEM software and expand the calibration range of the surface energy values with the equivalent surface energy value (4.0 J/m 2) as the intermediate value.
Expanding the surface energy value to both sides, we determine the value between 3 and 8 J/m 2 .

DEM Simulation Tests
The discrete element virtual calibration test is completely consistent with the test device described in Section 2.1. Particle factory generation was compiled using the commercial DEM simulation software with an API open-source interface; the filling rate was 0.56, and the generation time was 0.01 s. The funnel was filled with the particle group, and the large moment between the particles was eliminated. After the particles were allowed to stand for 3 s, the brush was rotated at 60 rpm. At this time, the particle group fell freely, with a gravitational acceleration of 9.81 m/s 2 (Figure 8), and finally formed a stable pile on the receiving table. The powder particle pile simulation time remained 20 s for the whole process, and the final 2 s was used for the excessive porosity when the particles piled up (this pile-up leads to the excessive formation of a simulated packing). The simulation test method based on that described in Section 2.3.1 was used to propose a coarse-grained modeling criterion based on the JKR model to characterize the actual accumulation phenomenon of the PU powder.

Sensitivity Analysis of the Particle Intrinsic Parameters
To simplify the calibration process, we analyze the effect of the particle Eigen parameters on the angle of repose. Therefore, in combination with the range of values of the DEM simulation parameters shown in Table 1, the Eigen parameters (density, Young's modulus, Poisson's ratio, and particle aspect ratio) of the particles were investigated by the simulation test described in Section 2.3.2. Thus, we can conclude that the angle of repose was affected.

Sensitivity Analysis of the Particle Intrinsic Parameters
To simplify the calibration process, we analyze the effect of the particle Eigen parameters on the angle of repose. Therefore, in combination with the range of values of the DEM simulation parameters shown in Table 1, the Eigen parameters (density, Young's modulus, Poisson's ratio, and particle aspect ratio) of the particles were investigated by the simulation test described in Section 2.3.2. Thus, we can conclude that the angle of repose was affected.  Figure 9a denotes that the Young's modulus is between 3.5 × 10 6 and 3.5 × 10 7 Pa and that the angle of repose gradually increases. Notably, almost no change is observed in the stacking angle when the Young's modulus becomes greater than 3.5 × 10 7 Pa. As described in Section 2.3.1, the computational time of discrete element simulation can be minimized by reducing the magnitude of Young's modulus. In this case, the Young's modulus is reduced by two orders of magnitude (set to 3.5 × 10 7 Pa). The mid-rest angle has almost no effect. In addition, the particle density has little effect on the angle of repose (Figure 9b), consistent with the results obtained from other studies [29]. Similarly, the particle Poisson's ratio increases from 0.3 to 0.7. Only a small reduction is observed in the static angle of repose of the particles (approximately 0.5). In the results of Bharadwaj et al. [23], the best aspect ratio problem for coarse particle granulation is proposed in Section 2.3.1, as depicted in Figure 9d. When the width ratio becomes 1.3, the angle of repose of the particles tends to be stable, consistent with the prediction of the static angle.
A set of thermodynamic PU powders with a particle size of 120 mesh (0.125 mm) was obtained through the aforementioned batch simulation test; the intrinsic parameters were Gp = 3.5 × 10 7 Pa, ρ P = 60 kg·m −3 , Vp = 0.4, particle aspect ratio = 1.3. on the angle of repose (Figure 9b), consistent with the results obtained from other studies [29]. Similarly, the particle Poisson's ratio increases from 0.3 to 0.7. Only a small reduction is observed in the static angle of repose of the particles (approximately 0.5). In the results of Bharadwaj et al. [23], the best aspect ratio problem for coarse particle granulation is proposed in Section 2.3.1, as depicted in Figure 9d. When the width ratio becomes 1.3, the angle of repose of the particles tends to be stable, consistent with the prediction of the static angle. A set of thermodynamic PU powders with a particle size of 120 mesh (0.125 mm) was obtained through the aforementioned batch simulation test; the intrinsic parameters were Gp =

BP Neural Network Database
Multidimensional sampling of the calibrated discrete element contact parameters was performed using a uniform distribution-based LHS method [30] to form a virtual test sample space. The distribution function domain of each stochastic input variable is equally divided into

BP Neural Network Database
Multidimensional sampling of the calibrated discrete element contact parameters was performed using a uniform distribution-based LHS method [30] to form a virtual test sample space. The distribution function domain of each stochastic input variable is equally divided into ∆X k i (k = 1, 2, . . . , N) in probability, and each subinterval is subjected to independent equal-probability sampling. During each deterministic calculation analysis, sampling is strictly guaranteed in each subinterval. To ensure that the extracted random numbers belong to each subinterval, the random number V i in the i-th subinterval must satisfy the following conditions: where i = 1, 2, . . . , N; V i is the random number of the first subinterval; and V is the random number with evenly distributed intervals. The LHS sampling design for n-dimensional random X i (i = 1, 2, . . . , N) variables involves the following two steps: (1) Divide each random variable into N equal-probability subintervals ∆X k i (k = 1, 2, . . . , N), and obtain the midpoint of each subinterval as a sample representative X k i of the corresponding random variable.
(2) For each random variable X i , extract a sample representative X k i and arrange it according to a random number, and arrange the samples of all the random variables according to random numbers, resulting in N random arrangements. Each sample contains a sample of all the random variables representing X k i .
In this study, the sampling parameters e pp , µ pp , e pw , µ pw , and Γ, which were to be calibrated, were sampled at 120 sample points for BP neural network training and verification. Thus, all the sampling points constitute a 5 × 120 matrix. Each row of the matrix represents a variable, and each column represents a set of test parameters, with the first 80 samples denoting the training data and the final 40 denoting the validation data. The virtual test data are presented in Table A1 of Appendix A.

BP Neural Network Structure Design
The BP neural network was proposed by Rumelhart and McClelland [31]; it is a multilayer feedforward neural network trained according to the error BP algorithm. It is extensively used in many fields as a supervised learning algorithm that maps a set of inputs to the accurate output [32,33]. This study is based on the Sigmoid weighting function and uses the mean square error (MSE) deviation function as the core content of the BP neural network. In the current study related to discrete element parameter calibration, the DEM parameter is the general input layer, whereas the angle of repose is the output layer. Based on the DEM batch simulation experiment presented in Section 2.3.2, this study establishes an approximate model of the relation between the microscopic parameters and the macroscopic angle response for each set of sample points. The topological structure of the BP neural network prediction model for the particle macro parameter is presented in Figure 10. In addition, the numerical fluctuation of static angle is measured because of the defects of the BP neural network itself. Therefore, the GA is used to adjust the weight and threshold of the network to reduce the influence of noise in the model. Finally, the behavior of the particles can be predicted for any set of DEM simulation parameters. The BP neural network prediction model is subsequently used to extend the finite numerical simulation database to "infinity."

Optimization of the BP Neural Network Performance on the Basis of the GA
While dealing with nonlinear problems, the BP neural network algorithm often converges to local solutions because of its serial search mechanism. To address this problem, the GA of "survival of the fittest in nature" is adopted [34]. Optimization to form a parallel search optimization method makes it easy to converge into the global sample space. While training the BP neural network, as depicted in Figure 11, the objective is to continuously optimize the weight and threshold of each layer such that the error between the final network output static angle value and the predicted output is sufficiently small in accordance with our expectations. The process of optimization is the process of constant "evolution" of weights and thresholds. During each iteration, the weights and thresholds are selected on the basis of the "natural environment," mutations, and crossovers and then leave excellent varieties, i.e., weights and offsets minimize output errors. Inferior varieties are eliminated, and, after multiple iterations, the "genes" of weights and thresholds improve, achieving

Optimization of the BP Neural Network Performance on the Basis of the GA
While dealing with nonlinear problems, the BP neural network algorithm often converges to local solutions because of its serial search mechanism. To address this problem, the GA of "survival of the fittest in nature" is adopted [34]. Optimization to form a parallel search optimization method makes it easy to converge into the global sample space. While training the BP neural network, as depicted in Figure 11, the objective is to continuously optimize the weight and threshold of each layer such that the error between the final network output static angle value and the predicted output is sufficiently small in accordance with our expectations. The process of optimization is the process of constant "evolution" of weights and thresholds. During each iteration, the weights and thresholds are selected on the basis of the "natural environment," mutations, and crossovers and then leave excellent varieties, i.e., weights and offsets minimize output errors. Inferior varieties are eliminated, and, after multiple iterations, the "genes" of weights and thresholds improve, achieving the training objectives. Furthermore, the GA-BP neural network is able to fit the data using fewer iterative steps. In this study, it can quickly find the input DEM parameters and optimize the target using the preset MSE values, thereby achieving global optimization faster.
While dealing with nonlinear problems, the BP neural network algorithm often converges to local solutions because of its serial search mechanism. To address this problem, the GA of "survival of the fittest in nature" is adopted [34]. Optimization to form a parallel search optimization method makes it easy to converge into the global sample space. While training the BP neural network, as depicted in Figure 11, the objective is to continuously optimize the weight and threshold of each layer such that the error between the final network output static angle value and the predicted output is sufficiently small in accordance with our expectations. The process of optimization is the process of constant "evolution" of weights and thresholds. During each iteration, the weights and thresholds are selected on the basis of the "natural environment," mutations, and crossovers and then leave excellent varieties, i.e., weights and offsets minimize output errors. Inferior varieties are eliminated, and, after multiple iterations, the "genes" of weights and thresholds improve, achieving the training objectives. Furthermore, the GA-BP neural network is able to fit the data using fewer iterative steps. In this study, it can quickly find the input DEM parameters and optimize the target using the preset MSE values, thereby achieving global optimization faster. In the BP neural network, the number of nodes in the input and output layers are determined, and the appropriateness of the number of nodes in the hidden layer is determined, which directly affects the performance of the network. In general, in case of a large sample space, a multilayer network is more accurate than a single-layer network; however, the time required for training considerably increases. An excessive number of hidden-layer nodes will increase the redundancy of the network, whereas the presence of very few hidden-layer nodes will increase the network error. No clear theoretical guidance is available for determining the number of hidden-layer nodes, which are generally determined by a combination of empirical formula methods and repeated test methods [35] and can be given as follows: where h is the number of hidden-layer nodes, m is the number of input layer nodes, and n is the number of output layer nodes, which is an adjustment constant between 1 and 10. The number of neurons in the hidden layer of the obtained Equation (6) is approximately 3-12, which is further explored below. Figure 12 denotes the relation between the number of neural network hidden-layer neurons determined by the BP neural network prediction model and the response-angle coefficient of determination R2. When the number of neurons in the hidden layer became seven, the coefficient of determination of the verification sample reached a maximum of 0.924. Therefore, seven hidden-layer neurons were selected to establish a BP neural network model for predicting the angle of repose. explored below. Figure 12 denotes the relation between the number of neural network hidden-layer neurons determined by the BP neural network prediction model and the response-angle coefficient of determination R2. When the number of neurons in the hidden layer became seven, the coefficient of determination of the verification sample reached a maximum of 0.924. Therefore, seven hidden-layer neurons were selected to establish a BP neural network model for predicting the angle of repose.

Prediction Results and Discussion
Among the 120 sets of sample spaces sampled by Latin hypercube, 80 sets of data were used in the BP neural network training model, whereas the remaining 40 groups were used in the post-training BP model to analyze the fitting accuracy of the BP model to the DEM model. The training data and forecast data used in this study are presented in Table A1 of Appendix A. The fitting accuracy of the BP model was evaluated by the coefficient R 2 determined in Equation (7). The closer the value of R 2 is to 1, the more accurate will be the BP model. In general, the approximate model and R 2 ≥ 0.9 are considered for engineering applications.
where y i is the first simulated value of the model, ŷ i is the first predicted value of the model, y i is the average of the simulated values, and N is the number of samples. The coefficient of Figure 12. Implied layer node performance test.

Prediction Results and Discussion
Among the 120 sets of sample spaces sampled by Latin hypercube, 80 sets of data were used in the BP neural network training model, whereas the remaining 40 groups were used in the post-training BP model to analyze the fitting accuracy of the BP model to the DEM model. The training data and forecast data used in this study are presented in Table A1 of Appendix A. The fitting accuracy of the BP model was evaluated by the coefficient R 2 determined in Equation (7). The closer the value of R 2 is to 1, the more accurate will be the BP model. In general, the approximate model and R 2 ≥ 0.9 are considered for engineering applications.
where y i is the first simulated value of the model,ŷ i is the first predicted value of the model, y i is the average of the simulated values, and N is the number of samples. The coefficient of determination R 2 based on the stacking angle of Equation (7) is 0.93. Figure 13 presents a comparison between the predicted and simulated values. determination R 2 based on the stacking angle of Equation (7) is 0.93. Figure 13 presents a comparison between the predicted and simulated values. As can be observed from Figure 13a, the evaluation coefficient is obtained as 0.924, which is calculated in accordance with Equation (7) by linearly fitting the relation between the DEM simulation value and the predicted value. However, Figure 13a also shows a small number of discrete points which can be attributed to the initial state and flow process of the sample powder particles to be tested during the simulation and to the poor reproducibility of the angle of repose caused by randomness. The predicted value of the BP model and the DEM simulation response angle are compared in Figure 13b. This discrete phenomenon can be observed only in a small amount, and the predicted value is close to the DEM simulation value, indicating that the established BP model is sufficiently reliable to be used for calibrating the DEM parameters. As can be observed from Figure 13a, the evaluation coefficient is obtained as 0.924, which is calculated in accordance with Equation (7) by linearly fitting the relation between the DEM simulation value and the predicted value. However, Figure 13a also shows a small number of discrete points which can be attributed to the initial state and flow process of the sample powder particles to be tested during the simulation and to the poor reproducibility of the angle of repose caused by randomness.
The predicted value of the BP model and the DEM simulation response angle are compared in Figure 13b. This discrete phenomenon can be observed only in a small amount, and the predicted value is close to the DEM simulation value, indicating that the established BP model is sufficiently reliable to be used for calibrating the DEM parameters.

BP Model Predictive Value of the Angle-of-Repose Search
To accurately calibrate the DEM parameters, the sample space was expanded to 300 sample points by LHS resampling and the mapping relation between the DEM parameters and the angle of repose was successfully established through the developed BP neural network (Figure 14a). The results denote a large dispersion and overlap in the output of the BP prediction model. To reasonably obtain a feasible solution that satisfies the mapping relation and to achieve a target value of the stacking angle mean of the actual test, which is listed in Appendix A, Table A2), the discrete outputs of the predicted output are compared and eliminated. The difference between the output value of the BP prediction model and that of the calibration experiment is less than 5% [19], and the condition is satisfied. These conditions can be expressed as: Subject to outputs = f (inputs), inputs min ≤ inputs ≤ inputs max . To quickly obtain a feasible solution that satisfies these conditions, a common loop structure algorithm is used to search for the distribution of the feasible domain space. The algorithm is presented in Figure 15. To quickly obtain a feasible solution that satisfies these conditions, a common loop structure algorithm is used to search for the distribution of the feasible domain space. The algorithm is presented in Figure 15.
The specific steps are as follows: (1) initialize the sample point i, and predict the output value AOR i = AOR 1 ; (2) identify whether the sample point i exceeds 300 sample points; (3) filter the predicted values of the BP prediction model according to the AOR Tol tolerance range and store them in a new sample space; (4) repeat steps (2) and (3) until 300 sample points have been searched, resulting in a new sample space as depicted in Figure 14b. To quickly obtain a feasible solution that satisfies these conditions, a common loop structure algorithm is used to search for the distribution of the feasible domain space. The algorithm is presented in Figure 15.

Mean-Shift Cluster Analysis
The mean-shift algorithm is a parameterless density estimation algorithm. The core idea of the mean-shift algorithm is to first calculate the offset mean of the current pixel, move the point to its offset mean, and continue to move as the starting point until the certain conditions are met and then stop iterating. Calculate the similarity between the cluster center candidate point and the candidate

Mean-Shift Cluster Analysis
The mean-shift algorithm is a parameterless density estimation algorithm. The core idea of the mean-shift algorithm is to first calculate the offset mean of the current pixel, move the point to its offset mean, and continue to move as the starting point until the certain conditions are met and then stop iterating. Calculate the similarity between the cluster center candidate point and the candidate area sample, and use the maximum value of the similarity function to find the target mean-shift vector, move the cluster center candidate point to the offset mean, and iterate until the cluster center candidate point No offset occurs. In the DEM parameter optimization, a feasible domain may comprise multiple subdomains exhibiting similar solution values, indicating similar particle macroscopic properties.
To effectively obtain the required DEM parameters, an appropriate algorithm is used to optimize the distribution and combination of feasible solutions. The classification and differentiation of data with approximate features using mean-shift clustering analysis [36] is a hill-climbing algorithm based on kernel density estimation, which is a clustering algorithm that does not require the number of clusters to be specified. Processing [37] and cluster analysis [38] are the commonly used cluster analysis methods in the field of data mining. The purpose of mean-shift clustering is to find the data sample points of the same cluster along the direction of increasing density. Further, the class with the highest frequency of sample point access is marked as the class of the current point set. The clustering results of this study for the optimized sample points (Figure 14b) are listed in Table 2.

DEM Parameter Verification Under Quasi-Static Conditions
Under quasi-static conditions, the verification of DEM input parameters is the most basic calibration method. The cumulative lift angle of polyurethane powder particles was verified by the baffle lift test commonly used for static angle of repose. The test apparatus is presented in Figure 16f; the square casing was made of 60 × 60 × 80 mm 3 Plexiglas. During the test, a sample material with a mass of 15 g was placed in the square shell. When the baffle was lifted at a constant speed of 0.1 m/s, the flow of the PU powder resulted in the formation of a pile; after the particles of the population were stabilized, the formed powder particles were inclined and horizontal. The angle between the two is the static angle of repose of the powder. The simulation settings were consistent with the experimental conditions, and five sets of calibration parameters optimized in Section 3.5.2 were used as the inputs. The test simulation results are denoted in Figure 16a

Dynamic Condition DEM Parameter Verification
In order to further verify the DEM parameter calibration results, a more realistic dynamic angle of repose of the PU particles was verified via a rotary drum test. The experimental setup is shown in Figure 17. The cylinder is made of Plexiglas with an inner diameter and length of 45 and 100 mm, respectively. Prior to the experiment, the particle surface was set to a level with a particle loading of 30%, a rotational speed of 60 rpm, and a clockwise rotation for 5 s. After standing for 2 s, the excessive porosity was eliminated; after the particles were stabilized, the angle between the surface of the particles and the horizontal plane was measured as the dynamic angle of repose. The five sets of calibration parameters in Table 2 are used as inputs, and the simulation settings are consistent with the experimental conditions. The test results are shown in Figure 18. Small errors indicate the validity of the calibration results.

Dynamic Condition DEM Parameter Verification
In order to further verify the DEM parameter calibration results, a more realistic dynamic angle of repose of the PU particles was verified via a rotary drum test. The experimental setup is shown in Figure 17. The cylinder is made of Plexiglas with an inner diameter and length of 45 and 100 mm, respectively. Prior to the experiment, the particle surface was set to a level with a particle loading of 30%, a rotational speed of 60 rpm, and a clockwise rotation for 5 s. After standing for 2 s, the excessive porosity was eliminated; after the particles were stabilized, the angle between the surface of the particles and the horizontal plane was measured as the dynamic angle of repose. The five sets of calibration parameters in Table 2 are used as inputs, and the simulation settings are consistent with the experimental conditions. The test results are shown in Figure 18. Small errors indicate the validity of the calibration results. 30%, a rotational speed of 60 rpm, and a clockwise rotation for 5 s. After standing for 2 s, the excessive porosity was eliminated; after the particles were stabilized, the angle between the surface of the particles and the horizontal plane was measured as the dynamic angle of repose. The five sets of calibration parameters in Table 2 are used as inputs, and the simulation settings are consistent with the experimental conditions. The test results are shown in Figure 18. Small errors indicate the validity of the calibration results.

Results and Discussion
In this section, mean-shift clustering analysis is used to obtain five sets of DEM parameter sets input into the baffle lift and drum device for dynamic and static angle verification, respectively. Table 3 gives the relative error between the response angles under different load conditions. It is found that the DEM parameter set measuring the dynamic angle of repose using the funnel flow test device in Section 2.1 does not meet the static angle of repose verification (22.18%, 20.29%, 17.57%, 21.34%, 16.11%). Such a large error is due to the fact that the powder is subjected to different conditions of the external load and the response-angle deviation is generated, which means that the material calibration test and verification examples under dynamic/static conditions need to be performed separately.

Results and Discussion
In this section, mean-shift clustering analysis is used to obtain five sets of DEM parameter sets input into the baffle lift and drum device for dynamic and static angle verification, respectively. Table 3 gives the relative error between the response angles under different load conditions. It is found that the DEM parameter set measuring the dynamic angle of repose using the funnel flow test device in Section 2.1 does not meet the static angle of repose verification (22.18%, 20.29%, 17.57%, 21.34%, 16.11%). Such a large error is due to the fact that the powder is subjected to different conditions of the external load and the response-angle deviation is generated, which means that the material calibration test and verification examples under dynamic/static conditions need to be performed separately. Whether under static (baffle lift test) or dynamic conditions (rotary drum test), the error between the comparison tests is lower than 5%, which is also in line with the error criteria mentioned by Ye et al. [19]. Simultaneously, we found that although the relative error verified by the baffle lift test is small, it is between 0.77% and 4.86%. However, the simulation angles (37.

Conclusions
The JKR model exhibiting adhesion characteristics was used for an ultrafine agglomerated powder (120-mesh PU powder) that is difficult to calibrate, and the sample material was characterized. A DEM model was used based on the improved FBS-104 powder measuring instrument. In the parameter calibration method, the GA-BP optimization algorithm was used during the calibration test to reduce the number of simulations required for optimization. The results of this study can be summarized as follows: • To achieve the DEM parameter calibration of the micron-scale material samples, the concept of cohesion number proposed by Bharadwaj et al. [23] for coarse-grained discrete element modeling can improve the calibration efficiency without affecting the calibration accuracy.

•
In the DEM parameter calibration process, the BP approximation model modified using the GA can easily converge into the global sample space to avoid local optimal solutions, which would result in large errors in the calibration results.

•
The study found that the rotating drum as a verification example of dynamic conditions can be used to verify the DEM parameter set of the funnel flow meter more accurately than the baffle lift test under static conditions and that the relative error is 0.63%-2.72%.

•
The intrinsic parameters of the particles were explored during the simulation. The small range of the intrinsic parameters of the particles did not strongly influence the angle of repose; however, reducing the shear modulus was observed to considerably reduce the calculation cost.

•
To avoid large dispersion and overlap between the predicted value and the simulated target value, a combination of the cyclic search algorithm and mean-shift cluster analysis can be used to find the optimal solution for each subclass in the sample space.

•
Although the ultrafine agglomerated PU powder exhibited poor reproducibility with respect to the angle of repose during the stacking test, the baffle lifting method verified that the predicted output value of the BP model after training agreed with the angle obtained in the actual stacking test.

•
During this calibration process, we observed that the DEM parameter calibration results are not a set of solution sets; each output value (angle) has multiple sets of mapping relations with the input DEM parameters. This also revealed the calibration work of the DEM parameter, which is a search process of the parameter set "best combination." Based on a test design method, AOR was used as the target response value to find a suitable mapping relation for application in this study.

•
In the future, a large number of variables (bulk density, time step) can be added as the response value of the output layer of the prediction model, and the influence of the DEM parameters on the angle of repose can be further investigated.  The average value of the static angle of the bulk accumulation experiment is 47.801 • .