Calibration of micromechanical parameters for DEM simulations by using the particle filter

The calibration of DEMmodels is typically accomplished by trail and error. However, the procedure lacks of objectivity and has several uncertainties. To deal with these issues, the particle filter is employed as a novel approach to calibrate DEM models of granular soils. The posterior probability distribution of the microparameters that give numerical results in good agreement with the experimental response of a Toyoura sand specimen is approximated by independent model trajectories, referred as ‘particles’, based on Monte Carlo sampling. The soil specimen is modeled by polydisperse packings with different numbers of spherical grains. Prepared in ‘stress-free’ states, the packings are subjected to triaxial quasistatic loading. Given the experimental data, the posterior probability distribution is incrementally updated, until convergence is reached. The resulting ‘particles’ with higher weights are identified as the calibration results. The evolutions of the weighted averages and posterior probability distribution of the micro-parameters are plotted to show the advantage of using a particle filter, i.e., multiple solutions are identified for each parameter with known probabilities of reproducing the experimental response.


Introduction
The Discrete Element Method (DEM) often brings forth new understanding that are difficult to acquire from sophisticated experiments [1].Notwithstanding the versatility, the DEM model of a granular material with appropriate micromechanical parameters generally requires a 'calibration' against the macroscopic experimental response.The conventional calibration employs 'one at a time' sensitivity analyses of the parameters to derive the micromacro relations.Depending on the contact laws and material types, various micro-macro relations are obtained.Among a handful of systematic calibration approaches, the design of experiments (DOE) methods are efficient in searching possible solutions in the multi-dimensional parameter space with a manageable number of DEM simulations and optimized outcomes.Hanley [2] applied the DOE to calibrate the DEM models of crushable agglomerates.The interaction between key parameters were effectively considered by the orthogonal arrays designed with the Taguchi methods.Nevertheless, the DOE requires the knowledge of interaction between parameters, which is generally unavailable for various types of granular materials.A significant limitation of the aforementioned approaches is that the calibration can only be conducted with instantaneous bulk behaviors (e.g., Young's modulus) rather than a complete history of stress-strain response.To e-mail: h.cheng@utwente.nldeal with the above-mentioned difficulties, the particle filter 1 (PF) and sequential importance sampling (SIS), which can jointly account for loading history, are selected for the current problem.The PF applies the recursive formula of the sequential Bayesian estimation and approximates the posterior probability density function (PDF) with the SIS algorithm.The proposed calibration approach for DEM simulations is expedient, because the PF is well-justified for nonlinear and non-Gaussian geotechnical problems.Moreover, both the PF and SIS can be easily implemented in the open-source DEM package YADE [3], where parallel DEM simulations are run as model trajectories for the Monte Carlo sampling.

DEM simulations of granular materials
The parameter values for conducting three-dimensional DEM simulations are sampled from the posterior probability distribution, given the sequentially obtained experimental data.YADE [3] models granular materials as packings of solid grains with simplified geometries and vanishingly small intergranular overlaps.It tracks the kinematics of each grain within the explicit time integration scheme, based upon the net force resulted from the intergranular forces.In the present work, the robustness of the PF for calibrating DEM models is examined by modeling a Toyoura sand specimen (50 × 50 × 100 mm) with a void ratio of 0.68 as a dense packing of polydisperse spherical grains.Among all the randomly generated dense packings, three (N g = 1000, 2000 and 5000) are selected, because of their relatively smooth and spherical contact orientation diagrams after the initial isotropic compaction.The Hertz-Mindlin (HM) and classical linear (CL) contact laws are employed to govern the intergranular mechanical behaviors along the normal and tangential directions.Rolling/twisting moment transfer is allowed between adjoining spheres to account for the effect of grain roughness.Note that these non-physical parameters are necessary to obtain bulk friction angles which are close to experimental values of dense sands, without extra computational effort to capture the kinematics related to irregular shapes.Both intergranular tangential forces and moments are assumed to be bounded by Coulomb type plastic limits.Both contact laws have five micromechanical parameters, namely, the Young's modulus E c and Poisson's ratio υ c of contacting grains, factors β m and η m which define rolling/bending stiffnesses and maximum moments, and intergranular friction angle μ.The DEM simulations of sand specimens under drained triaxial compression (DTC) are performed in a strain-controlled quasistatic manner.
3 Data assimilation for granular materials with the particle filter

State space model and state estimation
Given a specimen of granular soil and its equivalent DEM packing, the physical measurements on the specimen and the prediction resulting from the DEM packing can be described by a nonlinear and non-Gaussian state space model: where the state vector x t consists of stress ratio σ a /σ r , radial strain ε r and volumetric strain ε v at a discrete loading step t, whereas the observation vector y t was directly measured in the DTC experiments [4]; υ t and ω t are the system and observation errors whose PDFs follow normal distributions with zero means.f t denotes the operator of the nonlinear dynamic model (i.e., DEM model), whose current state depends on the entire loading history.In the current problem, the nonlinear observation operator h t is reduced to an identity matrix of size three.

Particle filter
The PF approximates PDFs via a set of particles (referred as an ensemble) and the associated weights.From the ensemble {x (1)  t−1|t−1 , x (2)  t−1|t−1 , ..., x (N) t−1|t−1 } and the weights {w (1)  t−1 , w (2)  t−1 , ..., w (N) t−1 }, the filtered distribution p(x t−1 |y 1:t−1 ) at time t − 1 is approximated as where N is the number of particles, δ is the Dirac delta function, and the superscript (i) indicates the ith particle x (i) t−1|t−1 and its associated weight w (i) t−1 .It should be noted that all weights are greater or equal to zero and the summation of them must be one.With the ensemble sampled using the Monte Carlo Method, the recursive formulas for approximating the one-step-ahead forecast and filtered distributions p(x t |y 1:t−1 ) and p(x t |y 1:t ) are simplified as follows: One-step-ahead forecasting with If the observation system is linear and p(y t |x t|t−1 ) of the ith particle reads (7) where m is the dimension of x t , R t is a predetermined covariance matrix for the observation error, and H t is the observation operator in matrix form.Given the previous weight w (i)  t−1 and current regulated weight w(i) t of each particle (Equation 6), new weight w (i)  t can be updated as

Sampling and filtering procedure
The loading step t ranges from one to one hundred.For the current simple problem with few parameters, a total number of 2000 particles N p is sufficient to achieve good accuracy.The readers are referred to [5] for further details.The proposed calibration procedure is summarized as follows: 1. Initialization: Generate particles {x (1)  0 , x (2) 0 , ..., x (N) 0 } from the initial distribution p(x 0 ) by initializing parallel DEM simulations with a low-discrepancy sequence in R n (n is the number of parameters to be identified).Predict the state of each particle x (i) t−1 based on the parallel DEM simulations run for the t − 1 loading step, at which the observation data y t−1 are obtained.

Filtering:
Calculate the weight of each particle w (i) t that expresses the 'fitness' of the prior particles to the observation data y t with Equation 8 and Equation 6.

Weight updating:
The filtered distribution p(x t |y 1:t ) is approximated by the particles and the updated values of the associated weights.Return to Step 2 with p(x t |y 1:t ), {x (1)  t , x (2) t , ..., x (N) t } and {w (1)  t , w (2) t , ..., w (N) t } for the prediction and filtering at the t + 1 loading step.

Calibration of micro-parameters for DEM simulations
A Halton sequence generated with the means and standard deviations in Table 1 is used as the particles (N p = 2000).Each particle, which represents a specific combination of micro-parameters, is assigned to the contact law (HM or CL) for the three initially prepared 'stress-free' packings (N g = 1000, 2000 and 5000).2000 DEM simulations of DTC tests are then conducted in a two-phase loading program: the packings are first isotropically compressed to a prescribed confining stress (σ c = 0.2 MPa), and then sheared under triaxial compression in a quasistatic manner.The state variables are extracted at each loading step where the measurements were made.In what follows, the respective evolutions of the weighted averages of each micro-parameter obtained from different contact laws and DEM packings are first shown to examine the efficiency of the PF.The posterior PDF at different simulation steps is also plotted to further demonstrate the features of the PF.

Weighted averages of micro-parameters
The values randomly generated for each parameter are weighted by the approximate posterior PDF at each loading step.The resulted evolutions of the weighted averages of the micro-parameters are illustrated in Figure 1.While updating the weights through the entire loading history, the PF eventually results in less pronounced variation of the posterior PDF, which can also be understood from the weighted averages converging into the true values.The combinations of micro-parameter values that are associated with higher weights at the converged steps are identified as the calibration results.Although, convergence is reached for both contact laws and all three DEM packings, the convergence rate appears to be lower with larger fluctuations in the cases of DEM models governed by the HM law (dashed curves).Note that the true values which the weighted averages converge to are very close for the three packings governed by the same contact law.The agreement is almost perfect in the cases where the CL law governs intergranular behaviors (solid curves).The converged weighted averages of μ and η m of the HM and CL laws (with moment transfer) differ slightly, whereas the differences between the calibration results for E c , υ c and β m of the two contact laws are more evident.This is because the same Coulomb type friction is applied on tangential forces and contact moments, despite that the normal and tangential stiffnesses are calculated from different contact laws.

Posterior PDFs of micro-parameters
The advantage of the PF in a parameter identification problem is the ability to approximate the posterior PDF in a multi-dimensional parameter space.By taking advantage of the PF and SIS, the updated distribution can be understood at each loading step, as shown in Figure 2. Subplots (a)-(e) and (f)-(j) show respectively the posterior PDF of the five micro-parameters of the HM and CL laws at five consecutive loading steps with an increment of 20 steps, given the triaxial response of Toyoura sand.
It can be seen that the initial distribution along the E c and μ axes gradually shifts from a uniform to a Gaussianlike distribution.It is surprising to find that the preferable values for E c and μ are already identifiable within the first 20 loading steps (ε a = 2%).The distribution associated with υ c appears to exhibit a bimodal character at the 40th loading step (ε a = 4%) which corresponds to the peak stress state.However, the weights are only assigned to few particles after the 60th steps.A further refinement of the range of υ c would help to better represent the bimodal distribution.Although suitable values for β m and η m are identified, the associated posterior distributions appear to be more scattered and do not evolve consistently, even in the early stage of the calibration process.The fact that the distributions can hardly be described by a Gaussian or mixtures of Gaussians might be due to the non-physical nature of these rolling/twisting parameters.
Similar to the converging weighted averages of the micro-parameters as shown in Figure 1, the posterior PDFs collapse for the DEM packings with different numbers of grains (N g = 1000, 2000 and 5000), governed by the same contact law.This finding verifies a scaling rule which facilitates the calibration of DEM models, that is, a minimal number of DEM grains can be used for fast calibration.For DEM simulations which require a large number of grains, same parameters can be used, as long as the initial void ratio and stress state are the same as in the simulations run for the calibration.For the current parameter identification problem, 2000 simulations were conducted for each DEM packing.It took less than two days to run 1000-grain quasistatic DEM simulations on a 16-core workstation.

Simulation results versus measurement data
The simulation results reproduced with the identified parameters (particles with the highest weights) are plotted with the experimental data in Figure 3. Good agreements are achieved, regardless of the different contact laws and numbers of grains applied in the models.While the PF is proved to be an effective tool in calibrating the DEM models, further study is needed to investigate the capacity of this method on the DEM modeling of granular materials at various stress states and more complex situations.

Conclusions
A new calibration method is proposed for the DEM simulations of granular soils.The PF and SIS algorithms are implemented to work with the open-source DEM package YADE.The micro-parameters for the DEM simulations of granular soils under triaxial compression are successfully calibrated with the proposed method.The posterior probability distribution is an important feature which is approximated by the PF.It helps to gain insights into parameter identification problem.It is found that the converged posterior PDFs collapse for different DEM simulations governed by the same contact law.A scaling rule is verified and can be adopted as a guideline for calibrating DEM models in a cost-effective manner.

Figure 1 .
Figure 1.Evolution of weighted averages of the microparameters versus loading step.

Figure 2 .
Figure 2. Evolution of posterior probability distribution associated with each micro-parameter.Subplots (a-e) are made with the HM law and subplots (f-j) with the CL law.

Figure 3 .
Figure 3.Comparison of experimental data and DEM simulation results reproduced with the identified parameters.

Table 1 .
Means and standard deviations of the particles.