Linearization-based methods for the calibration of bonded-particle models

In the work at hand, two methods for the calibration of the elastic material parameters of bonded-particle models (BPMs) are proposed. These methods are based on concepts of classical mechanics and enable a faster calibration compared to the conventional trial and error strategy. Moreover, they can be used to counter-check the consistency of the BPM. In the first method, the mathematical model is linearized and solved in a classical matrix-vector formulation similar to a finite element scheme. Further linearization in the second approach enables the direct calculation of the elastic stiffness tensor, reducing computational costs massively. To analyze the capabilities and limitations of both methods, they have been applied in three different case studies. Obtained results have shown that the new strategy allows us to significantly reduce the calculation time.


Introduction
For the numerical investigation of the mechanical behavior of different materials, the bonded-particle model (BPM) can be effectively applied. The BPM is an extension of the discrete element method (DEM) which was initially developed for modeling of granular materials. The DEM is a mesh-less approach, where the modeled material is represented as a set of interacting primary particles. In the case of the BPM, an additional type of elements, namely bonds, exists. Every two primary particles can be connected via bonds [24]. In recent years, the BPM has been widely applied to different materials and different length scales. The scales of the application are ranging from the indentation of silica glass with an indenter on the scale of micrometers [1] to modeling large-scale concrete samples in the dimensions of centimeters [21]. Also on the scale of a meter, the BPM has been applied, for instance, to simulate the soil displacement in tillage [16]. Furthermore, studies of the properties of microporous materials [11] can be found in the literature just like studies of the compression of titania agglomerates [14], the fracture of brick and mor-B Kolja Jarolin kolja.jarolin@tuhh.de 1 Institute of Solids Process Engineering and Particle Technology, Hamburg University of Technology, Denickestr. 15,21073 Hamburg, Germany tar materials [27], the crack coalescence of silicon carbide ceramics [12] or cyclic loading tests of pinewood chips [33].
One of the main advantages of the BPM compared to the finite element method (FEM) is the high flexibility in the specification of the material microstructure and the straightforward way to mimic material fracture. Due to the discrete nature, i.e., every single bond can have unique geometrical and material properties, it is possible to generate anisotropic multicomponent materials. Moreover, every bond can have a different rheological behavior with a unique set of parameters.
There are two main challenges related to the application of the BPM: -high computational efforts caused by a large number of discrete objects and small simulation time steps; -non-trivial estimation of unknown model parameters.
The first problem can be partially solved using modern computer architectures. The discrete nature of the approach even allows DEM models to scale very well on graphic processor units (GPUs). For instance, Govender et al. [9] show how a GPU implementation can be achieved. To deal with the second problem, no clear straightforward method exists and novel strategies need to be developed.
In the BPM, the modeled material is represented as a set of spherical primary particles. Rarely, it is possible to generate a BPM model with a microstructure which is almost identical to the reality [29]. Nevertheless, such applications do exist. For example, Jauffrès et al. [11] simulated realistic representations of the microstructure of partially sintered materials. However, in most cases, significant deviations between the real material microstructure and the BPM representation exist. This means that even with exact knowledge about the properties of the binding material, additional adjustments of model parameters are needed.
Most widely applied calibration strategies are based on the trial and error principle. The corresponding experiment, e.g., uniaxial compression or dynamic impact tests, is simulated with varied model parameters, and the obtained numerical solution is compared with the experimental data ( [5,28]). The main goal is finding an optimal parameter set that minimizes the deviations for different scenarios. The parameter variation can be performed manually or by applying design of experiments (DOE) methods. In recent years, different DOE-based methods have been proposed to calibrate BPMs. Yoon [34] used a response surface analysis and the Plackett-Burman design to estimate the parameters for bonded rock particles. Hanley et al. [10] applied the Taguchi method for DOE and used it for bonded disks in 2D and spheres in 3D. Chehreghani et al. [2] performed calibration with a response surface methodology and a central composite design approach. The adaptive moment optimization technique has been applied by Qu et al. [26] to calibrate a BPM. Also, strategies based on neural networks or genetic algorithms can be found in the literature [3,19,28].
However, parameter variation in DOE-based methods requires a large number of simulation runs. Although for the calibration of a lot of parameters only small deformations need to be simulated, the calibration process is very costly from a computational point. To overcome this problem and reduce costs for the calibration, we present new linearizationbased approaches.

Theory
The bonded-particle model is a meshfree approach where the physical object is represented as a set of interacting primary particles connected with bonds. During the simulation, particles as well as bonds are treated separately. Hence, the BPM is a very flexible approach. The model can be generally subdivided into three main sub-parts: Structural model Describes the internal structure of the modeled material. Hence, it is related to the spatial material distribution, positions, and radii of primary particles (nodes), connections between primary particles, bonds radii, etc.; Functional model Describes the response of the material to deformation caused by external events. In particular, these are forces and moments acting in the bonds between particles and between particles and walls. Moreover, the functional models contain specific criteria to describe the breakage of the bonds; Parameters of functional model Material-dependent parameters to describe the stiffness, strength, plasticity, softening, etc.
Depending on the type of the modeled material, the structural model can exactly describe real physical microstructure consisting of spherical particles and cylindrical bonds [5,29] or it can be used only as an approximate representation [6,14]. For example, to generate optimal structural models, Spettl et al. [29] have proposed a bonded-particle extraction approach. However, if the real structure cannot be directly reproduced, some sort of homogenization must be used. In this case, the particles are randomly placed in 3D space. Afterward, using a force-biased algorithm, the particles are iteratively shifted to minimize interparticle overlaps and bonds between them are generated. In Fig. 1, exemplary applications of this approach to represent homogeneous [6,14] or heterogeneous [7,32] materials are shown.
For the functional model, mainly three different models can be defined: one describing particle-wall interactions, one describing particle-particle interaction, and one describing interaction via bonds (see Fig. 2). Depending on the desired behavior, bonds can replace the particle-particle interaction or they can superimpose the model for particle-particle interactions. Both concepts, i.e., superposition and replacement, will occur in the different cases considered later. The methods derived in the work at hand can be used for both concepts.
Particle-particle and particle-wall interactions exist also in the classical discrete element method. In the case study of this work, both interaction types are modeled with the Hertz-Mindlin contact model based on Mindlin and Deresiewicz [17] and Tsuji et al. [31]. For the interaction via bonds, different functional models can be found in the literature, starting from linear elastic material behavior [25] ending with cohesive damage-plasticity models [20]. In this contribution, a linear elastic model proposed by Potyondy and Cundall [25] is mainly used. Moreover, since the proposed strategy for the estimation of unknown parameters does not consider material failure and is applicable for small deformations only, the breakage of bonds is not considered. The increments of the forces in normal ΔF n and tangential ΔF t directions as well as torsional ΔM n and bending ΔM t moments are calculated in each time step as: where L init and A b are the bond initial length and crosssectional area, E and ν are the Young's modulus and Poisson ratio, I and J are the moments of inertia of the bond, v r and w r are the translational and rotational relative velocities, Δt is the simulation time step. These calculations are performed for each single bond i separately. Because every single bond can have its unique material parameters, a general functional model can have N b · N p degrees of freedom (DOF), where N b is the number of bonds and N p is the number of parameters in the functional model. To reduce the possible degrees of freedom, the following assumptions are made: -all bonds belonging to the same material component have the same Young's modulus E and Poisson ratio ν; -after the generation of the structural model, the structure is supposed to be fixed and is not varied during the procedure of parameter calibration. This means that bond's initial length L init and cross-sectional area A b are excluded from the procedure; -to simplify the estimation procedure the Poisson ratio is considered to be fixed and is excluded from the variation procedure for this work.
In summary, the number of parameters N p for the degrees of freedom is reduced to one, because only the Young's modulus E remains. Furthermore, instead of the number of bonds N b the first factor for the global DOF is now the number of material components N cmp . Thus, the number of DOF is equal to N cmp · 1.

Method
The developed methods aim to effectively determine the elastic model parameters of a BPM representing an agglomerate or continuum for which the mechanical properties are known. The mechanical properties, like Young's modulus or the Poisson ratio, are describing how the object reacts globally. However, the model parameters influence the local reaction of the primary particles. Therefore, the problem of calibration is an inverse problem, as we need to calculate the causal factors, i.e., the local properties, from an observation of the global behavior. The global behavior is, for example, the reaction force caused by a displacement in an uniaxial compression test. In general, the relation from local to global properties is not invertible and an optimization problem needs to be solved. Since calculating the response of the system with the DEM for a lot of different parameter sets is very expensive from a computational point of view, we are proposing two methods simplifying the calculation of the response. The proposed calibration methods for the parameter estimation are based on the assumption of small deformations. The first approach is based only on the assumption of a linear motion of the primary particles resulting in the well-known direct stiffness method (DSM). In the second variant, the additional assumptions of linear forces and locally independent stiffnesses are made. Therefore, a direct stiffness tensor estimation (STE) of the entire model is possible and no displacements need to be calculated.

Linearization of movements (direct stiffness method)
Analogous to the small deformation theory, we assume only small relative motions of primary particles and only a small deformation of the entire investigated object represented by the BPM. Hence, for a solid object, the individual particle motions are also small and can be linearized with an insignificant error. Applying this linearization to the structural and functional model of the BPM results in the direct stiffness method (DSM). This method originates from the calculation of trusses and is often used to derive the stiffness matrix used in the FEM from a mechanical perspective. The DSM is well described in the literature (e.g., Okereke and Keates [22], Pavlou [23]). Therefore, we will only give a brief introduction.

Application of DSM for solid bonds
To explain the application of the proposed approach, we will use it for bonds described with linear functional models given in Eq. (1)-(4) ( [6,25]). Assume a structural model, which consists of N primary particles (nodes) and only bonds transfer forces. Since we assumed the motion to be linear, the displacement of a particle can be described by a simple threedimensional vector. Furthermore, we are going to neglect the rotations of single particles, because the local rotation should be insignificant for small deformations of continuum objects.
Combining the displacements in all three dimensions for all particles into one vector gives: The relative motion of these particles causes the deformation of bonds and results in the formation of response forces. Since each particle can be connected to more than one neighboring particle, a set of response forces will act on it. For example, the force f nm acting on particle n due to a relative displacement u nm between particle n and particle m is a function f of this displacement Assuming small deformations, the function f is nearly linear for most bond models. In the local coordinate system of the bond, where the x-axis is in the direction of the vector connecting both particles, the relationship can then be written as a simple matrix-vector multiplication where K loc nm is the local stiffness matrix. For most cases, it will be a diagonal matrix of the normal stiffness k n and the tangential stiffness k s : If we sum all forces acting on each particle and write them in vector notation, similar to the previous formulation for the displacement, we obtain The relationship between all forces and displacements of each particle can also be written as a matrix-vector multiplication where K is the global stiffness matrix and K ∈ R 3N ,3N .
To assemble the global stiffness matrix, the local stiffness matrices must be rotated into the global coordinate system via a rotation matrix R with Subsequently, the rotated local stiffness tensor can be added to the global stiffness tensor according to the indices n and m.
By applying a certain load or prescribing the displacement of specific particles, the degrees of freedom can be reduced and a solution exists for This system can be easily solved via typical matrix inversion algorithms. Thus, the DSM allows solving the system directly for the entire displacement field under an external load. It replaces integrating the equations of motion, needed in the DEM, with solving the matrix equation (9).

Extension of DSM for particle-particle interactions
In the next step, we want to extend the approach to treat particle-particle interactions and to consider nonlinear functional models. One example of such a model is a Hertzian contact law which is widely used in DEM calculations ( [35]). From the perspective of the DSM, the Hertzian contact law is a nonlinear function h(u) of the displacements u. Thus, Eq. (7) can be extended to As we assume only small motions, the list of contact partners for each primary particle is considered to be nearly the same during the deformation. Therefore, we can speed up computations and determine the neighboring particles only once in advance and form a list of potential contacts. Nevertheless, the condition, if neighboring particles are actually in contact, is still during the computation. Furthermore, we assume no-slip condition and that the entire tangential displacement between any two contacted particles is accumulated in the tangential overlap. Note that we also assume that the direction vectors are constant, i.e., the normal and tangential direction between two particles is constant.
Finally, a fixed size nonlinear vector equation can be solved for the displacement u. A possible algorithm for solving the equation set is, for example, a Jacobian-free Newton-Krylov method (Knoll and Keyes [13]). Since the iterative numerical technique is applied for solving of Eq. (10), a convergence criterion must be defined. As a necessary criterion, a sufficiently low residual should be defined. For a calibration process, one can also define a criterion based on the change of the stiffness, since it might converge faster than the global residual.

Extension of DSM for particle-wall interactions
For many applications, the interaction between solid walls and the investigated object plays an important role and might be significant during calculations of the reaction force. For example, in the case of a curved surface or if the particles are very soft, giving a fixed displacement for the outer particle layers which are in the contact with geometries might not be sufficient. In these cases, an application of a moving imaginary wall is proposed.
Since we assume the wall to be ideal, i.e., no friction exists, only normal forces are acting between the wall and the particles. Therefore, the wall forces can be calculated from a nonlinear function of the displacements. The vector of all forces is then added to the equation system.

Linearization of forces (stiffness tensor estimation)
In the previous section, we showed how the linearization of movements allows to replace the time integration needed in the BPM with solving a matrix equation. In the following, we extend this technique by linearizing also forces and interactions. The resulting approach can directly estimate the stiffness tensor of the entire model. Hence, we refer to this method as stiffness tensor estimation (STE). The mathematical derivation of the method is based on Misra and Chang [18] and will be explained in the next subsections. In the following, Einstein notation for subscript indices will be used. The subscripts variables are in the range {1, 2, 3}, and subscripts appearing twice imply the summation of that term over all the values of the index. Superscripts will still be used to refer to the global structural model.

Main concept
Consider Hooke's law for a single primary particle n σ n i j = C n i jkl n kl (11) and the stress-strain relationship for a spherical particle n [30] with where m goes over all particles m interacting with the original particle n, d nm i is the branch vector connecting the particle centers and f nm i is the force between the particle n and the interacting particle m. V n , in general, is the volume of the particle, but in case of the BPM, it might also be the volume of the sub-domain the particle is representing. Assuming now a linear force-displacement relationship for each bond or for each particle-particle interaction gives where u nm j is the relative displacement and K nm i j is the stiffness of the bond or interaction between particles n and m. If we can then find a linear equation for the displacement u j and the global strain n kl u j = A jkl n kl , then, by combining Eqs. 11-14, it is possible to calculate the local stiffness tensor C n i jkl with By, for example, simple averaging the particle values we get an estimate of the object's stiffness tensor C i jkl . The derivation of the tensor A jkl is explained in the next section.

Derivation of linear relationship tensor A jkl
The following derivation is mainly based on Misra and Chang [18]. The first step in the derivation is to write the local displacement u nm i between any two particles n and m in terms of the overall average displacement, denoted with¯, and the local fluctuations from the overall average denoted with˜, d nm j is the connecting vector of the particles centers, ω is the rotational displacement, e i jk is the three-dimensional Levi-Civita symbol, r k is the particle radius, and φ n i is defined as The average strain tensorε ji is expressed as Following this notation and considering the dependencies in Eqs. 13 and 17, the balance of forces f = 0 for the particle n with neighboring particles m can be written as The same equation is true for the particle m with neighbors p. By changing the indices accordingly, we, therefore, get Typically, three different interactions are needed for m K nm i j to be invertible. For most models, this might only cause problems at or close to the boundaries. A simple solution for such cases is to neglect these particles, as their effect on the overall properties should be neglectable. If they are relevant for the case or too many exist, the method should not be applied, because it assumes a continuous and not a discrete reaction.
The same transformation as for particle n can be achieved for the particle m, i.e., The displacement fluctuations are thereforẽ φ n j = Γ n jkl +Γ n jkl ε kl and (28) For the displacement, we can write and therefore where δ jl is the Kronecker delta. Finally, from Eqs. (31) and (15) Therefore, the basic procedure to compute each stiffness tensor is:  (27)). The computation of all the terms follows the same procedure as written above for the simple example.

Incorporating wall interaction
As mentioned before, when calibrating very soft materials the stiffness of the interaction between the wall and particles can significantly reduce the total stiffness. Thus, in the calibration process, it is necessary to estimate the influence of the wall. One possible approach is assuming a simple spring model. Based on the target displacement, it is possible to calculate the resulting force and the necessary overlap to achieve the force. In the next step, the stiffness of the contact can be evaluated and the total stiffness is calculated as three springs in series. Furthermore, in the DEM a small overlap between the particles and the wall exists. In the DEM, this overlap results in the reaction force. Thus, the displacement of the wall is slightly bigger than the displacement of the particles in contact with the wall. For very soft materials, the overlap can be big enough to change the stiffness slightly. This has been considered when calculating the response of the system.

Case studies
To analyze the capabilities and limitations of the proposed strategies, a case study was conducted. In three different examples, the BPM was used to analyze the mechanical behavior of complex-structured materials under uniaxial compression. An overview of the analyzed case studies is given in Table 1 and Fig. 3.

Case 1
The first study is based on a previous work, where the mechanical behavior of biopolymer aerogels was analyzed with the DEM (see Dosta et al. [6]). In the contribution, three different alginate aerogel samples with varied CaCl 2 concentrations were produced in a cylindrical shape. Afterward, they were examined with quasi-static uniaxial compression tests. It was investigated how the variation of concentrations influenced the material stiffnesses. The material was assumed to be homogeneous, and therefore, randomly generated structures were used for modeling their behavior with the BPM. To obtain reliable simulation results, six different random structures for each concentration were generated and calibrated against uniaxial compression tests. The detailed description of the performed investigations can be found in Dosta et al. [6]. The final structural BPM model of the cylindrical aerogel sample consisted of about 10 5 particles and 4 · 10 5 bonds. A screenshot of the structural model is shown in Fig. 3.

Case 2
The second example is a set of bi-component agglomerates examined in Dosta et al. [7]. In this study, the analyzed materials consist of one type of primary particles and two types of solid bonds. Dosta et al. [7] investigated the mechanical behavior of the system in uniaxial compression tests for different stiffness and mixing ratios. Three different stiffness combinations (1:2, 2:1, and 1:4) with ten different mixing ratios varied from 10% to 100% were examined. For each case, ten different samples were randomly generated. The obtained structural models consist of 6000 particles and 29 000 bonds. A small cut of such agglomerate with two different types of bonds in the proportion of 40% and 60% is shown in Fig. 3.

Case 3
In the third example, the proposed calibration methods have been applied to a set of silica inverse opal structures. Such metamaterials show high specific strength and modulus and are widely used in different application areas. In this contribution, the macroporous silica structure is represented as a set of primary particles connected with bonds. The diameter of each pore is 756 nm, and the diameter of the primary particles in the BPM representation is 15 nm. The detailed description of this structure and its mechanical behavior can be found in [4]. Overall, samples with three different dimensions have been generated. This results in three Case 2 300 6 · 10 3 2.9 · 10 4 3.5 · 10 5 -2.1 · 10 6 Case 3 3 2 · 10 5 -8 · 10 5 2 · 10 6 -8 · 10 6 1.2 · 10 5 -7.1 · 10 5  BPM structural models consisting of about 2 · 10 5 , 4 · 10 5 , 8·10 5 } particles and {2·10 6 , 4·10 6 , 8·10 6 } solid bonds accordingly. One of the generated structures is shown in Fig. 3.

Performance analysis
To analyze the efficiency of the proposed calibration methods, all cases were explicitly solved with the DEM and the resulting stiffness of the whole sample was determined. Afterward, the proposed calibration techniques, i.e., the direct stiffness method (DSM) and the stiffness tensor estimation (STE), were applied to these cases. For the DSM, a fixed lower wall and a fixed displacement of the upper wall were used as boundary conditions to predict the final state of the DEM solution. For all cases, in the first step of the DSM, the pure bond system was solved iteratively with a target value of 1 · 10 −6 for the Euclidean norm of the residual. The resulting displacement field was used as initial solution for the nonlinear system. The nonlinear system was solved with a Jacobian-free Newton-Krylov solver. A residual smaller than 1 · 10 −3 has been chosen as abortion criteria. Based on the resulting force field, the reaction force was calculated. Afterward, the global stiffness of the sample was calculated with the displacement of the upper wall. For the STE, no boundary conditions or abortion criteria are necessary. It directly predicts the samples Young's modulus E. The stiffness was then calculated based on K = E · A/l. The cross-sectional area A and the length l of the sample were calculated from the body formed by the primary particles. For all three methods, DEM, DSM, and STE, the time needed to compute the final stiffness was measured.
Based on this procedure, two main performance characteristics have been analyzed: Speed-up factor is the ratio of the computation time for the explicit DEM calculations to the time needed for the calculations using the specific calibration technique, i.e., the DSM or STE Estimation error is the difference between the stiffness from the DEM calculations and the calibration techniques.
Both calibration approaches, the DSM and the STE were implemented in Matlab R2019a. The explicit DEM simulations were performed with the Musen simulation framework (see Dosta et al. [6]), which is implemented in C++. The calculations algorithms in Musen are parallelized for GPU using CUDA. The newly developed calibration methods were also implemented into Matlab using "gpuArrays" to see the basic performance trend for computing on the GPU. All calculations were carried out on an Intel Core i7-7820X processor with a quad-channel memory configuration and an NVIDIA GeForce GTX 1080 Ti graphic card.
To calculate the elasticity parameters in Case 1, a target strain of 2% was chosen. For this range, the material showed a linear behavior. Under large strains, i.e., strains exceeding 4%, aerogels reveal elastoplastic behavior [6]. Moreover, from the design of the structural model, which mimics high material porosity, it is known that for small strains almost no direct particle-particle interactions take place. Nevertheless, the DSM calculations considering nonlinearities (Eq. (10)) have been conducted.
The deviation of the stiffness due to the linearization and the solving process compared to the DEM simulation can be seen in Fig. 4. The figure shows that for Case 1 the DSM achieves an average deviation of 1% for all of the 18 samples. The linear version of the DSM has slightly higher deviations due to the linearized particle-wall interactions. The direct estimation using the STE technique results in a deviation with an average of 7%.
The calculation time needed to perform the direct DEM simulations and to apply the different variants of the proposed parameter estimation techniques is also shown in Fig. 4. The DSM implementation on the CPU is roughly 46 times faster than explicit DEM calculations. Due to the better optimization of the DEM code, the application of the DSM approach using GPU results in a smaller speed-up factor of 13. Neglecting particle-wall interactions increases both speed-up factors to 146 and 42 for the CPU and GPU implementation accordingly. The GPU version scales worse since the total computation time is already down to 4 s and the communication overhead using Matlab's gpuArrays becomes significant. Direct estimation with the STE method on the CPU takes only 1.5 s, while the duration of DEM calcula-tions is on average 1250 s. That results in the speed-up factor of 878. The version of the STE on the GPU takes also 1.5.
Contrary to the aerogel samples, in Case 2, the stiffness of the samples is strongly dominated by particle-particle and particle-wall interactions. Therefore, the reaction force is nonlinear, i.e., the stiffness is not constant over the displacement even for small strains. Nevertheless, we are applying the developed methods to test their limitations. The target strain is 2%. A residual smaller than 0.01 has been chosen as abortion criteria for the DSM.
As can be seen from Fig. 4, for this relatively small case, the DSM is about 40 times faster on the CPU. The speed-up factor is higher than in Case 1, because of the lower time step required for the DEM due to the stiffness. On the GPU, the DSM is only two times faster because the case is too small to benefit. The communication overhead becomes big enough that the CPU version is as fast or even faster than the GPU version in this case. A direct implementation in CUDA instead of Matlab should improve this, but optimizing for small cases was not considered here.
Due to the high nonlinearity of Case 2, the average estimation error for the DSM is about 6% and about 46% for the STE. The large speed-up factors of 2550 and 139 for the STE result from the small number of particles and the highly reduced complexity due to linearization. Since the assumptions of the STE are not fully full-filled anymore, the standard deviation increases and the maximum error reached 65%. However, for Case 2 the STE technique was still able to predict qualitatively the trend of how the stiffness changes with different bond stiffnesses. A major difficulty for the STE is to predict the tangential force of the particle-particle interaction. The cause is the nonlinear dependency of the stiffness in particle-particle interactions on the normal and tangential displacement, as well as the jump in stiffness in case of sliding. At the same time, although the forces are highly nonlinear, the assumption of linear motion stays accurate enough and the DSM can predict the reaction force with relatively small error. Moreover, in the given structural model, very few particles exist and, therefore, the reaction is rather discrete and not continuously distributed. Therefore, averaging becomes less accurate than, for example, for the samples of Case 1.
In Case 3, the reaction force of the functional model is nearly linear for the chosen deformation of 2%. From Fig. 4, it can be observed that despite the complex highly porous structure the DSM, as well as the STE, agrees very accurately with the results of the DEM. The estimation error is similar to Case 1. This demonstrates that both methods show consistent results for different scenarios. The speed-up factor on the CPU is also similar. In contrast, on the GPU the speed-up factor is a bit higher since Case 3 has more objects than Case 1. This results from the scaling of the methods, as we will discuss in the next subsection. The estimation error is shown as points in percent on the right axis. The error bars show the standard deviation for the cases with multiple samples Note that in case of the inverse opal structures another calculation step is needed for the STE. To get the stiffness from Young's modulus for the complex structure, some sort of mechanical model is needed. We used a simple model and assumed a set of springs in series. For this purpose, the domain was subdivided into multiple layers in the direction of the load. A layer height of two particle diameters was used to ensure sufficient particles inside the layer. The stiffness tensor of all particles in the layer was averaged based on the volume fraction. The stiffness of each layer was afterward calculated based on the cross-sectional area of the current layer. In the last step, the global stiffness was calculated as equivalent spring stiffness of multiple springs in series. As Fig. 4 shows, the approach results in relatively accurate results.

Scaling
Based on the third case, the scalability of the methods depending on the number of primary particles and solid bonds has been analyzed. In Fig. 5, the computation time on the GPU is plotted versus the number of particles for the different methods. The DEM computation time on the GPU increases by a factor of 2.2 for each model, while the number of particles increases by a factor of 2. In contrast, the DSM and the STE show an inconsistent scaling factor. In the first step, it is 1.9 and 1.7, respectively, and in the second step, it is 2.0 for both cases. The increase is probably caused by a constant computational overhead in the DSM and STE implementations. The overhead is becoming less relevant for larger cases, and the scaling factor converges to the true scaling. This suggests that the methods could become faster for smaller cases.
Since the scaling factor for the current implementation of the methods is smaller than the one of the DEM solver, the speed-up factor increases with an increasing number of particles. For example for the STE, it starts at 140 and ends at 190. Thus, the absolute time-saving increases for a higher number of particles. Note that already a constant speed-up factor would cause an increase of the absolute savings. As shown in Fig. 5, the methods reduce the computation time by about 950 s for the smallest model and 4600 s for the biggest model.

Three-dimensional distribution of elastic parameters
The STE gives, as a result, the contribution of each particle to the global stiffness tensor. To derive the stiffness of investigated material, we averaged over the entire volume of the sample. Instead, it is also possible to average only over a local sub-volume or the local particle neighborhood. As a result, a three-dimensional spatial distribution of the elasticity parameters inside the sample can be derived.
As an example, the cylindrical samples of Case 1 are used. As mentioned before, they consist of 10 5 particles and 4 · 10 5 bonds. For the evaluation, the deviation of the stiffness in each direction from the mean value is calculated with where i is the direction 1, 2, 3 or x, y, z, K local is the estimated stiffness at each particle position, and K is the arithmetic mean of all local stiffnesses. Since the method requires averaging to cope with fluctuations, the local stiffness of each particle is averaged over local volumes with a diameter of three particle diameters. The result is shown in Fig. 6. Deviations between −15% and 15% of the mean value occur locally. In the inner part of the sample, the stiffness is mainly higher than on the boundary. This is probably caused by the spatial heterogeneity of the coordination number. The particles situated in the outer layers have less neighbors as particles in the middle. Also, due to the relative coarse BPM representation some local minima and maxima can be identified.
The main goal of such a analysis is to verify whether the model's stiffness distribution corresponds to the expectation. Thus, the STE allows detecting defects in an early stage of the calibration process. Furthermore, the STE allows a fast calibration for anisotropic or other complex material behavior, like for samples with stiffness gradients.

Discussion
Calibration and validation of parameters for a bondedparticle model (BPM) is a time-consuming but crucial prerequisite for any further usage of the model. We presented an approach based on the direct stiffness method (DSM) and typical assumptions for small deformations. A comparable DSM-based approach for bond calibration was presented by Estay et al. [8]. They derive their method from the FEM and test the accuracy of the approach applying it for different particle diameter ratios. We extended this approach also to consider nonlinear particle-particle interactions and com-plex interactions with walls. Furthermore, we could show that our approach scales well for different sizes of structures and can be parallelized on GPUs. Since the method drastically reduces the computational costs, it allows us to efficiently calibrate large models on CPUs, which would need high-end GPUs to calculate with the DEM in a feasible amount of time.
The application of the proposed DSM approach for a bonded-particle model corresponds to performing a finiteelement analysis of a three-dimensional truss with a nonlinear stiffness matrix. Depending on the shape function and the mass matrix, one could define a structural model that can equivalently be solved with the FEM and BPM for small deformations.
To speed up the calibration process, we also developed an approach that can directly estimate Hooke's stiffness tensor of the model. Kruyt [15] showed already that a localadjustment field approach can be effectively used. Instead of applying a displacement field, we used simplifications shown by Misra and Chang [18] to allow a direct stiffness tensor estimation (STE).
As both approaches, DSM and STE, are based on continuum assumptions and are steady-state solutions, it is possible to use them to counter-check the BPM results. They can reveal problems with time-stepping or undesired wall influences. Furthermore, they can easily show how the load is distributed between bonds and particle interactions. Also, the STE can directly show local stiffness gradients or the stiffness distribution inside the object. Since we can predict such problems in the early stages, unnecessary calibration steps can be avoided. The direct estimation method also returns the stiffness for all directions. Thus, anisotropic behavior becomes easier to model as well. However, the stiffness tensor estimation is by design not capable of predicting accurately highly nonlinear behavior. Nevertheless, it still can show trends and give a rough estimation of mechanical response, which is its purpose. If higher accuracy is needed, the DSM must be applied.
Based on our results, we suggest starting a calibration with the stiffness tensor estimation. Afterward, one can use the direct stiffness method to further improve the accuracy of the predictions, especially if the model's reaction is nonlinear. In the last step, the DEM could be combined with an optimizer for fine adjustments and non-elastic parameters as shown in Dosta et al. [6]. x -direc on y -direc on z -direc on -15% -10% -5% 0% 5% 10% 15% Bonded-par cle model ( The deviation of the local stiffness in percent from the mean stiffness in each direction is shown (see Eq. (36)). The displayed model is one of the simulations of Case 1 (based on Dosta et al. [6]). The local stiffness is averaged over the local neighbors of each particle. On the right-hand side, an image of the structural model is shown for scale

Compliance with ethical standards
Conflicts of interest The authors declare that they have no conflict of interest.

Code availability Custom code
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.