Replicating cohesive and stress-history-dependent behavior of bulk solids Feasibility and definiteness in DEM calibration procedure

This paper presents a multi-step DEM calibration procedure for cohesive solid materials, incorporating feasibility in ﬁnding a non-empty solution space and deﬁniteness in capturing bulk responses independently of calibration targets. Our procedure follows four steps: (I) feasibility; (II) screening of DEM variables; (III) surrogate modeling-based optimization; and (IV) veriﬁcation. Both types of input parameter, continuous (e.g. coefﬁcient of static friction) and categorical (e.g. contact module), can be used in our calibration procedure. The cohesive and stress-history-dependent behavior of a moist iron ore sample is replicated using experimental data from four different laboratory tests, such as a ring shear test. This results in a high number of bulk responses (i.e. (cid:1) 4) as calibration targets in combination with a high number of signiﬁcant DEM input variables (i.e. > 2) in the calibration procedure. Coefﬁcient of static friction, surface energy, and particle shear modulus are found to be the most signiﬁcant continuous variables for the simulated processes. The optimal DEM parameter set and its deﬁniteness are veriﬁed using 20 different bulk response values. The multi-step optimization framework thus can be used to calibrate material models when both a high number of input variables (i.e. > 2) and a high number of calibration targets (i.e. (cid:1) 4) are involved.


Introduction
To simulate, design, and optimize processes and equipment for handling bulk solids, such as iron ore and coal, the discrete element method (DEM) is a suitable computational method.However, DEM simulations can only predict bulk level responses (e.g.shear strength) accurately if their input parameters are selected appropriately.To select the input parameters with confidence, the common procedure is to calibrate and to validate DEM simulations [1][2][3][4].The calibration can be done by finding an optimal combination set of DEM input parameters that replicates the captured bulk response [5].
Over the past decade, reliable DEM calibration procedures have been developed to model free-flowing bulk solids, such as iron ore pellets [1], glass beads [6], sinter ore [7], sand [8,9], and gravel [10,11].By setting multiple targets for the DEM calibration, more than a single bulk response can be considered.This prevents the ''ambiguous parameter combinations" problem in the DEM calibration procedure, which is discussed in detail in [11].For example, to calibrate DEM input variables for simulating iron pellets in interaction with ship unloader grabs, Lommen et al. [1] considered at least three different calibration targets.They replicated the static angle of repose using the ledge and free-cone methods; the penetration resistance of iron pellets was also replicated, using a wedge penetration test setup.
In general, DEM calibration is performed following the generic procedure shown in Fig. 1.To find an optimal combination of DEM input parameters that satisfies multiple calibration targets, optimization methods can offer a solution.Various optimization methods have already been applied to calibrate the continuous type of DEM variables successfully [6,7,10,12].Continuous DEM variables are numerical variables that have an infinite number of values between any two values [13].For example, the coefficient of static friction is an important continuous DEM variable during calibration [5].Richter et al. [10] concluded that surrogate modeling-based optimization methods are most promising for DEM calibration when continuous variables are included.
Categorical-type DEM variables have not yet been included in the calibration procedure when optimization methods are used.Categorical variables are finite numbers of groups or categories that might not have a logical order [13].For example, shape of particles is a DEM categorical variable that plays an important role during calibration [14].One can use design of experiments (DoE) methods to include categorical variables in the DEM calibration procedure.However, a high number of simulations might have to be run with no guarantee of finding an optimal set of DEM input parameters [10].Additionally, iron ore fines and other similar bulk solids (e.g.coal) have an irregular distribution of particle shape [15] as well as fine particle sizes [16].Modeling accurate particle shapes and sizes for cohesive bulk solids in DEM simulations thus leads to a computation time that is generally impractical for studying industrial bulk handling processes, such as flow in silo [16].
Furthermore, selecting an appropriate contact model from the available options is an important challenge in the DEM calibration.Applying optimization methods without choosing a proper contact model might, for example, lead to an empty solution space or inadequacy in meeting macroscopic bulk behaviors other than the selected calibration targets [7].A contact model generally includes multiple modules to calculate forces and torques between elements (e.g.particles).Fig. 2 schematically illustrates a contact spring-damper system between two particles, a and b.Here, three main modules are identified: contact force in the normal direction is denoted by f N , while f T and s R represent force in the tangential direction and rotational torque respectively.Contact modules can be selected independently of each other.For instance, a rolling friction module can be implemented in various ways to determine rotational torque between two particles [1,17,18].Therefore, each module of the contact model can be considered as a categorical variable in DEM calibration.
By contrast with free-flowing materials, cohesive bulk solids such as moist iron ore fines usually show a stress-historydependent and cohesive behavior [19].Their bulk responses, such as shear strength, bulk density, and penetration resistance, depend on the history of applied normal pressure on the bulk specimen [19][20][21].Thiss stress-history-dependent and cohesive behavior can be simulated by using contact models based on an elastoplastic adhesive spring [20,[22][23][24][25]. Orefice and Khinast [25] used a multi-stage sequential DEM calibration procedure to model cohesive bulk solids using a linear elasto-plastic adhesive model; the calibration was done by replicating a specific bulk response at each stage, starting with the angle of repose (measured using the funnel test) as the first calibration target.Three continuous DEM variables were included during the calibration; other DEM input parameters, continuous and categorical, needed to be kept constant during their calibration procedure.The multi-stage sequential calibration procedure might fail to meet the following criteria.
-Feasibility.Replicating all the selected bulk responses can be infeasible using chosen values for the input parameters that are constant during the calibration, such as a specific contact module.Therefore, considering the necessity of including multiple calibration targets, the calibration procedure can lead to an empty solution space for one or more than one of the calibration targets.-Definiteness or avoiding ''ambiguous parameter combinations" [11].To meet this criterion, a bulk response independent of the calibration targets needs to be simulated successfully using the calibrated set of DEM input parameters.Additionally, properly selecting all modules of the DEM contact model is a prerequisite.Otherwise, the calibrated set of input parameters might fail to capture a bulk response different than the selected calibration targets.
For example, the ''definiteness" criterion has been focused on in the automated calibration procedure developed by [26], which is based on a genetic algorithm to replicate stress-historydependent and cohesive behavior of bulk solids in the ring shear test.By introducing cohesive forces as well as elasto-plastic stiffness into the DEM calibration procedure, the number of DEM input variables and the number of required bulk responses increase [25,[27][28][29].For that reason, the abovementioned criteria become important in developing a reliable calibration procedure to simulate cohesive and stress-history-dependent behavior of bulk solids.As yet, however, no literature has addressed how to ensure that both criteria, feasibility and definiteness, are met in a DEM calibration procedure considering both continuous and categorical DEM input variables.Additionally, calibrating DEM input parameters is still a challenge when a high number (i.e.> 2) of variables in combination with a high number of bulk responses (i.e.> 2) is involved.
In this paper, we develop a reliable multi-step DEM calibration procedure to capture the cohesive and stress-history-dependent behavior of bulk solids.In each step of the calibration procedure, the variables space is narrowed down to be further optimized in the next step.The first step uses a feasibility analysis, based on Latin hypercube design (LHD), to choose a suitable contact model by efficiently searching for a non-empty solution space.This ensures that the calibration procedure meets the ''feasibility" criterion.The second step screens the significant DEM variables to quantify their influence on the selected bulk responses.This allows us to find an optimal combination of DEM input variables in the  Fig. 2. A contact spring-damper system between two particles, including normal, tangential, and rotational directions.
third step by applying a surrogate modeling-based optimization method.In the third step, we use a different set of calibration targets, compared to the first and second steps, to consider the ''definiteness" criterion.The final step is to verify the adequacy of the optimal combination in replicating the cohesive and stresshistory-dependent behavior for several bulk responses, such as bulk density, shear strength, and penetration resistance.

DEM calibration procedure: a multi-step optimization framework
In general, a calibration procedure aims at identifying an optimal combination of DEM input parameters, Ny , adequately similar to responses captured in physical laboratory or in-situ tests, Y ¼ y 1 ; Á Á Á ; y Ny [5].Ns is the number of DEM input parameters and Ny the number of calibration targets.Bulk responses such as bulk density and shear strength thus need to be determined first, using appropriate physical tests.This allows us to set calibration targets and to quantify the difference in bulk responses between simulated and physically determined values.To ensure that feasibility and definiteness criteria are satisfied for multiple calibration targets, a multi-step DEM calibration procedure considering categorical input parameters is proposed in Fig. 3.The following four steps are included: (I) feasibility; (II) screening of DEM variables; (III) surrogate modeling-based optimization; and (IV) verification.
To apply surrogate modeling-based optimization, the parameter space needs to be searched effectively to be able to approximate Y 0 .Accordingly, F(X) maps relationships between new calibration targets, Y = y 1 , . .., y My , and (significant) DEM variables.
Although the full factorial design can be used to create multivariate samples, all the possible combinations between significant DEM variables must be included.This leads to a high number of simulations needing to be done.Fractional factorial designs, such as Taguchi [30], Placket-Burmann [31], and Box Behnken [32] designs, can be used to generate multi-variate samples required for surrogate modeling without the need to create all the possible combinations of variables.For example, if a full factorial design is used for 4 input variables having 3 levels each, that leads to 3 4 = 81 combinations to run.Using the Taguchi (orthogonal) method, a fractional factorial design can be created by running only 9 or 27 possible combinations.The accuracy of the surrogate model is evaluated using the coefficient of determination, R 2 .This coefficient quantifies the surrogate model accuracy in representing variability of values obtained from DEM simulations.To ensure that the surrogate model converges to a verifiable X*, a minimum R 2 value of 0.75 is considered to be met for all calibration targets.Otherwise, more samples are used to train the surrogate model.
Next, the response optimizer searches for an optimal combination of input variables, X*, that jointly meets a set of calibration targets, Y.To find X* using the surrogate model, we use the response optimizer toolbox available in Minitab [33].
The mean of absolute relative differences is used to quantify error in the verification step.If y and y' represent measured bulk responses in the experiment and the simulation, respectively, then |e| mean is determined according to Eq. ( 1) for a number of bulk responses, N e .In the current study, an |e| mean 10% is considered an acceptable outcome during verification.
Therefore, in each step of the calibration procedure the variables space is narrowed down to be further optimized in the next step.In the final step, a verified parameter set is found by checking |e| mean .

DEM calibration targets: Y
In this study, DEM calibration targets are set to values reported from our comprehensive measurements campaign on cohesive iron ores [19].Bulk property variability of cohesive iron ores has been characterized using the following laboratory tests:  A, B, and C) are used in the current study to set DEM calibration targets.During the calibration procedure, two out of three influencing parameters, MC and r, are considered as sources of possible bulk property variability.Below we describe characteristics of the selected bulk solid sample as well the measured bulk responses.

Bulk solid sample
The bulk solid sample is a sinter feed type of iron ore from the Carajas mines, one of the largest iron ore resources on earth [34].The average density of the particles is 4500 kg/m 3  , with a standard deviation of 125 kg/m 3 .The median particle size, d 50 , is equal to 0.88 mm [35].The dry-based moisture content was determined according to the method described in [36], in which the sample is dried using a ventilated oven.This resulted in MC = 8.7%.An overview of measured properties of the sample is presented in Table 1.

Measured bulk responses
Table 2 displays physically measured bulk responses of the sample using the ring shear and ledge angle of repose tests when r pre 20 kPa and D MC = ±2%.Pre-consolidation or pre-shear stress, r pre , is a normal confining pressure that is applied initially.In the ring shear test, for example, a normal confining pressure of 20 kPa is applied initially during the pre-shear stage, and next a normal confining pressure of 2 kPa (r shear ) is applied.Fig. 4 shows the results of shear stress measurements in the ring shear test, including one pre-shear stage and one shearing stage.In general, r shear is smaller than r pre , which allows us to investigate a stresshistory-dependent bulk response, such as shear strength in the case of shear tests.The ledge angle of repose test has been conducted under no pre-consolidation stress, which represents the free-surface flow of bulk solids under gravity force.Maximum and minimum values of physically measured bulk responses are shown under D MC, up to ± 2%, compared to its as-received condition.By considering the maximum and minimum measured values of bulk responses, extreme values can be included in the feasibility evaluation step of the DEM calibration procedure.In other words, the feasibility is evaluated for a range of bulk response values.
According to the Mohr-Coulomb equation, the shear strength of bulk material s s is often approximated by Eq. ( 2) [37]: where tan(u) indicates the angle of internal friction.c is the shear strength of the bulk material when r n ¼ 0, thus it denotes the cohesion of the bulk material.Eq. ( 2) suggests that increasing the normal stress r n, decreases the contribution of c to the shear strength.Additionally, increasing r n results in a higher contribution of particle-particle friction to shear strength.
The wall friction was also determined in [19]; this was done using the ring shear test by applying small adjustments according Step I. Feasibility Step II.Significant DEM variables Step III.Surrogate modelling-based optimization Step IV.Verification • Searching for a feasible solution space that covers selected bulk responses.
• Using DoE techniques, variables space can be searched effectively using a minimum number of sampling points.
• X is feasible if a satisfactory coverage of solution space is reached.
• A sensitivity analysis to identify significance of DEM variables.
• One-variable-at-a-time (OVAT) is the most suitable DoE technique for this step.
• To ensure that the definiteness criterion is met, a different set of calibration targets is used in this step, compared to previous steps.
• F(X) maps relationships between new calibration targets and (significant) DEM variables, X.

Verified X *
• |e| mean is used to quantify error.
• The definiteness of X * is confirmed if bulk response(s) different than the calibration targets are simulated successfully.
Definiteness: Verify X * for various bulk responses (i.e.|e| mean ≤ 10%)  to [38].The test was done with a r pre equal to 20 kPa and then the wall friction was measured for eight different levels of r shear between 2 and 17 kPa.The wall friction measurements resulted in a wall yield locus with an average wall friction angle of 19°a nd a negligible adhesion strength of 0.1 kPa.Table 3 displays measured bulk responses of the sample using the consolidation-penetration test when r pre !65 kPa and D MC = 0%.This test procedure is designed to represent the penetration resistance of iron ore cargoes during ship unloading when grabs are being used [39].To consider the stress-history dependency, two levels of r pre are included in the calibration procedure, equal to 65 and 300 kPa, respectively.As the first bulk response parameter, accumulative penetration resistance [J] on the wedgeshaped penetration tool is determined by integrating the reaction force over penetration depth [40].The secondary measured bulk response in the test is the bulk density after removing r pre .For example, after removing r pre of 300 kPa, the bulk density was measured according to the procedure described in [39], which for this sample was equal to 2807 kg/m 3 on average for three test iterations.Therefore, bulk property variability of the cohesive iron ore sample has been determined under variation of confining pressure as well as moisture content.This provides a comprehensive set of measurement data to be used in the DEM calibration procedure (illustrated in Fig. 3).

Contact modules in normal and tangential directions: elastoplastic adhesive
The EDEM (v2020) software package is used to create and run simulations.To capture the stress-history-dependent bulk responses as well as cohesive forces, an elasto-plastic adhesive contact model built into the software package is used.This is formulated in [41] under the name Edinburgh Elasto-Plastic Adhesive (EEPA).For details, refer to [41] or [20,26,42].This model has been used successfully by [20,26] to simulate bulk responses of cohesive bulk solids.
A schematic diagram of the EEPA contact spring for the normal direction of f-d (force-overlap) is provided in Fig. 5.The contact spring consists of three main spring-based parts during loading and unloading, as well as a constant pull-off force, f 0 .
Part 1.The contact starts with the loading part, with spring stiffness of k 1 , when the distance between the centers of two approaching particles is smaller the sum of their radiuses.The non-linear mode of the contact module is used in the current study by setting the slope exponent value to 1.5.
Part 2. By reducing the contact force, unloading commences; during this process, the plastic deformation is replicated by switching the spring stiffness to k 2 .The plasticity ratio, k P , determines the ratio between k 2 and k 1 .
Part 3. As unloading continues, a minimum attractive (adhesive) force is reached that is denoted by f min .The limit is determined using Eq. ( 3) [41].
where Dc and a are surface adhesion energy [J/m 2 ] and contact radius [m], respectively.If the unloading of the contact spring continues, the f À d follows the adhesive path with stiffness of Àk adh .In this study, an exponent value equal to 1.5 is used for d in part 3, which is similar to the slope exponent value used in part 1.Therefore, during the calibration procedure the adhesion path (part 3) can be controlled by adjusting f 0 and Dc as DEM input variables.The tangential stiffness of the contact model is varied as a multiplier, k t,mult , of the initial loading stiffness.

Simulation setups
DEM simulation setups are created representing the physical laboratory tests in the geometry scale of 1:1.

(A) Ring shear test
The ring shear test device used in [19] to characterize the shear strength of the iron ore sample is the same as the device used in our earlier study [26].For that reason, the same simulation setup and procedure is applied in this study.Fig. 6a and b illustrate components of the ring shear test in laboratory and simulation   Accumulative penetration resistance at 70 mm depth when r pre = 300 kPa W 70,300 J 121 5 Bulk density after applying r pre = 65 kPa q b,65 kg/m 3 2668 65 Bulk density after applying r pre = 300 kPa q b,300 kg/m 3 2807 14 environments, respectively.In the simulation setup, we use cylindrical periodic boundaries to simulate a quarter of the shear cell (Fig. 6b).This allows us to reduce the computation time by 50% with no undesirable influence on the simulation accuracy [26].
(B) Ledge angle of repose test A ledge test method, according to [1], for measuring the static angle of repose, a M , of the cohesive iron ore sample was used in [19].The test setup and its procedure are also referred to by other names in literature, such as ''shear box" [44] and ''rectangular container test" [8].Fig. 7a and b show the test box dimensions, including the slope formed after failure, in laboratory and simulation environments, respectively.The container is 250 mm high, but it has been filled only to the flap opening's height at 200 mm.After opening the flap, the bulk solid can thus flow out from the container.Once a static angle of repose is created, a M is quantified by applying the linear regression technique to fit a line on the slope of bulk surface.
(C) Uni-axial consolidation-penetration test Fig. 8 shows three main components of the consolidationpenetration test: a container, a lid, and a wedge-shaped penetration tool.The lid's surface area is equal to the container's sectional area.The wedge-shaped tool is 200 mm long, which allows it to create a plane contact with particles.
The procedure of the simulated consolidation-penetration test is illustrated in Fig. 9.
-First, the container is filled with DEM particles.A stable situation is reached when the maximum velocity of the particles is smaller 0.1 mm/s.-Second, the lid is moved downward with a constant velocity of 10 mm/s to create a consolidated situation.This is continued until the desired pressure on the lid is reached (i.e.65 and 300 kPa).-Third, the sample is unloaded by moving the lid upward with a velocity of 10 mm/s.-Finally, the wedge-shaped tool is moved downward with a velocity of 10 mm/s, similar to the laboratory test procedure [19].

Initial sampling strategy for step I (Feasibility) using LHD
The initial sampling aims at evaluating the feasibility of capturing calibration targets using selected DEM input constants and variables.This allows us to select a suitable solution, including levels of categorical variables and constants.Two simulation setups, ring shear and ledge angle of repose tests, are used in step I, feasibility.This means that the shear flow in two different test setups is simulated for r pre of up to 20 kPa.Three different bulk responses, s pre=20 , s 2;20 , and a M (angle of repose), are analyzed using DEM simulations for various combinations of input parameters.
During a calibration procedure, DEM input parameters, X ¼ x 1 ; Á Á Á ; x Ns , are divided into two groups: input variables and constants.Level input variables are varied in a range to meet calibration targets (Fig. 1).Levels of DEM input constants are chosen based on available literature, if applicable; otherwise, their level is selected based on rational assumptions, as recommended by [25], or by the direct measurement method, as discussed in [5].For example, modeling the actual shape and size distribution of a cohesive iron ore sample leads to a computational time that is impractical [45,46].For that reason, a simplified representation of particle shape and size can be used to develop a DEM simulation of cohesive iron ore.This technique has been applied successfully by [20,26,47] to model bulk solids that have fine particles with irregular shape distribution.Nevertheless, the rotational torque between particles needs to be considered; according to [48], two options are possible: (a) introducing a certain level of non-sphericity in particle shape; and/or (b) suppressing the rotational freedom of particles.In this study, option (b) is applied, as -compared to using multispherical particles -it does not have a negative influence on the computational time.The rotational freedom of particles can be suppressed artificially by either introducing a rolling friction module [17] or restricting the rotation of the particles [1,26,49].Both techniques are included as a categorical variable in step I, feasibility.The rolling friction module is implemented according to [18].This implementation was classified as ''rolling model C" by [17], so we refer to the rolling friction module as RC in this article.Restricting the rotation of particles is done by applying a counterbalance torque in each time-step necessary to prevent rotational movement.This leads to an increase in the particles' resistance to rotational torque.Restricting the rotation of particles has been used successfully to resemble realistic material behavior [1,24,40,48].Additionally, the number of input variables is reduced because, when using the restricted rotation (RR) technique, rolling friction coefficient does not play a role in rotational torque.

DEM input variables when RC option is used
Table 4 displays DEM input variables when the RC option, rolling friction module C, is used.Based on the available literature, the coefficient of static friction between particles, m s,p-p , is probably the most influential parameter on the internal shear strength of bulk solids [7,14,28,[50][51][52][53][54][55][56][57][58][59][60].Coefficient of rolling friction is also usually considered as an influential variable on shear flow [5].To calibrate the shear flow of cohesive bulk solids, [61] found that a range of 0.2 to 1.0 is reasonable for coefficients of static friction and rolling friction when rolling model C is used.Particle shear modulus determines the stiffness of the contact spring.Therefore, G, particle shear modulus, is included as a continuous DEM variable in our investigation.A range between 2.5 and 10 MPa is used for G, which covers values used by other researchers modeling cohesive bulk solids using the same elasto-plastic contact model [20,26].
Constant pull-off force (f 0 ) and surface energy (Dc) are included in the calibration to control the magnitude of adhesive forces in the contact spring.f 0 is varied between À0.0005 and À0.005 N, and Dc between 5 and 50 J/m 2 .These ranges are expected to be sufficient to capture a realistic shear flow based on the DEM calibration done in [26].

DEM input variables when RR option is used
Table 5 displays DEM input variables when the RR option, rotation restricted, is used.First, based on our simulation results reported in [26], the ranges of coefficient of static friction and surface energy are changed, compared to the values in Table 4.By restricting the rotation of particles, their mobility decreases and so lower restrictive forces (e.g.cohesive and friction) can be used during the calibration procedure, compared to the case when the RC option is used.The coefficient of static friction is varied between 0.2 and 0.4, while the surface energy variation is between 2.5 and 25 J/m 2 .Second, ranges of other input variables are similar to the case when the RC option is used.

DEM input constants
Table 6 presents other DEM input parameters that are kept constant during initial sampling for step I, feasibility.Particle density is set to 4500 kg/m 3 , similar to the measured value (Table 1).As discussed earlier, the representation of particles' shape and size is simplified.Spherical particles are used and the mean particle diameter value is set to 4 mm including a normal particle size distribution with a standard deviation of 0.1.In addition to a reasonable computation time when spherical particles are used, the coarse graining principles for the elasto-plastic adhesive contact model [46] can be applied during the calibration procedure to further minimize the computation time.For example, the ledge angle  of repose simulations are done using coarse grained particles with a scaling factor of S P = 2.25, as per [46].Constant pull-off force and surface energy are scaled with factors of S P 2 and S P to maintain comparable bulk responses with the unscaled simulation.For further details of particle scaling rules, please refer to [46].
The tangential stiffness multiplier, k t,mult., is recommended as 2/3 [62] for non-linear elastic contact springs.According to [63], to maintain simultaneous harmonic oscillatory positions between normal and tangential elastic springs, a value of 2/7 is   recommended.However, no recommendation was found in literature to select k t,mult.when a non-linear elasto-plastic normal spring is used.For that reason, a range of k t,mult.bounded to 0.2 to 1 was used in the ledge angle of repose simulation.Within this range, no significant influence on the simulation stability and simulated bulk responses were found, and therefore k t,mult.is set to 0.4.As suggested by [26], if a negligible adhesion strength is measured in the wall friction test, the Hertz-Mindlin (no-slip) contact model [64] can be used to describe interaction between particles and geometry.The sliding friction coefficient between particles and wall geometry, m s,p-w , is therefore determined directly by Eq. ( 4), which results in m s,p-w = 0.37 for the measured average angle of the wall yield locus u x of 19 (Section 2.2).The rolling friction coefficient between particles and wall geometry has a negligible influence on simulated shear stress [65], and therefore m r,p-w is set to 0.5.

Initial samples
Using design of experiments (DoE) techniques, parameter spaces -including their levels and possible combinations -can be searched effectively using a minimum number of sampling points.A Latin hypercube design (LHD) is constructed in such a way that each of the parameters is divided into p equal levels, where p is the number of samples.Based on the U P criterion [66], the location of levels for each parameter is randomly, simultaneously, and evenly distributed over the parameter spaces, maintaining a maximized distance between each point.The LHD is constructed according to the algorithm developed in [67], which satisfies the U P criterion for up to 6 parameters.This allows us to include up to 6 DEM input parameters in a feasibility evaluation.Fig. 10 displays levels of the 5 continuous DEM input variables at S P = 1 when the RR option, restricted rotation, is used.Forty different samples are created using the LHD to simulate ring shear and ledge angle of repose tests.Similarly, using the LHD, 40 different samples are created for the 6 continuous DEM input variables (based on Table 4) at S P = 1 when the RC option, rolling friction module C, is used.
In total, 160 simulations are run during step I, feasibility, which cover 2 categorical variables and 6 continuous variables.

Results
In this section, first the simulation results of the initial samples (step I) are presented.Then a feasible solution is chosen to continue the calibration procedure when executing its next steps.Additionally, new samples are created at the beginning of each new step to meet its specific objective.Step I: feasibility Fig. 11 displays the simulation results of the 40 initial samples when the RC option, rolling friction module C, is used.Three different bulk responses are quantified: -shear stress in the pre-shear stage, s pre=20 ; -shear stress in the shearing stage, s 2:20 ; and, -average angle of repose in the ledge test, a M .Thus, N y is equal to 3 in step I, feasibility evaluation.Simulation results are also compared with the maximum and minimum values that were measured in the laboratory environment (shown in Table 2).For example, s exp.max and s exp.min are shown using blue and red dashed lines respectively.
Using the RC option, a range of s pre=20 bounded to 6.2 and 12.3 kPa is captured.This shows that the 40 samples created using LHD could vary s pre=20 by around 100%.The maximum simulated s pre=20, 12.3 kPa, is around 25% lower than s exp.min .This means that simulating a comparable s pre=20 is probably infeasible using the RC option.To confirm whether this conclusion is limited to the selected ranges of the 6 DEM input variables, additional simulations using extreme values of DEM input variables are conducted.Extreme values are selected outside the selected ranges shown in Table 4.For example, using sample 32, which produced s pre=20 = 12.3 kPa, an additional sample is created by increasing particle shear modulus, G, to 100 MPa.This leads to only a marginal increase in simulated s pre=20 .Even though the angle of repose, a M , is simulated in a range of 43°to 90°, simulating comparable bulk responses is infeasible in the ring shear test.Therefore, according to Fig. 11 we can conclude that an empty solution space is reached when the RC option is used.
Fig. 12 displays the simulation results of the 40 initial samples when the RR option, rotation restricted, is used.The same list of bulk responses as in Fig. 11 analyzed here, and therefore the feasibility is evaluated for N y = 3.
First, a range of s pre=20 bounded to 13.9 and 26.6 kPa is simulated; this covers both s exp.max and s exp.min .Second, a range of s 2:20 bounded to 2.5 to 6.5 kPa is simulated.This range covers both s exp.max and s exp.min .Third, a range of a M bounded to 60 and 90 is simulated; this covers the maximum and minimum values measured in the laboratory environment.Therefore, according to Fig. 12, a non-empty solution space is reached when the RR option is used.However, no sample satisfies all three calibration targets jointly.For example, sample 39 seems to be an optimal parameter set, however the simulated bulk responses compared to s exp,max(pre=20) , s exp,max(2:20) and a exp,max have errors, |e|, of 1.13%, 22.53% and 5.88% respectively.By establishing mathematical relationships between input variables and each calibration target, such errors can be minimized.For that reason, the RR option is used in the next steps as a feasible solution to be optimized further.
Step II: significant DEM variables A one-variable-at-a-time (OVAT) technique is used to create samples that allow us to investigate the direct effect of each DEM variable, x j , on simulated bulk responses by running a limited number of simulations.
Table 7 displays the samples created for this step, including 6 DEM input variables at the reference particle scale (S P = 1), when the RR option is used.This results in 60 samples in total, to be simulated in the ring shear and ledge angle of repose tests.When one variable is changed, the others are maintained at the displayed reference values.Reference values are based on one of the samples that was used in step I.In addition to 5 DEM input variables that were included in step I, the tangential stiffness multiplier, k t,mult., is also varied in this step.This allows us to check whether k t,mult.has any significant influence on the selected bulk responses.A similar list of bulk responses including s pre=20 , s 2:20 , and a M is analyzed in step II.Furthermore, larger ranges for the DEM input variables, compared to the previous step, are used to create samples.This allows us to run a comprehensive sensitivity analysis showing relationships between the DEM input variables and the selected bulk responses.Fig. 13 displays isolated effects of the 6 DEM input variables at S P = 2.25 on the simulated angle of repose.Since the ledge test box is performed in a rectangular container (as shown in Fig. 7), a M would be always equal or smaller than 90 .By varying coefficient of static friction, the maximum possible angle of repose, a M = 90 , being reached when m s,p-p !0.6.As expected based on the Mohr-Coulomb theory (Eq.( 2)), there is a positive strong correlation between m s,p-p and a M , as shown in Fig. 13a.A higher particle-particle friction results in a higher shear strength when normal pressure and cohesion strength are constant.By contrast, there is a negative correlation between G and a M , as can be seen in Fig. 13b.By increasing G from 1 to 128 MPa, a M decreases by around 20 .By increasing G, a lower contact overlap, d, is created.This is expected to result in lower forces in the adhesive branch of the contact spring (part III).Increasing G to higher values has negligible influence on a M .The ledge angle of repose simulations using k P equal to 0 and 0.99 result in unstable simulations, in which the stable situation (as discussed in Section 2.4) is not reached.As shown in Fig. 13c, by increasing k P from 0.1 to 0.5, a M decreases by around 20 , and further increasing k P has a negligible influence on a M .There is a strong positive correlation between Dc and a M , showing a non-linear trend near the extreme values (Fig. 13e).According to the Mohr-Coulomb theory (Eq.( 2)), a higher cohesion strength results in a higher shear strength.
Constant pull-off force and tangential stiffness multiplier are found to have negligible effects on a M in the investigated range, as shown in Fig. 13d and f, respectively.Coefficient of static friction, particle shear modulus, surface energy, and plasticity ratio are significant DEM variables influencing the angle of repose.Fig. 14 displays the results of the OVAT-based sensitivity analysis for simulated s pre=20 .There is a strong positive correlation between m s,p-p and s pre=20 .According to the Mohr-Coulomb theory (Eq.( 2)), the higher angle of internal friction of bulk material results in a higher shear strength when normal pressure and cohesion strength are constant.A linear trend seems to exist between these two parameters.The other 5 DEM input variables, compared to m s,p-p , have a weaker influence on s pre=20 .Particle shear modulus and surface energy have positive correlation values with s pre=20 .
The surface energy contributes in the cohesion strength of bulk material (denoted by c in Eq. ( 2)), thus contributing in the shear strength too.Fig. 15 displays the results of the OVAT-based sensitivity analysis for simulated s 2:20 .Coefficient of static friction has a strong positive correlation with s 2:20 , similar to its correlation with s pre=20 .The surface energy plays a more important role in s 2:20 , compared to s pre=20 .Increasing surface energy, Dc, from 0 to 25 J/m 2 causes an increase of more than 200% in s 2:20 .According to the Mohr-Coulomb theory (Eq.( 2)), at relatively low vertical pressure values, the cohesion strength, c, has a higher contribution to the shear strength, compared to shear flow at high vertical pressure values.
As expected, based on the results of the ledge of repose simulations, G has a negative correlation with s 2:20 .This is probably due a lower normal overlap created in the contact spring by increasing the value of G. Contact plasticity ratio, k p , also has some level of influence on s 2:20 , but not in a predictive manner .
In conclusion, only one input variable, k t,mult., has a negligible influence on the investigated bulk responses.Therefore, all the other 5 input variables are included in the surrogate modelingbased optimization in the next step.
Step III: surrogate modeling-based optimization In this step, first the Taguchi method is used to create multivariate samples to include variations of 5 significant DEM input variables when the RR option is used.Second, relationships between each calibration target and the DEM input variables are mapped to create F(X).This is done using the multiple linear regression technique.As discussed in Section 2.1, to consider the definiteness criterion, calibration targets are modified by excluding the ledge angle of repose test and by including W 80,65 and W 70,300 measured in the consolidation-penetration test.This means that four calibration targets are included in step III, and therefore M y = 4. Additionally, the maximum values of shear strength (shown in Table 2) are used as calibration targets in the simulation of a ring shear test.Third, an optimal set of DEM input parameters is found; these jointly satisfy the four selected calibration targets.
Table 8 presents the levels of the 5 significant DEM input variables at S P = 1 that are used to create multi-variate samples.Given the adequate simulated bulk responses in step I, m s,p-p is bounded to 0.2 and 0.4.For the same reason, levels of G are set to 2.5, 5, and 7.5 MPa.Three levels are selected for G to capture any possible non-linear relationship between G and the DEM calibration targets.k p is bounded to 0.2 and 0.6.This range is expected to be enough to capture a wide range of plasticity in the contact spring.Two other parameters, f 0 and Dc, which control cohesive forces in part III of the contact spring, are confounded.In other words, their levels are varied simultaneously in a way that allows us to minimize the number of samples.Thus, 4 coded variables are used in the Taguchi design to create samples.In total, 18 samples are created using the Taguchi method.
As investigated in [1], the reaction force on the wedge-shaped penetration tool is affected by the particle scaling factor.For that reason, the consolidation-penetration simulation is calibrated only for level of particle size (S P = 2.25), which is similar to the particle size used in the ledge angle of repose simulations.
Next, the matrix of simulated bulk responses, [Y'], including 4 different bulk responses for 18 samples, is created.This matrix is used to map relationships between DEM variables, X, and simulated bulk responses, Y'.Details of F(X) are presented in Table 9, including coefficients of the DEM variables in linear regressions fitted on simulated bulk responses, Y'.Cte.stands for the constant term in the regression model.Remarkably, in all the fitted linear regression models the coefficient of static friction has the highest level of significance.Values of coefficient of determination, R 2 , are also presented; in all the regression models, these are higher than 0.75.
Therefore, the multiple linear regression model is found to be adequate for us to continue with response optimization.If insufficient values of R 2 are reached in this step of the calibration procedure, either a higher number of training samples or more advanced surrogate modeling techniques can be used.
Fig. 16 presents an optimal set of DEM input variables that jointly satisfies four different calibration targets in step III with a composite desirability, d composite , equal to 0.61.Composite desirability, d composite , represents the geometric mean of individual desirability values, d, as shown in Eq. ( 5) and Eq. ( 6), respectively.
where f X ð Þis the predicated bulk response using the linear regression, and y is the target bulk response that is measured physically.y 0 min and y 0 max respectively represent the lowest and highest simulated values of a specific bulk response among all samples in step III.Each row in Fig. 16, except the top one, represents a specific simulated bulk response with its maximum possible d value obtained by finding an optimal set of DEM input variables.For example, the last row represents the response optimization for shear strength in the pre-shear stage, s pre=20 .For this bulk response, the physically measured value, y, is equal to 19.4 kPa.
Using the mapped relationship between DEM variables and y', simulated bulk response, a combination of variables is found that is predicted to lead to f(X*) = 18.7 kPa.This means that the outcome predicted in the simulation of a ring shear test using the current solution, shown in red, is a s pre=20 equal to 18.7 kPa, with d = 0.80.

Verifying the calibration procedure
This section discusses verification of the calibration procedure, step IV.First, we need to verify whether the outcome of surrogate modeling-based optimization is adequate.This is done by running simulations using the optimal set of DEM input parameters and comparing simulated bulk responses to predicted values, f(X*).Second, |e| mean is used to compare simulated bulk responses -using the optimal set -with all the calibration targets, corresponding to the maximum values in Table 2 and the target values in Table 3. Third, the entire yield locus in the ring shear test, including 1 level of r pre and 4 levels of r shear , is compared between the calibrated simulation and experiment.Fourth, the wall friction test as an independent bulk response is verified for various stress states.
First, ring shear and consolidation-penetration tests are simulated using the optimal set found Fig. 16.In Table 10, four different simulated bulk responses are compared with values predicted using the surrogate-based optimization.The relative difference is 10% in all cases, and therefore the adequacy of the multiple linear regression technique together with the response optimizer is confirmed for our DEM calibration problem.If large differences between y' and f(X*) had been captured, a higher number of samples or more advanced regression techniques could have been used to minimize the relative difference.
Second, |e| mean is used to compare simulated bulk responses -using the optimal set -with all the calibration targets, corresponding to the maximum values in Table 2 and the target values in Table 3.In other words, bulk density, shear strength, ledge angle of repose, and accumulative penetration resistance values are verified here.The shear stress in the pre-shear and shearing stages is simulated with |e| equal to 1% and 12.5% respectively.Bulk density values in loose and pre-sheared conditions, q b,0 and q b,20 , are simulated with |e| equal to 5.8% and 1.4%.On average, a relative deviation of 7% is captured in a ring shear test including four calibration targets.In the consolidation-penetration test, four different calibration targets are evaluated, including accumulative penetration resistance and bulk density values measured at two different pre-consolidation levels.In the consolidation-penetration test, accumulative penetration resistance parameters, W 80,65 and W 70,300 , are simulated with |e| smaller than 10%.Additionally, bulk density values at two different levels of r pre , 65 and 300 kPa, are simulated with negligible |e| values (smaller than 1%).This confirms that, using the elasto-plastic adhesive contact model, the calibration procedure was successful in capturing history-dependent behavior of the cohesive iron ore sample in terms of penetration Fig. 16.Finding an optimal set of DEM input variables that jointly satisfies calibration targets using response optimization.

Table 10
Comparing simulated bulk responses using the optimal set with predicted values of surrogate modeling-based optimization.resistance and bulk density.Finally, the ledge angle of repose, which was not used during the surrogate modeling-based optimization, is replicated with |e| = 7.1%.Therefore, considering simulated bulk density values in four different stress states and a M , the definiteness criterion is met using the optimal set of calibrated parameters, X*.Third, the entire yield locus is verified for the ring shear test conducted with r pre=20 .Fig. 17 compares the results of the ring shear test simulation using the optimal parameter set.Comparable shear stress values are measured in both simulation and experiment, with |e| mean = 6.7%.This verifies that the calibration procedure is able to replicate shear strength in various stress states and is able to capture the non-linear yield locus.Finally, wall friction measurements as a bulk response independent of the calibration targets are compared in Fig. 18, including 8 different stress states.The simulated wall yield locus shows a linear trend that replicates experimental values, with |e| mean = 5.5%.Since the Hertz-Mindlin (no-slip) contact model (without adhesive forces) was used to model particle-wall interactions, this linear trend could be expected.This finding is similar to the conclusion of [26], obtained by modeling a cohesive coal sample in a wall friction test.

Conclusions
This paper has successfully established a reliable and novel DEM calibration procedure by incorporating two important criteria: feasibility and definiteness.The DEM calibration procedure was applied successfully to model cohesive and stress-historydependent behavior of moist iron ore based on an elasto-plastic adhesive contact module.The definiteness of the calibrated parameter set has been verified using 20 different bulk response values in four test cases, such as ring shear, consolidation-penetration, and wall friction tests.
-The established calibration procedure can be used to calibrate material models when a high number of DEM input variables (e.g. 6) as well as a high number of calibration targets (i.e.> 2) are involved.-Both continuous and categorical variables can be used in step I, feasibility.Using the Latin hypercube design (LHD) method, it has been shown how a categorical DEM variable (i.e.rolling friction module) can be used during calibration.-During the calibration procedure, significant DEM variables can be screened using the one-variable-at-a-time (OVAT) method in step II.For ring shear and ledge angle of repose simulations, coefficient of static friction between particles (l s,p-p ) was found to be the most significant DEM variable.In general terms, this outcome is consistent with findings by other researchers [5].
Particle shear modulus (G), surface energy (Dc), and contact plasticity ratio (k P ) were the other significant variables when the elasto-plastic adhesive contact module was used.-In the current study, we have shown that surrogate modelingbased optimization is applicable when a high number (i.e.!4) of DEM input variables is involved.-The combination of Taguchi and multiple linear regression techniques was successful in the surrogate modeling-based optimization, with coefficient of determination values>0.75for all the calibration targets.
Further research is recommended to focus on, firstly, validating the calibrated model of the cohesive iron ore in simulating an industrial process (e.g.ship unloader grabs) where all the bulk responses (discussed in Section 4) play a role.Secondly, future researchers should apply the calibration procedure established here to other applications where high numbers of input variables and bulk responses are present.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Main components of a generic DEM calibration procedure.

(
A) Schulze ring shear test; (B) ledge angle of repose; and, (C) consolidation-penetration test.Additionally, three influencing parameters related to bulk properties were varied in the laboratory tests: (1) iron ore sample; (2) moisture content, denoted by MC; and (3) vertical consolidation pressure, denoted by r.The results obtained in the laboratory tests listed above (

Fig. 3 .
Fig. 3. Main steps of the DEM calibration procedure considering feasibility and definiteness criteria.

Fig. 4 .
Fig. 4. Schematic shear stress measurements in ring shear test, including pre-shear and shearing stages.

Fig. 7 .
Fig. 7.The ledge test box to determine angle of repose including dimensions: a) laboratory environment [19], side view; b) simulation environment, cross-sectional view.

Table 4
DEM input variables to model interaction between particles when RC option is used.

Fig. 10 .
Fig.10.Forty different samples for RR option at S P = 1, are created using Latin hypercube design for 5 variables.

Fig. 11 .
Fig. 11.Shear strength and angle of repose values captured in 40 samples when RC option is used: a) s pre=20 ; b) s 2:20 ; c) ledge angle of repose (a M ).

Fig. 12 .
Fig. 12. Shear strength and angle of repose values captured in 40 samples when RR option is used: a) s pre=20 ; b) s 2:20 ; c) ledge angle of repose (a M ).

Fig. 13 .
Fig. 13.Isolated effects of 6 DEM input variables at S P = 2.25 on the average angle of repose: a) coefficient of static friction; b) particle shear modulus; c) contact plasticity ratio; d) constant pull-off force; e) surface energy; f) tangential stiffness multiplier.

Fig. 14 .
Fig. 14.Isolated effects of 6 DEM input variables on the shear stress in the pre-shear stage (s pre=20 ): a) coefficient of static friction; b) particle shear modulus; c) contact plasticity ratio; d) constant pull-off force; e) surface energy; f) tangential stiffness multiplier.

Fig. 15 .
Fig. 15.Isolated effects of 6 DEM input variables on the shear stress in the shearing stage (s 2:20 ): a) coefficient of static friction; b) particle shear modulus; c) contact plasticity ratio; d) constant pull-off force; e) surface energy; f) tangential stiffness multiplier.

Table 5
DEM input variables to model interaction between particles when RR option is used.

Table 7
Sampling for step II, finding significant DEM variables.
Table 11 compares 9 different simulated bulk responses with their target values, which were measured physically using the laboratory tests.Four parameters in the ring shear test are compared, indicating shear strength and bulk density.

Table 8
Levels of DEM input variables at S P = 1 in step III: surrogate modeling-based optimization.

Table 9 F
(X) when i 2 M y ; mapped relationships between DEM variables and simulated bulk responses.

Table 11
Verification of calibration procedure; comparing simulated bulk responses with their calibration targets.