i. Indirect calibration

Methods for the calibration of DEM simulations can be divided into two broad categories: direct (or ‘bottom-up’) calibration and indirect (or ‘topdown’) calibration. In the former, the parameters of the particle contact model are directly estimated from particles’ microscopic properties measured us ing particle-scale experiments, such as measuring the coefficient of friction of a single particle, to give an example. In the latter, one or more bulk properties of the particulate are measured, and the microscopic properties of the individual particles back-computed from these measurements to determine the contact model parameters. This can be achieved either by using scaling laws mapping macroscopic to microscopic properties (⁶) or, as we will focus on in this article, by using ‘digital twins’ examples of which may be seen in Fig. 1 of the powder characterisation devices used to determine the set of model parameters that provide the same bulk behaviour exhibited in the ‘real’, physical system. In addition to direct and indirect calibration, we must also consider input - the process of determining the particle properties (e.g. morphology and density) we wish to emulate in simulation.

Figure 1 From top left to bottom right, digital twins of a Schulze shear cell, a Granutools GranuDrum, an FT4 powder rheometer, and a Granutools GranuHeap. Though in these representative images the particles depicted are spherical (unlike those used in a majority of real industrial processes), the models are fullycompatible with aspherical particle models. The opensource models of all four systems may be downloaded the University of Birmingham Positron Imaging Centre Github repository (https://github.com/uob-positronimaging-centre). 

When building a DEM simulation, the sizes and shapes of the DEM particles have to be defined. There exists a wide array of methodologies through which to estimate the sizes, and particle size distributions, of particles. These methodologies include sieving (⁷), image analysis (⁸), laser diffraction, dynamic light scattering and other light-scattering methods (⁹, ¹⁰), the electrozone method (¹⁰, ¹¹), and time-of-flight measurements (¹²), each of which carries its own distinct strengths and weaknesses which, for the sake of brevity, we will not discuss at length here.

A crucial limitation common to many of the aforementioned methods, however, is that the size of a particle is quantified using a single, scalar value, typically a sphere- or circle-equivalent sphere diameter (¹³), thus omitting crucial details regarding the particle’s geometry. While this may be sufficient for highly-spherical particles, for aspherical particles, key properties such as elongation (aspect ratio) (¹⁴), angularity (¹⁵) or nonconvexity (¹⁶) - which may significantly affect particle dynamics and thus should be represented in any rigorous DEM model - are not meaningfully captured. Nonetheless, widely-available (twodimensional) image analysis methods can capture averaged statistics relating to these properties, and the full, three-dimensional geometries of particles may be accurately captured using 3D scanning methodologies (¹⁷,¹⁸).

The accurate measurement of particle geometry is, however, only half of the battle, as the implementation of complex shapes in DEM simulations usually requires a rework of the contact detection algorithm, as well as other parts of the DEM model, relating to the mechanical forces and heat transfer between particles (¹⁹-²⁶). While this is not necessarily the case for ‘clump’ (multi-sphere/gluedsphere) methods (²⁷), such methods are generally considered not to be suitable for modelling highly angular particles. Implementing an accurate particle size distribution for spherical particles is, by contrast, mercifully straightforward - though simulating wide PSDs with high fidelity can prove computationally challenging (²⁸).

In the following sections we discuss the direct and indirect calibration of contact model parameters in further detail, outlining their major strengths and weaknesses, the practicalities of their application, and their general suitability for use in industrial DEM.

ii. Direct calibration

Direct calibration, while conceptually simple, can, in practice, be highly challenging. For the measurement of certain properties there exist wellestablished, standardised procedures, which may be simply performed using commercially-available systems. For example, particle density can generally be measured to a high degree of accuracy using gas pycnometry (²⁹).

However, other particle properties are distinctly more difficult to measure precisely than the sizes and shapes of particles. The coefficient of restitution (CoR), ε, provides an infamous example. Intuitively, one may expect this parameter to be fairly easy to measure - drop a particle from a known height, measure how high it bounces (off a solid surface or another particle) and use the ratio of bounce height to drop height to measure the energy lost, and thus the value of ε. In reality, however, there are numerous complicating factors to consider, of which we give here just a few of the most significant: the particle must be dropped in such a manner that it has no ‘spin’ either before or after collision, as in most set-ups energy lost to the particle’s rotational modes cannot be distinguished from energy lost due to restitution; the collision surfaces must be clean, dry and as close to perfectly smooth as possible so as to avoid imperfect contacts or energy loss to other sources; particles must be suitably large, dense and dry that cohesion effects - which cannot be easily distinguished from the effect of the CoR - do not influence results; a large enough number of repeat experiments must be performed to ensure adequate statistics. For aspherical particles, measurements are more complex still (³⁰): though the shape of a particle does not influence the inherent CoR of a material, if the particle used in an experiment such as those described above is not a perfect sphere (or, alternatively, a perfect cube) it is extremely difficult to perform said experiments without energy loss to the rotational modes, which may be difficult to decouple from energy loss due to restitution.

Despite these difficulties, several researchers have made excellent attempts at characterising and measuring the restitution coefficients of diverse materials (³¹-³⁵), including the particularly challenging case in which particles or surfaces are wetted (³⁶,³⁷). However, the time-intensive nature of the experiments, the lack of commercially-available equipment with which to perform them, and the aforementioned limitations on particle size, make such measurements unsuitable for many industrial applications.

The coefficient of friction is somewhat easier, though still by no means trivial, to measure directly. Commonly-used approaches for such mea surements include microtribometry (See Fig. 2) (³⁸-⁴⁰) and simple planar angle of friction tests (⁴¹,⁴²). The former involves the use of a robotic arm to apply a pre-defined normal force (Fn) to the material of interest (typically a single particle), which is glued or otherwise attached to said arm, whilst translating it horizontally across a surface and measuring the tangential force (Ft) produced, thus allowing the (sliding) frictional coefficient, µs, to be calculated based on Coulomb’s law (Ft ≤ µsFn). The latter test is conducted by affixing a monolayer of particles to a pair of flat plates which are placed atop one another. They are then tilted until the plates slide apart allowing a sliding friction coefficient to be estimated as µs = tanθ. Both methods can be used to measure both particle-particle and - where a suitable sample of the relevant wall material is also available - particle-wall sliding friction coefficients.

Figure 2 A microtribometer system (left) and its application measuring the particle-wall coefficient of friction of a coffee bean by measuring the force required to slide it across a flat, horizontal surface under a known normal load. 

In the context of DEM calibration, there remain questions as to whether the frictional coefficients measured using the above techniques correspond precisely to the coefficients implemented in DEM, and thus whether they can be used safely for direct calibration. Certainly, prior work comparing a ‘real’ and simulated planar angle of friction tester has shown that there is not a one-to-one correspondence (⁴¹). The discrepancies observed between the experimental and numerical cases likely arise from differences between the perfect contact assumed in DEM and the real-life situation, where microscopic surface roughness will typically yield an imperfect contact between any two macroscopic objects. This discrepancy could potentially be reconciled by measuring the surface roughness of said objects to determine an effective contact area though - as is so often the case with DEM calibration - this is not a trivial task.

The coefficient of rolling friction, µ _r - a parameter relevant to spherical particles, akin to sliding friction but, as its name suggests, applied to the relative rotational motion of particles - may also be measured using a similarly simple test. The ramp rolling friction test (⁴³,⁴⁴) involves rolling a particle down a smooth ramp of height h and onto a smooth, level substrate. The coefficient of rolling friction may then be determined simply as the ratio /, where d is the distance travelled across the substrate before coming to rest. Similarly to the previously mentioned CoR tests, the results of the ramp rolling friction tests become difficult to interpret when adhesive effects are likely to be significant. As such, while such tests can be useful for measuring the properties of larger particles, they become less useful for smaller (and otherwise more adhesive) particles, limiting their value in many industrial applications. The tests are also once again reliant on the use of highly clean, smooth surfaces. The cohesive and adhesive properties of particles can also be measured in a relatively straightforward manner using, for example, atomic force microscopy (⁴⁵-⁴⁸); once again, however, such measurements require careful preparation and choice of equipment (⁴⁵,⁴⁷), as well as expertise in the operation thereof.

While the above certainly does not represent a comprehensive list of direct measurement techniques, it nonetheless demonstrates a clear running theme: that such measurements, though conceptually simple, require an extremely (in some cases almost unfeasibly) careful execution in order to produce meaningful results. It thus illustrates an important point, namely that in order to fully calibrate a material using direct measurement techniques, one requires both a significant amount of specialised, often not-commercially-available equipment, and a significant amount of time and resources to execute the measurements correctly and produce adequate statistics.

It is finally worth noting the common - though somewhat contentious - suggestion that due to the key simplifications associated with numerical models (DEM included), direct calibration cannot be expected to accurately reproduce the bulk behaviours of the systems modelled. To put it differently, the advantage of a precise evaluation of microscopic particle parameters can be lost in the simplicity of the contact model - even if we put in the effort to rigorously and accurately measure the frictional coefficient and CoR, the resultant simulation often has little correlation to the real system. Certainly this is not true in all cases - for example, in the 2003 work of Xu et al. (⁴⁹), good, quantitative agreement was achieved between low-energy microgravity experiments, kinetic-theory-based models, and hard-sphere DEM simulations calibrated using exclusively direct methods. However, for less idealised systems (notably dense, high-stress systems where considerable deformation may be expected and thus the fundamental assumptions of DEM begin to break down) the efficacy of direct-calibration methods has not, to the authors’ knowledge, been rigorously proven. Further research into a) the effectiveness of direct calibration methods under diverse system conditions, and b) the development of more detailed models of the physical mechanisms underlying particle interactions, would both represent highly valuable pursuits.

From the previous section, it is apparent that while the full calibration of a DEM simulation using exclusively direct methods is possible under suitable conditions, for most industrial applications - where particles are often small, cohesive and geometrically complex, time and manpower are finite, and a wide range of distinct materials is typically used - such an approach is unfeasible. Indirect calibration provides an alternative approach which overcomes some of these limitations, but nonetheless introduces complexities of its own.

A common approach to indirect calibration in industry can be broadly outlined as follows:

Perform one or more characterisation tests on the particulate of interest, thus determining one or more macroscopic responses.

Perform multiple simulations of the characterisation process using a digital twin of the characterisation tool(s) used, implementing various combinations of values for the microscopic parameters to be calibrated.

For each simulation performed, determine the macroscopic response(s) measured in i., and compare the equivalent simulated value obtained to the experimental measurement.

For a suitable set of experiments and simulations, the set of microscopic parameters producing the closest agreement between simulation and experiment can be taken as the ‘true’ values.

A simple example would be to use an angle of repose tester to measure the experimental angle of repose (AoR), φ _exp. , of a given particulate, perform a set of simulations with a DEM model of the same angle of repose tester using various values of, for example, the sliding (µ _s ) and rolling (µ _r ) coefficients of friction, and then determine the values of µ _s and µ _r for which the difference between the experimental and simulated values of the angle of repose, σ _φ = |φ _exp. − φ _sim. |, is minimal. Assuming all other particle properties are already known, the combination of µ _s and µ _r yielding a minimum in σ _φ should, hypothetically, represent the true values of these parameters. In reality, however, this is not necessarily the case - see, as an example, Fig. 3, which shows the value of σ _φ for 256 permutations of µ _s and µ _r . It is immediately obvious from this image that there exist numerous distinct permutations yielding extremely low σ _φ , suggesting that there are multiple possible combinations of µ _s and µ _r capable of recreating the experimentallymeasured angle of repose. Such an outcome should perhaps not be surprising, as in the current case we are, in essence, attempting to solve a problem with two free parameters using only a single objective. To put it in less abstract terms, as the AoR is dependent on both the sliding and rolling coefficients of friction, it makes intuitive sense that a decrease in µs might be compensated for by an increase in µr in order to maintain a single value of φ, thus giving rise to the ‘L’-shaped region of possible solutions observed.

Figure 3: Example of the variation with µ _s and µ _r of the absolute deviation σ _φ = |φ _exp. −φ _sim. | between the measured and simulated dynamic angle of repose for a simple, spherical, non-cohesive powder. 

Figure 4 In addition to simple angle of repose measurements, the free surface of a rotating drum (i.e. a dynamic angle of repose tester) in the rolling or cascading regimes can be characterised by a third-order polynomial, providing additional information which can be used to back-calculate microscopic particle properties. 

The problem outlined above is not insurmountable, however. The yield loci produced by shear cells (⁵⁰-⁵²), or the torque vs. depth profiles of the Freeman FT4 powder rheometer (⁵³), for example, produce distinct functional forms dependent on the frictional, cohesive and geometric properties of the particulates examined. As such, they each provide more than one objective/constraint, allowing them to potentially determine multiple free parameters simultaneously. Indeed, even for the case of a (dynamic) angle of repose tester, one may solve such problems by considering not only a single-valued angle of repose, but the full free-surface profile of the particulate flow (see e.g. Fig. 4). For highly multi-dimensional optimisation problems - that is to say, if we wish to simultaneously calibrate a significant number of distinct parameters using an exclusively indirect approach - it will likely become necessary to use multiple characterisation devices in order to provide the necessary number of objectives/constraints.

Increasing the number of parameters to indirectly calibrate also introduces additional complexities beyond simply the need to perform more characterisation tests. In the example depicted in Fig. 3, we wish to determine the values of two DEM parameters, meaning that we must, in effect, solve an optimisation problem in a two-dimensional (µ _s µ _r ) parameter space. Even in this comparatively simple, two-dimensional example, we have performed N _s = 16² = 256 simulations. If we wanted to calibrate also the restitution coefficient, for example, and to explore the new, three-dimensional phase space to the same degree of detail as the original 2D case, then we would need to perform instead N _s = 16³ = 4,096 simulations. To calibrate also the cohesive energy density would thus require N _s = 65,536. If we additionally wanted to separately calibrate particle-wall and particleparticle values for µ _s and µ _r then we would need to run a staggering N _s = 16,777,216 simulations which, even for an extremely well-resourced company, would likely prove an insurmountable challenge. As such, it is clear that in order to maximise the value of indirect calibration approaches, we need to utilise more intelligent schemes than the simple, brute-force approach exhibited here. Considerable increases in efficiency can be achieved, for example, by using standard Design of Experiment approaches (⁵⁴) to explore the parameter space more efficiently, and/or by reducing the size of the parameter space to explore using prior knowledge of the system (for example, if the material is known to be free-flowing, there is no need to calibrate the cohesive energy density, thus immediately reducing the dimensionality of the problem).

The development of intelligent methods and algorithms for the calibration of DEM models is, in fact, an active area of research. Researchers have proposed a variety of methods utilising, for example, advanced design of experiments (DoE) combined with simple optimisation strategies (⁵⁵-⁵⁷), artificial neural networks (58-60), Bayesian inference (⁶¹-⁶⁵) and evolutionary algorithms .

Even with an extremely well-designed DoE, however, there remains an additional problem to consider: to what extent do the tests truly characterise the material, and how relevant is a given characterisation test to the conditions of the real system of interest? For example, an angle of repose tester is likely to provide meaningful information regarding a particle’s frictional properties, but is largely insensitive to the Young’s modulus of a particle. Similarly, a shear cell may be expected to provide valuable insight into a particulate’s cohesive properties, but no meaningful information regarding its restitution coefficient. As such, it is important to assure that the characterisation tools used mimic, as far as possible, the system we wish to ultimately model - if we wish to model a high-stress environment, we should expose particles to high stresses in our calibration experiments; if we wish to model a rapidly-flowing system, a dynamic characterisation test should be employed.

Moreover, with indirect calibration the set of parameters of the contact model are determined such that the macroscopic response of the system matches that obtained from external characterisation of the real material. This allows one to overcome, to some extent, the limitation of the direct calibration related to the simplified particle shapes and size distribution used. Indeed, the information on the complex shape of the particles of the real material will be included and distributed into the calibrated model parameters. However, this is at the price of a loss of physical meaning of these parameters.

A final issue worth mentioning in relation to indirect calibration is that of precisely determining and modelling cohesion. Cohesion, and specifically the modelling thereof, is a matter of significant importance in the majority of particle-handling industries, yet one about which our present knowledge is decidedly lacking; indeed, while we have a reasonable understanding for the case of spheres, beyond this we know very little with any certainty. The ‘macro’ cohesion experienced by a given system is usually the result of an interplay between numerous particle properties and external forces, including (but not limited to) van der Waals, liquid bridging and electrostatics, as well as potentially geometric effects (e.g. interlocking) (⁶⁶,⁶⁷). Determining the individual contributions of these various factors through indirect calibration methods would be, at best, highly challenging. To account for this issue, DEM practitioners often use simpler (and thus less physical) force models (⁶⁸) to ‘lump together’ these contributions; while these simplified models are capable of successfully reproducing the dynamics of a given system, they nonetheless limit the physical insight one may gain. Nonetheless, in industry - where the primary aim is typically to accurately model a given system as opposed to obtaining a deep understanding of its underlying physics - such an approach is typically satisfactory.