Avoiding ambiguity in DEM in-situ calibration for dry bulk materials

Minerals Engineering 145 (2020) 106094

Contents lists available at ScienceDirect

Minerals Engineering journal homepage: www.elsevier.com/locate/mineng

Avoiding ambiguity in DEM in-situ calibration for dry bulk materials

T

Stefan Kirsch TU Dresden, Germany Robert Bosch Packaging Technology, the Netherlands

A R T I C LE I N FO

A B S T R A C T

Keywords: Bulk good handling Discrete element method Model calibration Vertical filling Numerical uncertainty

Discrete element simulation (DEM) of the behavior of particulate solids experiences a continuous growth in popularity. However, simulations often require an iterative approach for the determination of the material dependent contact parameters, which often requires a significant portion of the time and budget of a study. Several more efficient procedures for model calibration have been suggested in recent times, but most of them struggle to offer a unique solution to the calibration problem. In this study, a calibration procedure that seeks to avoid ambiguity by attempting to reproduce different time signals instead of a single target parameter is presented. The calibration was performed in a drop test and the search for a suitable parameter set was performed on a meta model obtained from a regression of samples of the parameter space. The calibrated models were then validated by comparison to experimental data. The results show that not all time signals used are equally suitable for the task and that there is some inconsistency in the predictions of the models. However, under the right circumstances, the models are capable of reproducing industrially relevant process parameters. The method can be easily adapted to other industrial processes where time dependent data is available.

1. Introduction Prediction and control of the flow behavior of bulk solids is an important competence for several industry fields, including agriculture, packaging and mining. In the past decades, an increase in computational power and expanding knowledge in Discrete Element Methods (DEM) have allowed virtual representation of many processes, reducing cost for prototyping and allowing additional insight. Some processes that have been described numerically with the DEM include bulk transfer via conveyors (Grima et al., 2011), chutes and funnels (Cleary and Sawley, 2002; Grima et al., 2010), comminution (Morrison and Cleary, 2008; Cleary and Morrison, 2011; Djordjevic et al., 2003) and separation (Harzanagh et al., 2018). The validity of DEM simulations strongly depends on the choice of contact parameters which need to be specified by the user. While there are some methods to measure contact parameters, most applications rely, at least partly, on iteratively tuning of the parameters to an experimental reference (model calibration) (Barrios et al., 2013; Coetzee and Nel, 2014; Coetzee, 2016; Coetzee, 2017; Katterfeld et al., 2019). This procedure can however be very time consuming, since every parameter set tested requires one DEM solver call. Rackl and Hanley (2017) tackled this issue by partly performing the model calibration on a computationally cheap surrogate model. They obtained the model from a Kriging regression of a set of predetermined solver runs. This

approach was also shown to handle uncertainty in the solver behavior when the model was calibrated in a realistic process with random variation (in-situ calibration) (Kirsch, 2018). A shortcoming of most calibration approaches is that they rely on a single scalar experimental reference, which can lead to several ambiguous parameter combinations that provide an equally good fit with the reference (Wensrich and Katterfeld, 2012; Roessler et al., 2019). It has been shown that this can be avoided by calibrating the model to a set of target parameters, which significantly reduces the ambiguity of the solution (Roessler et al., 2019). For this study, it was reasoned that calibrating the model parameters to an experimental time-variant parameter, instead of a single scalar target parameter, would reduce ambiguity in a similar way. The best parameter set is then the one that provides the best match between simulation and experiment over the entire time range. For this study, vertical filling was chosen as an example process for in-situ calibration. Several time signals were extracted from a drop experiment and the calibration was repeated with all of them for two different realistic example goods. Finally, the fidelity of the calibrated models was evaluated in a validation step with two different experimental scenarios with the bulk good.

E-mail address: [email protected]. https://doi.org/10.1016/j.mineng.2019.106094 Received 11 July 2019; Received in revised form 14 October 2019; Accepted 15 October 2019 0892-6875/ © 2019 Elsevier Ltd. All rights reserved.

Minerals Engineering 145 (2020) 106094

S. Kirsch

Nomenclature CoD CoP d → e ∊ Fn Fτ Fτ , el fmeta fsim Knl Knu LHC l μ μd μs PB PCC PP

p RR r ri

(a) (b)

coefficient of determination coefficient of prognosis particle diameter local residual of regression restitution coefficient normal force tangential force elastic tangential force meta model function solver function normal contact stiffness for loading normal contact stiffness for unloading Latin Hypercube particle length friction coefficient dynamic friction coefficient static friction coefficient between particles and boundary Pearson’s correlation coefficient between particles number of dimensions rolling resistance seed for random numbers radius of particle i

(c) SSEPrediction

SST Δs → si sn sτ t Δt t ftl tltl R (tres ) → u → û → w

→ ŵ → xc → xp → x p, opt → y

2. Process simulation

The discrete element method (DEM) is capable of simulating a wide range of processes involving granular matter (Zhu et al., 2008). The most popular implementation is the force-based rigid-body approach, where the particles are represented by perfectly stiff bodies (Schneider, 2012). As long as no contacts take place, the particles’ trajectories follow Newtonian motion over a set of discrete time steps with length Δt . Via the location of the particles’ center of mass and their extension in space, particle contacts can be detected. In the simplest case of spheres, a contact between two particles is given if

As an example process for the in-situ calibration, a vertical filling process, which is commonly found in the packaging industry, was chosen for this study (Frank and Holzweißig, 2016; Kirsch and Philipp, 2018; Kirsch, 2018). A typical design is shown schematically in Fig. 2f. The goal of the process is to drop a defined amount of bulk material into a container that must not be sealed before the last particle has passed the container seam (in this case a bag made from plastic film). Each fill cycle includes the events presented in Fig. 3. The time stamps t ftl and tltl are the residence times of the first and last particle to leave the end of the tube. The difference between them is the residence time range R (tres ) and determines how much time is necessary for the filling of a bag and thus determines the maximum bag production rate. Since the trajectories of individual particle’s are highly prone to randomness, R (tres ) varies between cycles. Thus, the time allowed for filling must be chosen so that it can accommodate the variation.

(1)

(2)

Fτ (t ) = min (|Fτ , el |, μFn (t )) With Fτ , el = Fτ (t − Δt ) − K τ Δs τ

(3)

estimate for solver response initial conditions DEM contact parameters optimal DEM contact parameters meta model response

2.2. Reference filling process

→ → where s1 and s2 are the locations of the particles and r1 and r2 are their radii. The DEM allows particles to overlap, which mimics deformation upon contact. The repulsive forces between particles at a given time step t depend on the current magnitude of overlap. Various contact models are used in the different DEM implementations available. For the present study, we used the commercial software package Rocky DEM’s implementation of the linear hysteresis model. This model considers elastic and plastic deformation for the calculation of contact forces. This is achieved by a higher stiffness constant during unloading, i.e. the portion of the contact event during which the particles separate (see Fig. 1) (Walton et al., 1986; Luding, 1998). The laws for the normal and the tangential forces Fn and Fτ at time step t are given in Eqs. (2) and (3) (ESSS, 2019). Fn (t ) = Knl sn (t ) Loading Fn (t ) = Fn (t − Δt ) + Knu Δsn Unloading

estimate for experimental observation solver response

loading and unloading respectively. Their ratio determines the energy dissipation of contacts, which depends on the restitution coefficient ∊. sn and s τ are the accumulated overlaps in normal and tangential direction at the respective time step t with the change in overlap Δs = s (t ) − s (t − Δt ) since the previous time stamp.

2.1. Particle contact model

||→ s1 − → s2 || − r1 − r2 < = 0

visible area covered by the bulk good over time location of the virtual center of mass of visible bulk surface location of portion front (c)1 and tail (c)2 estimated amount of solver variation captured by meta model total variation of solver response change in overlap location of particle i normal particle overlap tangential particle overlap DEM time step or time stamp length of DEM time step residence time of the first particle to leave the tube residence time of the last particle to leave the tube residence time range experimental observation

Fig. 1. Force-overlap relationship in the linear hysteresis model.

μ is the friction coefficient. Knl and Knu are the contact stiffnesses for 2

Minerals Engineering 145 (2020) 106094

S. Kirsch

to the partially random initial orientation of the particles in the collection bucket before the drop. To account for this in the simulations, analogously to Kirsch (2018), different randomized particle beds were generated in Rocky DEM. For each solver run, one of them was then chosen randomly as the initial state. This results in a noisy solver i.e.the w depends not only on the chosen input parameters, solver response → but also on a random seed r:

→ u ̂= → w = fsim (→ xc , → xp, r )

(5)

2.4. Ambiguity in DEM model calibration For DEM model calibration, the usual approach is to chose a single scalar reference parameter which is obtained from a reference trial. Then the contact parameters → x p are varied to determine the parameter set → x p, opt with which the simulations reproduce the reference the closest. Typically, reference parameters are the dynamic or static angle of repose (Gröger and Katterfeld, 2006; González-Montellano et al., 2012; Marigo and Stitt, 2015; Liu et al., 2016; Rackl and Grötsch, 2018) or discharge parameters (Kirsch, 2018; Roessler et al., 2019; Markauskas and Kačianauskas, 2010). However, it has been shown that several parameter combinations may allow to reproduce the experimental reference value. This could happen for instance, if several parameters influence the value of the target parameter in a similar way, and thus can compensate for each other. (Roessler et al., 2019; Wensrich and Katterfeld, 2012). It has been shown that a unique solution can however be found by using several reference parameters obtained from separate calibration tests. Furthermore, a consideration of the measurement uncertainty might be required to obtain the full range of possibly suitable parameter combinations. (Roessler et al., 2019). Additionally, it is important to consider whether the chosen calibration tests or, more specifically, the chosen calibration target parameters are sufficiently sensitive to the parameters to be calibrated (Rackl and Hanley, 2017; Coetzee, 2016). Sensitivity can be judged by varying the input parameters → x p of the DEM simulation and tracing the w . Then, the linear correlation between changes in the solver response → individual input and an output parameter can be expressed via the Pearson coefficient of correlation PCC. Large values of |PCC| indicate a strong influence of the respective parameter on the solver behavior (Rogers and Nicewander, 1988). It is, however, not necessary, that the target parameter is sensitive to all calibration parameters. Instead, it should be made sure that the productive process to be modeled (in this case the vertical filling process) is sensitive to the same input contact parameters as the calibration test. All insensitive parameters can then be chosen arbitrarily within certain bounds. More specifically, this means that both calibration test and actual process should share the same sensitivity pattern to the input parameters. A straightforward way to ensure identical sensitivity is to calibrate the process in the actual process itself. In the case of the vertical filling process, the range of the residence time R (tres ) determines overall cycle duration, so it is the most obvious target parameter. However, in theory, every parameter that can be measured could also be used as a calibration target.

Fig. 2. Schematic view of filling in a vertical tubular bagger. Adapted from Kirsch (2018)).

Fig. 3. Timeline of one cycle of the vertical filling process and the relevant time parameters.

2.3. Numerical uncertainty The DEM simulations aim to approximate an experimental obw . The solver response depends on servation → u with the solver response → x c (such as equipment geothe input parameters for initial conditions → w is the metry and particle location) and a contact parameter set → xp. → output of the solver function fsim and is thus an estimate for the ex→ perimental outcome u . (Kirsch, 2018; Kleijnen, 2017)

→ u ̂= → w = fsim (→ xc , → x p)

2.5. Calibration on a meta model In the meta model approach, the parameter space is sampled with n parameter sets. There, the solver is executed yielding the outputs → wi = 1,2, … , n . In the present case, the solver response is the deviation between the DEM simulation’s prediction at sample i and the experimental reference. In the next step, the samples are used as anchor points for a y = fmeta (→ x p) is in turn an approxregression. The regressive function → → ̂ imation w for the solver response. (Kleijnen, 2017)

(4)

Note that Eq. (4) represents a deterministic solver, that produces a conx c and → stant output, as long as → x p are constant. The industrial filing process exhibits some stochastic variation, due

→ y =→ w ̂ = fmeta (→ x p) with → x c = const . 3

(6)

Minerals Engineering 145 (2020) 106094

S. Kirsch

camera aimed at the falling tube at 1024 frames/s. The reference videos were processed with the Matlab® Image Processing ToolboxTM, so that the bulk good was separated from the background (Fig. 5). This allowed for the automatic extraction of positional data over time. In this way, three time signals were obtained from every reference video:

The meta model is much cheaper to evaluate than the DEM solver, which significantly speeds up the search for the parameter set at which the lowest deviation from the experimental reference is achieved (Rackl and Hanley, 2017; Barton and Meckesheimer, 2006). x p) and the meta The difference between the solver response fsim (→ x p, i ∈ → x p in the parameter space is called x p) at any point → model fmeta (→ the local residual → ei .

→ ei = fsim (→ x ) − fmeta (→ x)

(a) The visible area covered by the bulk good (b) The position of the virtual center of mass over the visible particle surface (c) The position of the bulk portion’s front (c)1 and tail (c)2, i.e. the respective position of the leading and the trailing particle

(7)

The quality of the meta model, i.e. its accuracy in predicting the solver w , can be estimated by correlation measures between the response → solver response wi and the meta model response → yi . If linear correlation is to be examined, Pearson’s correlation coefficient PCC can be used (Rogers and Nicewander, 1988). Correlation measures for more sophisticated regression models include the coefficient of determination CoD for polynomial fit (Montgomery et al., 2003). If, on the other hand, meta model quality should be judged independently of the regression method, the from the CoD derived coefficient of prognosis CoP can be used. The CoP was introduced by Most and Will (2008, 2011).

CoP = 1 −

SSEPrediction SST

The area covered by the bulk good was normalized to the maximum area that could be covered if the entire tube was filled. The center of mass was normalized to the tube’s height (one being the top and 0 being the bottom of the tube). The center of mass exhibits large variation between runs as soon as the majority of the bulk has left the tube at the bottom end. It was therefore cut off after the time stamp where the point wise-standard deviation between runs showed a stark increase. For the Jelly Beans, this was at 610 ms and for Penne this was at 640 ms. The portion front (c)1 and tail (c)2 were also normalized. Note that the time stamps where the bulk portion’s front and tail reach the bottom of the tube are equivalent to the time stamps t ftl and tltl respectively discussed in Section 2.2. Consequently, their difference is equivalent to the residence time range R (tres ) . The latter has already been shown to be a suitable parameter for DEM model calibration (Kirsch, 2018). The drop experiments were each conducted at least 10 times. The final curves were obtained by point-wise averaging the curves from the individual runs. The setup for the funnel discharge tests was developed by Kirsch (2018) and is shown schematically in Fig. 6. The experiments were carried out by filling the transparent polycarbonate funnel with bulk good and flattening the fill line horizontally. Then, the slider at the bottom was removed and the discharge was filmed at 100 frames/s. The footage was processed analogously to the drop test in Matlab®, so that the bulk good was separated from the background (Fig. 7). The visible surface of the bulk good was tracked over time, normalized to the first value before discharge was initiated. The experiment was repeated three times and the resulting curves were averaged point-wise. For the simulations, the experimental setups were reproduced virtually in the DEM environment Rocky DEM. For the drop test, the

(8)

SST is the total variation of the solver response → x p . SSEPrediction is an estimate of how much of the variation is captured by the meta model and is calculated in a series of cross-validation steps. In each step, the error between the solver response wi at the parameter set → x p, i and the meta x p, i ) is calculated for the case that sample i was not used for model fmeta (→ the regression. The CoP approaches one if the meta model still makes accurate predictions, when individual anchor points are omitted from the regression. (Most and Will, 2008, 2011) A noisy solver (see Section 2.3) has some non-trivial implications for the problem formulation of the calibration and the notion of an optimal parameter set → x p, opt . A straight-forward requirement of the optimal parameter set is the one that reproduces the mean solver response with → x p = const . the closest (Kleijnen, 2017). The meta model attempts to predict the expected value of the solver response (i.e. the mean over an infinite number of runs). A benefit of the meta model approach is that this way, some noise can be ”smoothed out” (Barton and Meckesheimer, 2006). Adding samples allows for a more detailed approximation of complex solver behavior, while further reducing noise (Kleijnen, 2017). A criterion that can be used to determine whether enough anchor points were created to reach satisfactory meta model fidelity, is the CoP. However, the global CoP tends to underestimate local meta model fidelity in the case of non-uniform sampling density (Liu et al., 2018). A suitable stop criterion is the → stagnation of the residuals e after adding anchor points over several iterations (Kirsch, 2018). In a similar fashion, we reasoned that stagnation of the CoP (rather than it’s current value) indicates that the maximum meta model quality has been reached. 3. Parameter identification 3.1. Reference trials Experiments were carried out in two different test setups to obtain reference data of the true bulk good behavior. The first one was a simplified drop test setup as developed by Kirsch and Philipp (2018). The second one was a funnel discharge setup, variations of which are commonly used for DEM model calibration (Frank and Holzweißig, 2016; Roessler et al., 2019; Coetzee, 2017). The simplified drop setup is presented in Fig. 4. In order to allow observation of the bulk’s behavior, the setup was constructed from transparent polycarbonate. The sample bucket’s bottom consisted of servo-driven flaps to release the bulk good. The drop tests were filmed with a high-speed

Fig. 4. Experimental and measurement setup of the drop test with high speed camera. Measures shown in mm. Adapted from Kirsch and Philipp (2018). 4

Minerals Engineering 145 (2020) 106094

S. Kirsch

t = 0.0s

t = 0.5s

Fig. 5. Screenshots of videos of the experimental setup and the result of the binarization. The asterisk marks the current virtual center of mass.

t = 1.0s Fig. 7. Screenshots of videos of the experimental funnel setup. The detected surface of the bulk good is shaded in lighter gray.

3.2. Parameter variation The solver output was defined as the deviation between the experimental reference curve and the simulated curve. Both curves are affected by measurement uncertainty due to randomness, so it is arguably not desirable to aim to perfectly reproduce the reference curves. Instead, the root mean square error (RMSE) between the curves was used as the error function. This allowed to put heavier weight on large local deviations while putting less weight on small ones. The calibration’s goal is thus to find the parameter set → x p, opt where the RMSE between the simulated and experimental curves becomes minimal. With this approach, it was reasoned to be able to achieve a strong physical foundation for parameter sets where the error function becomes low, despite the randomness intrinsic to the drop process. Three different objective functions were obtained from the time signals measured, all of which were formulated as separate optimization problems. The RMSE of the location of the particle tail and front were summarized into a sum optimization criterion with weight 1. A x p were considered for the calibratotal of seven parameters x p,1 … 7 ∈ → tion. A list of the parameters and their considered ranges are shown in Table 2. The parameter variation and optimization where conducted in the design of experiment and optimization package Optislang® (Dynardo

Fig. 6. Experimental and measurement setup of funnel discharge with high speed camera. Measures shown in mm. Adapted from Kirsch and Philipp (2018). Depth = 100 mm.

measured motion profile of the flaps was included. Of every simulation, a video was exported and processed analogously to the videos of the experiment. The drop tests were carried out with Jelly Bean candy and Penne pasta as two different sample bulk goods. The funnel discharge test was performed with the Jelly Beans. For the virtual representation of the materials, the parameters presented in Table 1 where used in Rocky DEM. The Young’s modulus was chosen high enough so that impractical particle overlap well above 1% during collisions was prevented (Paulick et al., 2015). The shapes of the particles were approximated by spherocylinders (pill shapes) to save numerical cost. The virtual representations of the particles are illustrated in Fig. 8. The dimensions of the virtual particles were chosen so that the volume was represented accurately. The particle density, however, was adjusted so that the virtual bulk good’s bulk density was identical to the real good. The experiments and simulations of the drop test were conducted with 700 grams of the respective bulk good. The funnel discharge experiments and simulations were conducted with 2000 grams of bulk good.

Table 1 Constant parameters. Parameter Jelly Beans

Pasta (Penne)

Material

Value

Young’s Modulus

107 Pa

Bulk density

919. 0 kg/m3 8.5 mm 16.1 mm

Particle diameter d Particle length l Young’s Modulus Bulk density

Polycarbonate

5

Particle diameter d Particle length l Young’s Modulus

107 Pa 367.3 kg/m3 8.9 mm 39.4 mm

107 Pa

Minerals Engineering 145 (2020) 106094

S. Kirsch

performed for each of the two sample goods. Consequently, we obtained six different meta models specific to one objective function and one bulk good. In total, about 950 solver calls were performed per meta model. After the fifth refinement iteration, the objective function was minimized on the final meta model. Optimization was performed with the Evolutionary Algorithm available in Optislang® (Dynardo GmbH, 2019). The process was performed six times on the respective meta models with the respective objective function. The resulting parameter sets were then validated by using them as the input for the DEM solver. Validation was performed by simulating both experimental setups discussed in 3.1. To average out variability, the validation runs were performed with 50 and 10 solver calls for drop test and funnel discharge respectively. The point-wise mean curves obtained were once again compared to the reference curves from the experiments. Additionally, in the drop test, the accuracy of the prediction of the mean drop time range R (tres ) was evaluated for each of the six cases.

Fig. 8. Comparison between the real and the virtual particles.

Table 2 Calibration parameters → x p and corresponding range. Parameter

Material

Friction

Particles - Particles Particles - Boundary

Restitution Rolling resistance

Particles - Particles Particles - Boundary –

Case

Symbol

Range

Static

μs, PP

0 …1

Dynamic

μd, PP

0 …1

Static

μs, PB

0 …1

Dynamic

μd, PB

0 …1

∊PP ∊PB RR

0.1 …1 0.1 …1 0 …1

– – –

GmbH, 2019). The tool allowed to automate most of the following tasks. The parameter sets from which the meta model was constructed were determined by adaptive sampling, which will be described in the following: After a relatively coarse exploration sampling phase in the whole parameter space, further samples were added only in regions of low predicted errors. This allowed creating higher sample density in the more promising areas of the meta model, and thus allowed a finer distinction between small changes (Liu et al., 2018). The refinement phase was implemented similar to Kirsch (2018) according to Fig. 9. Before every iteration, the parameter space was sampled within the current set bounds with Latin Hypercube sampling (LHC) (McKay et al., 1979; Helton and Davis, 2003). During the exploration phase, the solver was evaluated once at every sample for more than 200 sample points. Following this, refinement was performed over several iterations with the adaptive meta modeling functionality available in Optislang® (Most, 2019). Each refinement iteration consisted of the following three steps: Pre-optimization, Sample addition and Meta model construction. In the first step, a set of parameter combinations on the current meta model was chosen to serve as starting points for the pre-optimization. This yielded a set of local optima which are clustered into regions where the meta model predicts high DEM model fidelity. In the third step, the overall parameter space is constrained to the promising regions where further samples are then added with LHC sampling. Following this, the solver is executed at the new samples. Lastly, a temporary meta model is constructed by Kriging regression from the existing anchor points and the current CoP is calculated (Most, 2019; Jones et al., 1998). During refinement, the solver was executed three times per sample to average out some of the variability of the solver behavior and increase meta model fidelity in the promising regions (Kleijnen, 2017). The meta model approach is capable of ”smoothing out” the noisy solver response to obtain a better estimate for the expected value (i.e. the mean of infinite runs). However, if only a low number of samples are available for regression, there is a considerable risk of over-fitting outliers. The more samples are used for the regression, the lower the chance of over-fitting. Since the CoP is sensitive to outliers, it is a suitable measure for the current quality of the model in noisy environments (Most and Will, 2008). In order to remain economical with solver calls, the refinement phase was canceled after five iterations. After this, the average CoP either stagnated or did not improve by more than 2% per iteration. The locations of the samples added during refinement depend on the optima found in the pre-optimization steps, which in turn depend on the chosen objective function. The sampling and refinement was therefore performed for each of the three objective functions (a)–(c), each yielding a specialized meta model. The entire process was

4. Results and discussion 4.1. Uniqueness of the optima Tables 3 and 4 show the values of the contact parameters for Jelly Beans and Penne identified as optimal according to the respective target parameter (a)–(c) used for calibration. Table 5 lists the result of the sensitivity study for Jelly Beans as Pearson’s correlation coefficients PCC of the contact parameters. The values of the parameters with the greatest linear correlation, i.e. the greatest impact on the respective target parameter, are printed in bold. For the calibrated parameter values, we find that the different calibration targets yield different results for optimal parameter combination. For Jelly Beans, the calibration runs however agree on the value of μd, PB = 0 and RR = 0 . Table 5 indicates that residence time range R (tres ) , i.e. the performance parameter for the industrial process, is mostly sensitive to the dynamic friction coefficients μd, PB and μd, PP . In fact, this is true for all calibration targets (a)–(c) as well as for the funnel discharge. However, while R (tres ) is correlated to μd, PB by 0.8, for (a), (b) and funnel discharge, we find a roughly 30% lower sensitivity to μd, PB . Instead (a), (b) and funnel discharge are mostly sensitive to μd, PP by >0.85. In contrast, the portion tail (c)2 shares an almost identical sensitivity pattern with

Fig. 9. Workflow for DEM input contact parameter calibration (Adapted from Kirsch, 2018). 6

Minerals Engineering 145 (2020) 106094

S. Kirsch

4.2. Model validation

Table 3 Calibrated parameters → x p for Jelly Beans.

Assessing the meta models alone does not allow for a final conclusion on whether the parameter sets identified by the optimizer are actually capable of making accurate predictions. This is due to two

Jelly Beans Values

(a)

(b)

(c)

μs, PP

0.410

1

0.801

μd, PP

0.106

0

0.493

μs, PB

1

0.744

0

0.7

0

0

0

0.6

0.208 0.521 0

0.100 0.747 0

0.517 0.579 0

0.5 d,PB

1.6

[-]

μd, PB

∊PP ∊PB RR

(a) RMSE surface [-] (Jelly beans)beans) (a) RMSEvisible visible surface [-] (Jelly

0.4 4 1.

0.3

(b)

(c)

μs, PP

0.208

0.463

1

μd, PP

0.124

0

0

μs, PB

0.231

0.978

0.872

μd, PB

0.201

0.229

0.397

∊PP ∊PB RR

0.644 0.499 0

1 0.802 0

0.356 0.720 0.571

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

[-]

d,PP

(b) RMSE of mass [-] (Jelly beans)beans) (b) RMSEcenter center of mass [-] (Jelly

0.7 0.6 0.5 d,PB

[-]

(a)

2

1

Values

1.

0.3

0.1

Penne

0.4

0.7

0.2

5 0.

Table 4 → Calibrated parameters x p for Penne.

0.4 1.2

Table 5 Sensitivity between input parameters → x p and respective process target parameter for Jelly Beans. Highlighted values are the ones with the greatest linear correlation.

0.3 1

0.

0.3

0.1

2

0. 1

Jelly Beans

0.8

0.5

0.2

0

Funnel

(a)

(b)

(c)1

(c)2

μd, PB

0.80

0.42

0.61

0.59

0.63

0.81

μd, PP

0.40

0.90

0.86

0.87

0.34

0.38

μs, PB

0.25

0.19

0.04

0.19

0.16

0.23

RR μs, PP

0.23 0.15

0.22 0.31

0.55 0.32

0.28 0.25

0.18 0.17

0.22 0.13

∊PB ∊PP

0.02 0.00

0.07 0.02

0.20 0.41

0.09 0.36

0.03 0.03

0.00 0.04

0

0.2

0.4

0.6

0.7

0.8

[-]

d,PP

(c) portion frontfront [-] (Jelly beans)beans) (c)11RMSE RMSE portion [-] (Jelly

0.07 0.0 6

0.6

0.0 5

0.5

[-]

R (tres )

0.4 0.3

3

0.2 0.1 0

0.0

0.015

0.0 25

0.1

0.2

0.02

0.3

0.04

0.4 d,PP

0.7

[-]

0.7

0.8

0.9

0.7

0.9

d,PB

0.6

[-]

0.6

0. 55

0.5

0.4

0.4 0.3

0.3 0.2

0.2

0.15

0.1

0.1 0

0.5

(c) 2 RMSE portion tail [-] beans)beans) (c) RMSE portion tail(Jelly [-] (Jelly 2

0.6 0.5

0.04

0.7

R (tres ) . As discussed in Section 2.4, the higher the correlation between a calibration target parameter, the higher the certainty at which a given input parameter can be determined. With this reasoning, (c)2 would be a better choice over (a) and (b) in determining the most important input parameter for R (tres ) . Figs. 10 and 11 show 2D-projections of the meta models of the error functions (i.e. the RMSE between the curves from the simulation and the experiment). The projection was constructed by setting the parameters with the greatest mean impact (i.e. μd, PP and μd, PB ) as variable and leaving all other input parameters constant. The constant parameters were set to their respective optimum i.e. where a sample predicted the lowest RMSE. The areas separated by isolines are shaded according to value. The region with he lowest values are colored black. A single continuous region with the lowest error gives a first indication for a unique minimum, and thus an unambiguous solution. We find such a region in all meta models, except for Penne with case (c)1 (portion front). However, the location of the portion front (c)1 and tail (c)2 are considered in a single objective function. Since the predicted error for the portion front is overall much lower than for the portion tail, the location of the optimum is dominated by the latter. Consequently, this should not be an issue for the uniqueness of the optimum, since the meta model (c)2 has a single low-error region.

5 0.03

5 0.01

d,PB

∣PCC∣

0.05

0.1

0.2

0.3

0.4 d,PP

0.5

0.6

0.7

0.8

0.9

[-]

Fig. 10. 2D-projection of meta models for Jelly Beans. 7

Minerals Engineering 145 (2020) 106094

S. Kirsch

(a) RMSE surface [-] (Penne) (a) RMSEvisible visible surface [-] (Penne)

0.7 0.6 0.5 0.4

1

2

2.5

3

0.3

1.5

d,PB

[-]

1.2

0.2

2.5

0.1 0.1

0.2

0.3

0.4

0.5

d,PP

0.6

0.7

0.8

[-]

(b) RMSE of mass [-] (Penne) (b) RMSEcenter center of mass [-] (Penne)

0.7

0.7

0.6

0.4

0.8

0.35

0.6

0.1

0

0.1

0.2

0.3

0.4

0.5

d,PP

0.7

Fig. 12. Comparison between snapshots from an experimental video (left) and a video obtained with the DEM model (right) for Jelly Beans calibrated to (c) portion tail and front.

0.5

0.2

0

0.7

0.3

0.25

d,PB

[-]

0.5

0.6

0.8

add a bias to the calibration. In order to further investigate if the identified optima are physically sound, the predictions of the solver at the optimized parameter set → x p, opt were compared to the experimental references. Firstly, we assess the simulation’s capability to reproduce reference from the drop test used for calibration. As an example for visual comparison, screenshots of the experiment and simulation with the model calibrated to (c) are given in Figs. 12 and 13. Time stamps t1 to t3 are identical to the ones in Fig. 7. Time stamp t4 was chosen, so that the tail of the portion of Penne is visible. Since both experiment and simulation are subject to randomness, the location of individual particles should

(c) 1 RMSE portion frontfront [-] (Penne) (c) RMSE portion [-] (Penne) 1

0.02 5

0.6

5 0.0

[-]

0.5 d,PB

0.7

[-]

0.4

0.15

0.1

0.3 0.2

0.2 0.1 5

0.

5 0.0

0.1

1

0.1

0.2

0.3

0.4

0.5

d,PP

0.8

0.6

0.7

0.8

[-]

(c) 2 RMSE portion tail [-] (c) RMSE portion tail(Penne) [-] (Penne) 2

0.7

d,PB

[-]

0.6

0.7

0.5 1

0.4

0.5

0.3

0.8

1

1.4

0.2 0.1

0.2

0.7

0.1

0.2

0.3

0.5

0.4 d,PP

0.8

0.5

1

0.6

0.7

[-]

Fig. 11. 2D-projection of meta models for Penne.

reasons: Firstly, it is not possible to visualize the entire topology of the meta models in 2D. In fact, the effect of only two input parameters is visualized, while the remaining five input parameters are held constant. Secondly, even a perfect mathematical match is not necessarily a good physical match due to possible systematic errors in measurement, that

Fig. 13. Comparison between snapshots from an experimental video (left) and a video obtained with the DEM model (right) for penne calibrated to (c) portion tail and front. 8

Minerals Engineering 145 (2020) 106094

S. Kirsch

not be considered. For Jelly Beans, we find a good reproduction of the experimental reference regarding the location and local density of the bulk. For Penne, we find a slightly less faithful reproduction of the experiment. The overlaid curves of reference used for calibration and the respective simulations are presented in Figs. 14 and 15. We see a good reproduction of the experimental curves for the Jelly Beans in cases (a) and (c). In cases (b), both for Jelly Beans and for Penne the predictions become notably less accurate after the prediction cut-off. This could be caused by the fact that the time range after the cut-off exhibited large variation in both the experiment and the simulation, so a more accurate prediction could require a larger sample size. For Penne, we note a considerable deviation between reference and simulation in the model calibrated to the visible area (a). This could be the result of a faulty parameter set due to a non-unique optimum. If we re-assess the meta model projection of case (a) in Fig. 11, we find that the predicted error is generally larger than in the other cases. This means that the solver was not capable of providing a good match with the experiment, even in the ”best” regions. This could be an indication for undetected systematic errors in the measurements. Secondly, we assess the calibrated model’s capability to reproduce the standard funnel discharge test. The comparison of the discharge curves is presented in Fig. 16. In contrast to the results of the drop test in Fig. 14, the prediction fidelity is much less consistent. In particular, we find that the DEM model calibrated to (a) provides a very good match with the experimental curve, while overestimating the amount of bulk good remaining in the funnel at the end of discharge. In case (b), the DEM model over-predicts the mass flow rate. This leads to a faster and more complete discharge than found in the experiment. Finally, in case (c), the DEM model significantly under-predicts the mass flow rate, resulting in a longer discharge duration. In this case, we also find the largest deviation between the real and simulated discharge curves. A possible reason for this inconsistency depending on the calibration target is that the curves used for calibration exhibit a different sensitivity profile to the input contact parameters (see Table 5). While |PCC| between funnel discharge and μd, PP is 0.90, for (c)1 and (c)2, it is only 0.34 and 0.38 respectively. |PCC| between μd, PP and (a) and (b) are almost at the same level as for the funnel discharge. This means, calibration to (c), gives a less exact prediction for μd, PP than (a) and (b). Since the funnel discharge is, however, mostly dominated by μd, PP , it seems plausible that the DEM model calibrated to (c), therefore, gives a less accurate prediction here. However, μd, PP is also the second most important parameter in determining R (tres ) , so if calibration to (c) does indeed yield a faulty estimate for it, it could also affect the prediction of R (tres ) . As a last step, we assess the fidelity of the calibrated models regarding the previously mentioned parameter that is the most relevant for industrial filling: the residence time range R (tres ) . The box plots comparing the reference and the simulation are presented in Figs. 17 and 18. We find that the prediction for R (tres ) is much better for Jelly Beans than for Penne, with deviations of the mean of (a) 7.2%, (b) 0.4% and (c) 0.6% respectively. The deviations for Penne are (a) 9.5%, (b) 22.8% and (c) 2.4% respectively. A straightforward reason for this could be that the sphero-cylinder shape representation idealized the actual shape of Penne much more than that of the Jelly Beans. However, we do find a much more accurate prediction for R (tres ) of Penne for model (c). You could interpret the difference in prediction quality as an indication that the visible particle surface and the center of mass are not suitable calibration targets for Penne. This, again, can have two possible reasons: Firstly, several ambiguous optima that are not visible in the given 2D-representation of the meta models in Fig. 11 could be present. Thus, we might have found one that reproduced the curves ”by coincidence”, but is not actually physically sound. Secondly, a systematic error in measurement of surface and center of mass could bias the calibration. This would mean that even if a perfectly physically sound parameter set → x p were used, the simulation would not accurately

reproduce the curves. This seems especially plausible since the real particles are hollow and can thus be seen through in the video measurement (depending on orientation), whereas the virtual particles are solid. When assessing case (c), where the simulations were calibrated to portion fronts and tails, we find a very good prediction of the average residence time range R (tres ) for both Jelly Beans and Penne. You could infer that (c) is the best calibration criterion for the examined cases. With the opposite logic as before, this would mean that the optima in (c) are unambiguous and that there were no significant systematic measurement errors. It should be noted that criterion (c) is the most closely related to the residence time range R (tres ) . In the extreme case, when calibration and validation criterion are identical, a mathematical optimum for that criterion could have no physical meaning. This is especially true if several optima exist or if unknown systematic errors

Fig. 14. Comparison between experimental and simulated curves of the drop test with Jelly Beans. 9

Minerals Engineering 145 (2020) 106094

S. Kirsch

Fig. 16. Comparison between experimental and simulated curves of the funnel discharge test with Jelly Beans.

Fig. 15. Comparison between experimental and simulated curves of the drop test with Penne.

are present in the measurement of the reference or in the solver itself. This issue was discussed in detail e.g. by Pernot and Cailliez (2017). In our case, an unphysical choice of parameters could compensate systematic errors to produce a good match optimum for both criterion (c) and R (tres ) . The solution would then not be transferable, i.e. predicting the bulk good’s behavior in a different scenario would not be possible (Coetzee, 2017). However, considering the entire curves rather than only a single value on them adds a substantial amount of information to the optimization criterion. With this added information, we reasoned that a coincidental match is sufficiently unlikely. In summary, criterion (c) seems be the strongest, since it allowed not only a good reproduction of the experimental time signals but also the prediction of a separate ”arbitrary” criterion. This is despite the fact that, when model (c) for Jelly Beans was validated in funnel discharge, we saw some indication that the second most important parameter for

Fig. 17. Comparison between ranges of residence time R (tres ) of experiment and simulation for Jelly Beans.

R (tres ) , i.e. μd, PP , might have not been determined correctly. In fact, for the bulk good whose particles were represented more accurately, criteria (a) delivered a near equally good result as (c) in predicting R (tres ) . There is thus some indication, that a more accurate shape representation could make the calibration more robust for the use of other calibration criteria.

10

Minerals Engineering 145 (2020) 106094

S. Kirsch

found for the product with the less faithful DEM shape representation. The positions of the bulk portion front and tail over time seemed to be the best criteria for model calibration. Using them in a sum criterion lead to accurate reproduction of the references, even when the shape representation was less accurate. The findings could be relevant for industrial applications, where material parameters are often unknown and time dependent process data is available to calibrate DEM models without the need of a separate test. Declaration of competing interest The research and the publication of this work was performed while the author was employed at Robert Bosch Packaging Technology B.V.. All funding, materials and equipment used was also made available by Robert Bosch Packaging Technology B.V.. Apart from the approval on the study design chosen by the author, and providing personnel to assist in data collection, they played no further role in the analysis and interpretation of the presented data.

Fig. 18. Comparison between ranges of residence time R (tres ) of experiment and simulation for Penne.

Acknowledgments

4.3. Outlook

The author would like to thank Robert Bosch Packaging Technology B.V. for funding the research project. Further thanks go to their department for engineering and development for their constant support and helpful feedback. Lastly, the author would like to thank the employees of TU Dresden’s chair of Processing Machines and Processing Technology for their helpful support. Special thanks goes to the head of the chair for the academic supervision of the project.

The presented method has delivered promising results for industrial application and can be theoretically employed independently of the bulk good relevant for the given process. However, the calibration procedure requires a large amount of solver calls, which is partly due to the high dimensionality p = 7 of the problem. In some cases, numerical cost can increase further if either the individual solver calls become more expensive or more runs are required. The former could occur due to a bulk good with a large particle size distribution, i.e. smaller and thus more particles being present. The same is true, if complex particle shapes are required, as they are typically represented by assemblies of spherical sub-particles (Coetzee, 2016), increasing the number of virtual particles. Secondly, the dimensionality of the calibration problem can increase further, if a (wet) cohesive ore material is to be modeled. In this case, a cohesive model (Tomas, 2007) might be necessary to describe all relevant effects. This introduces at least one more parameter for cohesion to be calibrated (Roessler and Katterfeld, 2019). However, the necessary number of samples to maintain sampling density is proportional to the square root of dimensions p (Altman and Krzywinski, 2008). Thus, more samples, i.e. solver runs, then become necessary. In both cases, if computational cost becomes impractical, the dimensionality p can be reduced by eliminating from calibration parameters that play a minor role. This can be done based on a sensitivity study, as it was conducted for Jelly Beans. The eliminated parameters could then be estimated instead.

References Altman, N., Krzywinski, M., 2008. The curse(s) of dimensionality. Nat. Methods 399–400. https://doi.org/10.1038/s41592-018-0019-x. Barrios, G.K., de Carvalho, R.M., Kwade, A., Tavares, L.M., 2013. Contact parameter estimation for DEM simulation of iron ore pellet handling. Powder Technol. 248, 84–93. https://doi.org/10.1016/j.powtec.2013.01.063. discrete Element Modelling. . Barton, R.R., Meckesheimer, M., 2006. Chapter 18 metamodel-based simulation optimization. In: Henderson, S.G., Nelson, B.L. (Eds.), Handbooks in Operations Research and Management Science. vol. 13. Elsevier, pp. 535–574. https://doi.org/10.1016/ S0927-0507(06)13018-2. http://www.sciencedirect.com/science/article/pii/ S0927050706130182. Cleary, P.W., Morrison, R.D., 2011. Understanding fine ore breakage in a laboratory scale ball mill using DEM. Miner. Eng. 24 (3), 352–366. https://doi.org/10.1016/j.mineng. 2010.12.013. special issue: Comminution. . Cleary, P.W., Sawley, M.L., 2002. DEM modelling of industrial granular flows: 3D case studies and the effect of particle shape on hopper discharge. Appl. Math. Model. 26 (2), 89–111. https://doi.org/10.1016/S0307-904X(01)00050-6. . Coetzee, C.J., 2016. Calibration of the discrete element method and the effect of particle shape. Powder Technol. 297, 50–70. https://doi.org/10.1016/j.powtec.2016.04.003. Coetzee, C.J., 2017. Review: Calibration of the discrete element method. Powder Technol. 310, 104–142. Coetzee, C.J., Nel, R., 2014. Calibration of discrete element properties and the modelling of packed rock beds. Powder Technol. 264, 332–342. https://doi.org/10.1016/j. powtec.2014.05.063. Djordjevic, N., Shi, F., Morrison, R., 2003. Applying discrete element modelling to vertical and horizontal shaft impact crushers. Miner. Eng. 16 (10), 983–991. https://doi.org/ 10.1016/j.mineng.2003.08.007. . Dynardo GmbH, 2019. Methods for multi-disciplinary optimization and robustness analysis. ESSS, 2019. Rocky DEM technical manual. Frank, M., Holzweißig, J., 2016. Simulation-based optimization of geometry and motion of a vertical tubular bag machine. Sächsische Landesbibliothek. . González-Montellano, C., Fuentes, J., Ayuga-Téllez, E., Ayuga, F., 2012. Determination of the mechanical properties of maize grains and olives required for use in DEM simulations. J. Food Eng. 111 (4), 553–562. https://doi.org/10.1016/j.jfoodeng.2012. 03.017. http://www.sciencedirect.com/science/article/pii/S0260877412001598?_ rdoc=1&_fmt=high&_origin=gateway&_docanchor=&md5= b8429449ccfc9c30159a5f9aeaa92ffb. Gröger, T., Katterfeld, A., 2006. On the numerical calibration of discrete element models for the simulation of bulk solids. In: 16th European Symposium on Computer Aided Process Engineering and 9th International Symposium on Process Systems

5. Conclusion We set out to investigate whether DEM model parameter identification can be carried out effectively in-situ when time signals are used as calibration references. Using time signals over a single scalar reference was assumed to reduce ambiguity with a higher chance of a physically sound parameter set. A total of six model calibrations were performed, using two different bulk materials and applying three different reference criteria for each. The calibration procedure with the drop test was implemented with an adaptive Kriging regression approach, capable of handling randomness in the noisy solvers. The optimization of model parameters was performed on the resulting meta models. We found indication that the chosen approach avoids ambiguous solutions. The DEM models were validated in the drop process and in a standard calibration test. The calibration targets’ sensitivity to input parameters seems to be relevant for model transferability. The models calibrated to target parameters, with a sensitivity profile similar to the validation process, showed the best match with the experiment. Also, some notable deviation between reference and simulations was 11

Minerals Engineering 145 (2020) 106094

S. Kirsch

Methods (DEM)’07. . Most, T., 2019. Personal communication on adpative meta modeling in Optislang. Most, T., Will, J., 2008. Metamodel of optimal prognosis – an automatic approach for variable reduction and optimal metamodel selection. Weimar Optimization and Stochastic Days. http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid= 56BFA6867DBA96D6C200BCACD14EA530?doi=10.1.1.411.1767&rep=rep1& type=pdf. Most, T., Will, J., 2011. Sensitivity analysis using the metamodel of optimal prognosis. Weimar Optimization and Stochastic Days. . Paulick, M., Morgeneyer, M., Kwade, A., 2015. Review on the influence of elastic particle properties on DEM simulation results. Powder Technol. 283, 66–76. https://doi.org/ 10.1016/j.powtec.2015.03.040. Pernot, P., Cailliez, F., 2017. A critical review of statistical calibration/prediction models handling data inconsistency and model inadequacy. AIChE J. 63 (10), 4642–4665. https://doi.org/10.1002/aic.15781. . Rackl, M., Grötsch, F.E., 2018. 3d scans, angles of repose and bulk densities of 108 bulk material heaps. Scientific Data 5https://doi.org/10.1038/data.2018.102. . Rackl, M., Hanley, K.J., 2017. A methodical calibration procedure for discrete element models. Powder Technol. 307, 73–83. https://doi.org/10.1016/j.powtec.2016.11. 048. . Roessler, T., Katterfeld, A., 2019. DEM parameter calibration of cohesive bulk materials using a simple angle of repose test. Particuology 45, 105–115. https://doi.org/10. 1016/j.partic.2018.08.005. https://www.sciencedirect.com/science/article/pii/ S1674200119300033. Roessler, T., Richter, C., Katterfeld, A., Will, F., 2019. Development of a standard calibration procedure for the DEM parameters of cohesionless bulk materials – part I: Solving the problem of ambiguous parameter combinations. Powder Technol. 343, 803–812. https://doi.org/10.1016/j.powtec.2018.11.034. . Rogers, J.L., Nicewander, W.A., 1988. Thirteen ways to look at the correlation coefficient. Am. Stat. 42 (1), 59–66. https://doi.org/10.2307/2685263. . Schneider, B.J., 2012. Polygonale diskrete Elemente zur Modellierung heterogener Materialien (Phd thesis). Universität Stuttgart. Tomas, J., 2007. Adhesion of ultrafine particles–a micromechanical approach. Chem. Eng. Sci. 62 (7), 1997–2010. https://doi.org/10.1016/j.ces.2006.12.055. . Walton, O.R., Braun, R.L., 1986. Viscosity, granular temperature, and stress calculations for shearing assemblies of inelastic, frictional disks. J. Rheol. https://doi.org/10. 1122/1.549893. http://eds.a.ebscohost.com/eds/detail/detail?vid=2&sid= ac521e56-c49e-476d-bddf-178c01016107. Wensrich, C., Katterfeld, A., 2012. Rolling friction as a technique for modelling particle shape in DEM. Powder Technol. 217, 409–417. https://doi.org/10.1016/j.powtec. 2011.10.057. . Zhu, H., Zhou, Z., Yang, R., Yu, A., 2008. Discrete particle simulation of particulate systems: a review of major applications and findings. Chem. Eng. Sci. 63 (23), 5728–5770. https://doi.org/10.1016/j.ces.2008.08.006. .

Engineering. Elsevier BV, pp. 533–538. https://doi.org/10.1016/s1570-7946(06) 80100-8. Grima, A., Fraser, T., Hastie, D., Wypych, P., 2011. Discrete element modelling: Troubleshooting and optimisation tool for chute design. In: BELTCON 16th International Materials Handling Conference. volume 16. Grima, A., Wypych, P., 2010. Discrete element simulation of a conveyor impact-plate transfer: Calibration, validation and scale-up. Aust. Bulk Handling Rev. 64–72. Harzanagh, A.A., Orhan, E.C., Ergun, S.L., 2018. Discrete element modelling of vibrating screens. Miner. Eng. 121, 107–121. https://doi.org/10.1016/j.mineng.2018.03.010. . Helton, J., Davis, F., 2003. Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems. Reliab. Eng. Syst. Saf. 81, 23–69. . Jones, D.R., Schonlau, M., Welch, W.J., 1998. Efficient global optimization of expensive black-box functions. J. Global Optim. 13 (4), 455–492. https://doi.org/10.1023/ A:1008306431147. Katterfeld, A., Coetzee, C.J., Donohue, T., Fottner, J., Grima, A., Álvaro Ramírez-Gómez,, Ilic, D., Kačianauskas, R., Necas, J., Schott, D., Williams, K., Zegzulka, J., 2019. Calibration of DEM parameters for cohesionless bulk materials under rapid flow conditions and low consolidation. White Paper. https://doi.org/10.13140/RG.2.2. 26318.31048/1. Kirsch, S. 2018. DEM model calibration for vertical filling: Selection of adequate trials and handling randomness. In: Schwarz, H. (ed.), 15th Weimar Optimization and Stochastic Days. DYNARDO (Dynamic Software and Engineering) GmbH. . Kirsch, S., Philipp, A., 2018. Simulation of vertical filling processes of granular foods for typical retail amounts. In: 9th Conference Processing Machines and Packaging Technology Dresden. . Kleijnen, J., 2017. Simulation optimization through regression or kriging metamodels. Discussion Paper. . Liu, F., Zhang, J., Li, B., Chen, J., 2016. Calibration of parameters of wheat required in discrete element method simulation based on repose angle of particle heap. Nongye Gongcheng Xuebao/Transactions of the Chinese Society of Agricultural Engineering. Liu, H., Ong, Y.S., Cai, J., 2018. A survey of adaptive sampling for global metamodeling in support of simulation-based complex engineering design. Struct. Multidiscip. Optim. 57 (1), 393–416. https://doi.org/10.1007/s00158-017-1739-8. Luding, S., 1998. Collisions & contacts between two particles. Phys. Dry Granular Media 350, 285–304. . Marigo, M., Stitt, E.H., 2015. Discrete element method (DEM) for industrial applications: comments on calibration and validation for the modelling of cylindrical pellets. KONA Powder Particle J. 32, 236–252. https://doi.org/10.14356/kona.2015016. . Markauskas, D., Kačianauskas, R., 2010. Investigation of rice grain flow by multi-sphere particle model with rolling resistance. Granular Matter 13 (2), 143–148. https://doi. org/10.1007/s10035-010-0196-5. McKay, M.D., Beckman, R.J., Conover, W.J., 1979. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (2), 239–245. https://doi.org/10.2307/1268522. Montgomery, D.C., Runger, G.C., 2003. Applied Statistics and Probability for Engineers. John Wiley & Sons Inc. Morrison, R.D., Cleary, P.W., 2008. Towards a virtual comminution machine. Miner. Eng. 21 (11), 770–781. https://doi.org/10.1016/j.mineng.2008.06.005. discrete Element

12