A hybrid particle-geometric scaling approach for elasto-plastic adhesive DEM contact models

ThecomputationtimeofDiscreteElementMethod(DEM)simulationsincreasesexponentiallywhenparticlesize isreducedorthenumberofparticlesincreased.ThiscriticalchallengelimitstheuseofDEMsimulationforindus-trial applications, such as powder ﬂ ow in silos. Scaling techniques can offer a solution to reduce computation time. In this paper, we have developed a hybrid particle-geometric scaling approach with a focus on Elasto-Plastic Adhesive contact models. It established relationships between particle scaling factors and DEM contact input parameters. The isolated effects of varying particle size and geometric dimensions on bulk properties were also evaluated using uniaxial consolidation, static angle of repose, and ring shear tests. This paper shows how the particle scaling can be applied together with geometric scaling to incorporate two important aspects of bulk materials, their Elasto-Plastic behaviour and their cohesive forces.


Introduction
The Discrete Element Method (DEM) is a particle-based computational method to model the flow of granular materials and their interaction with equipment. To predict the outcome of the flow process accurately, the DEM parameters need to be chosen carefully. To select the input parameters with confidence, the common procedure is to calibrate and to validate DEM simulations [1][2][3][4]. Three different calibration approaches for choosing DEM parameters can be named, as follows: The first approach, direct measuring, requires measuring the input parameters directly at particle or contact level [5]. Accurate measurements of the micro-properties do not necessarily lead to a successful prediction of the bulk flow properties [2,3]. Furthermore, modelling the actual shape and size of particles leads to a computation time that is impractical for industrial applications, such as silo flow [4,5], transfer chutes [10][11][12], and ship unloader grabs [1], where the approximate number of particles required is greater than 10 7 .
In the second approach, in-situ calibration, field experiments on a specific industrial process, at a scale of 1:1, are used to either calibrate or re-calibrate a DEM simulation replicating the real process. For example, Ilic et al. [11] have used the in-situ calibration approach in a qualitative way for modelling the accelerated flow in transfer chutes. Using this calibration approach, the shape and size of particles can be modelled in a simplified way that are different than those actually handled in full-scale operations [12]. In the in-situ calibration approach, a sufficient number of experiments must be conducted to avoid ambiguity of DEM parameter set [13]. Additionally, a disadvantage of the in-situ calibration approach is that the calibrated DEM parameter set is dependent on the design, and it might fail to simulate processes different than the in-situ calibration experiments.
In the third approach, bulk calibration, a laboratory experiment or series of experiments is/are first conducted to measure the bulk properties that are relevant to the application under consideration [5]. Next, the input parameters for DEM simulations are calibrated by minimising the mismatch between output simulations and laboratory measurements. In general, to produce comparable bulk responses, these calibration simulations replicate a laboratory setup and procedures at a scale of 1:1 [14].
Although the bulk calibration approach can use a less detailed representation of particle shape and size, an important challenge remains the huge size of equipment used in bulk terminals, compared with particle size and the scale of calibration laboratory experiments. For instance, a Schulze Ring Shear Tester used to measure the bulk properties of particles smaller than 6 mm has an internal volume of less than 10 −2 m 3 [15]. Equipment in bulk terminals, such as ship unloader grabs and silos, by contrast, have volumes that are 10 3 to 10 4 greater than laboratory devices. This challenge limits the use of DEM simulations for industrial applications [9,10,16].
Various approaches are used to reduce the computation time of DEM simulations. They can be categorised into two main groups. In the first group, computational techniques are used to speed up simulations, whereas particle size and geometry are kept constant. For instance, in [8] reducing the stiffness of contact springs successfully led to more rapid simulation by allowing a larger integration time step to be chosen. [17] also proposed using a more efficient DEM solver to reduce the computation time.
In the second group, either geometric size or particle size or both is/ are scaled. Table 1 provides an overview of the scaling techniques applicable to DEM simulations. Scaling can be done at either global or local level. At global level, scaling is applied to the entire simulation domain, influencing all the particles. For example, particle size can be scaled up in the entire domain by a constant scaling factor s. At local level, scaling is only applied to a specific region or a specific group of particles. For instance, in the Scalping or Cut-off technique [6,15], finer fractions of particles are omitted by replacing them with the larger particle size fractions.
The "philosophy" of Local Particle Refinement, which uses upscaled particles outside the area of interest, is similar to Local Mesh Refinement [19]. This technique was successfully applied in [20] to model cone penetration into free-flowing materials using DEM. However, the main challenge when using local scaling techniques is that the speed of current DEM solvers depends mainly on the smallest particle size used in the simulation. The reduction in computation time is therefore limited by the critical integration time step.
In Exact Scaling, both particle size and geometric dimensions are scaled by the same scaling factor s. An example can be found in [21], which investigates the upscaling of the uni-axial confined, uni-axial unconfined and cone penetration test for cohesive elasto-plastic soils. An important uncertainty with the Exact Scaling method is that both micro-properties at particle scale (e.g. particle mass), and macroproperties at bulk scale (e.g. bulk volume, porosity) are varied simultaneously. For instance, to create comparable initial stress states during scaling, Janda and Ooi [21] suggested to reduce the gravity with the same scaling factor that is used to up-scale particles and geometry sizes. Since the scaled system should have the same energy density as the original (unscaled) system [22], altering gravity during scaling is not recommended. Schott et al. [8] also proposed the Geometric Downscaling of excavation equipment, which did not result in confirmed scaling rules.
Another technique, referred to as Coarse Graining (CG), substitutes larger grains for the original DEM particles, thus allowing for a lower number of particles in simulations. [1] successfully applied this technique, using the calibrated DEM parameters in a largescale simulation of grabs and free flowing bulk materials. [23] investigated the coarse graining of the JKR contact model [1] combined with Hertz-Mindlin [24] in a shear tester, but no scaling rules for modelling cohesive bulk materials were established. A more successful study [9] investigated the Coarse Graining of an Adhesive Elasto-Plastic in the uniaxial consolidation process through a trial and error approach, and it found that the constant cohesion force of the contact spring (constant pull-off force) should be scaled up by the square of the particle scaling factor.  1. Comparing the idea behind the hybrid particle-geometric scaling approach (a) with Exact Scaling (b) and particle upscaling (c).
One or more combinations of the scaling techniques mentioned can be used to reduce the computation time of DEM simulations. For example, [14] used three different scaling techniques, Exact Scaling, Coarse Graining and Scalping, to simulate the angle of repose of a free-flowing bulk material. In practice, however, raw bulk materials and powders, such as moist iron ore fines and coal usually show cohesive elastoplastic behaviour. Bulk responses of this type of materials, such as shear strength, bulk stiffness, and bulk density, depend on the history of applied normal pressure on the bulk specimen [25][26][27][28]. This stresshistory dependent behaviour can be simulated properly by using contact models that are based on an elasto-plastic adhesive spring [4,[29][30][31]. Thus the question is how combinations of scaling solutions can be applied to a different contact model incorporating the behaviour of cohesive elasto-plastic materials.
This paper, therefore, develops a hybrid particle-geometric scaling approach with the focus on an adhesive elasto-plastic DEM contact model. This hybrid approach combines particle scaling with geometric scaling in a sequential manner. Fig. 1 compares the idea behind the hybrid approach ( Fig. 1a) with Exact Scaling (Fig. 1b) and particle upscaling (Fig. 1c) techniques. In hybrid scaling, only particle properties (e.g. particle size, particle interaction parameters) or geometric properties (e.g. dimensions) are varied at a time, which is the main novelty over Exact Scaling. As discussed earlier, the main uncertainty with the Exact Scaling is that both microproperties at particle scale, and macro-properties at bulk scale are varied simultaneously. In hybrid scaling, the geometric size is scaled each time after applying particle upscaling, so basically larger particles can fit in the simulation setup. This creates a novelty for hybrid scaling over particle upscaling, in which up-scaled particles might not fit properly in the simulation setup. ϕT Computation is the ratio of the computation time of the scaled simulation to the computation time of the reference simulation.
Our particle scaling rules for cohesive elasto-plastic materials is based on extending the CG principles described in [1], which has been develop to scale up free-flowing elastic materials. In Section 2, we establish the relationship between the particle scaling factor and the contact settings for an Elasto-Plastic Adhesive model. Sections 3 and 4 go on to investigate both the decoupled influence of scaling of spherical particles and geometric scaling on bulk properties in the quasi-static and dynamic regimes. Following simulation setups are used that their internal volume are considerably smaller than equipment used in solid bulk terminals, such grabs, silos and hoppers: • Uniaxial confined consolidation test • Ledge angle of repose test; also referred as the shear box [14] • Schulze ring shear test By applying hybrid particle-geometric scaling on the above listed simulation setups, their internal volume can be increased to a level comparable to equipment in solid bulk terminals. This allows to create DEM simulations of the mentioned industrial processes that have practical computational time using scaled up particles.

Particle scaling rules for an elasto-plastic adhesive contact model
The Coarse Graining technique substitutes coarse grains s times larger than the original particles for the original particles with radius R. In general, a higher scaling factor s leads to a lower computation time of DEM simulations. According to [22], the scaled system should have the same energy density as the original (unscaled) system. Using the same particle density maintains the same potential and kinetic energy densities through CG [1]. The contact stiffness and damping should be also scaled precisely, to maintain the same energy losses between the scaled system and the original one.
Lommen et al. [1] establishes the CG principles for the Hertz-Mindlin contact model in both the normal and tangential spring directions. Fig. 2 shows the normal direction of the spring-damper system for both eight original particles and the substituted coarse grain with s = 2. If spherical particles are used, by maintaining the same Young's modulus in two systems, the contact stiffness and damping are identical to the equivalent stiffness and damping of the original system respectively.
To model cohesive elasto-plastic materials, we use an Adhesive Elasto-Plastic contact model developed by [32], which is implemented in the EDEM 2019 software package. This contact model is referred as EEPA (Edinburgh Elasto-Plastic Adhesive) in the software environment. Fig. 3 is a schematic diagram of the non-linear mode of the EEPA contact spring in the normal direction. The contact spring in the normal direction consists of four different parts: • Constant pull-off force f 0 (N): an ever-present adhesive force that is added to other normal forces. • Branch I, the loading spring: the contact follows the path k 1 when two particles are approaching each other. Choosing n = 1.5 makes the loading spring equivalent to the Hertz-Mindlin contact spring in the normal direction [32]. • Branch II, the unloading and reloading spring: due to plastic deformation, upon unloading the contact spring switches to the unloading and reloading spring k 2 . The specific overlap corresponding to zero force during unloading is described as the plastic overlap δ p . This plastic overlap is tracked and updated as the stress history-dependent parameter of the contact spring. The stiffness of the contact spring in Branch II, k 2 , is a function of k 1 and the plasticity ratio λ P , and it is equal to k 1 1−λ P . The plasticity ratio, λ P , controls the ratio between stiffness in branch II (k 2 ) and stiffness in branch I (k 1 ), which shows the influence of plasticity ratio at contact scale. This means by increasing the plasticity ratio, a higher level of plastic overlap occurs during contact. • Branch III, the adhesive spring: if unloading continues beyond δ p , an adhesive (negative) force is created that is limited to the maximum adhesive force f min . Once this point is reached, the adhesive spring is activated, whose stiffness is k adh . If unloading continues, two particles separate if at δ N = 0 the absolute normal force is larger than the constant pull-off force f 0 .
Eq. (1) shows the mathematical formulation of the sum of hysteretic force in the normal spring: If n = 1.5 is used, the stiffness of the contact spring in Branch I, k 1 , is calculated as follows: where E ⁎ and R ⁎ are the equivalent Young's modulus and radius of two particles that are in contact. By using λ P = 0, the contact spring is converted to an elastic spring [33]. The plastic overlap is calculated in Eq. (3): The stiffness of the adhesive spring, Branch III, is calculated as follows: where δ min is the corresponding overlap with the maximum adhesive force f min . Since the minimum overlap, δ min , is the intersection between Branch II and Branch III, then it can be considered equal . The power value for Branch III is x. Similarly to JKR theory [34], the maximum adhesive force as formulated in Eq. (5) is a function of the contact patch radius a (m) and the adhesion surface energy Δγ (J/m 2 ).
Fig . 2 shows that the equivalent stiffness of the original system k eq , which consists of s 2 pairs of series springs, can be derived using Eq. (6) for n = 1.5. Using the fact that in series springs, δ series = s δ N , and based on the definition of k 1 (Eq. (2)), k eq during loading is determined.
Maintaining the same equivalent Young's modulus for the coarse system yields: Therefore, according to Eq. (8) in Branch I, the stiffness of the coarse system is equal to the equivalent stiffness of the original system. The same scaling rule is applicable to the tangential stiffness k t [1].
Eq. (9) rewrites the normal force in unloading and reloading, Branch II, for the coarse system. First, k 2 ′ and δ P ′ are replaced by k 1 1−λ P and λ P 1 nδ N respectively. Fig. 3. The relationships of the EEPA contact spring in the normal direction [32].
Next, according to Eq. (8), k 1 ′ is replaced by ffiffi s p k 1 . Also, since the coarse system is equivalent to s 2 pairs of series springs, δ N ′ = s δ N . If the same plasticity ratio λ p is maintained for the coarse system, for n = 1.5: Therefore, according to Eq. (10), during unloading and reloading, the stiffness of the coarse system is equivalent to the original system consisting of s 2 pairs of series springs. The next step is to find scaling rules for cohesive forces. Fig. 4 compares the constant pull-off springs between the original system and the coarse system with a scaling factor of 2. If the constant pull-off force is scaled up by s 2 , according to Eq. (11) the force in the coarse system (F A ) is equal to the equivalent force in the original system (F B ). The equivalent constant pull-off forces in the original and coarse systems are therefore equal.
Similarly to Eq. (6), the equivalent stiffness of the original system in Branch III, k eq.III , is derived in Eq. (12): If the Surface Energy Δγ is scaled by the factor s for the coarse system: Furthermore, since δ N ′ = s δ N, the contact radius a is proportional to the particle radius, which is scaled up by the scaling factor s. The minimum attractive force in the coarse system f min ′ is therefore: Thus, according to Eqs. (4), (10), (11), (12), (14) and (15), in Branch III, for n = 1.5, the stiffness of the coarse system is proven to be equal to the equivalent stiffness of the original system. It is remarkable that changing the value of x, the power value of Branch III, leaves the conclu- using the fact that this point is the intersection between Branch II and Branch III.

Geometric scaling
Geometric dimensions are linearly scaled, which means that all the dimensions are scaled by the same factor, referred to in this article as S Box . The scaling of geometric kinematics is also by the same scaling factor. This allows to maintain a constant shear strain rate during geometric scaling.

Hybrid experimental plan for particle and geometric scaling
In this section, we design an experimental plan to firstly investigate the influence of the proposed particle scaling rules on bulk responses in quasi-static and dynamic regimes. Secondly, to evaluate the decoupled effect of scaling from both a particle and a geometric perspective, we develop a hybrid particle-geometric experimental plan. To isolate the effects of S Box from S P , only one of them is varied at a time. In other words, when geometry dimensions are varied, particle properties (e.g. particle size, particle density, contact settings) are maintained constant. This allows to create intersects between two levels of S Box , as illustrated previously in Fig. 1 .
Following test cases are used in this study to study behaviour of bulk material from various aspects: (1) Uniaxial confined consolidation simulation, which captures the stiffness of bulk material under vertical consolidation stresses.
Using this case, we evaluate whether bulk stiffness is scaled properly. (2) Ledge angle of repose simulation, which is widely used to develop calibrated DEM simulations [1,18,35,36]. This test cases evaluates the performance of hybrid particle-geometric scaling when modelling the free-surface flow of cohesive bulk materials. (3) Schulze ring shear simulation, which is used to model shear flow under the effect of consolidation stresses. Shear tests are commonly used to characterize bulk materials (freeflowing and cohesive) [37], and to calibrate DEM simulations [7,[38][39][40][41][42][43].
The transition from a quasi-static to a dynamic regime can be characterized using the Inertial Number, I, which has been used in [44] to study the shear flow. The ratio between inertial forces and confining pressure can be expressed using (17).
where _ γ is the shear strain rate, R P is the particle radius, ρ P is the particle density, and P is the confining pressure. I ≤ 0.01 has been characterized as a quasi-static regime [44].
DEM input parameters Main constant DEM parameters to model particles and their interaction are listed in Table 2a. These parameters are selected based on the calibrated DEM simulation of cohesive elasto-plastic coal in ring shear test that is developed in [45]. The calibrated DEM parameters have been verified in terms of shear strength and bulk density for various levels of consolidation pressure. A spherical shape is used to model DEM particles. A normal distribution of particle size with a standard deviation of 0.1 is used. A particle density of 1350 kg/m 3 has been used in [45]. In the current study, particle density is increased to 4500 kg/m 3 that allows to simulate density level of heavier bulk solids, such as cohesive iron ore [46]. DEM input variables are listed in Table 2b. Main DEM variable in all tests is particle size. A reference particle radius of 5.5 mm has been used, corresponding to S P = 1. The rotational freedom of particles can be suppressed artificially by either introducing a rolling friction model [47] or restricting rotation of particles [1,48]. Restricting rotation of particles has been done in [45] by applying a counterbalance torque in each time-step necessary to prevent rotational movement. This leads to increasing the resistance of particle against rotational torque. Restricting rotation of particles has been successfully used to resemble realistic material behaviour [1,22,30,45,46]. A restricted rotation option as the reference value is used to consider the rotational torque between particles. The plasticity ratio is another DEM variable in this study. A reference value of 0.75 is used to enable elasto-plastic behaviour of the contact spring.

Simulation setup
A virtual uniaxial confined consolidation simulation setup is used to evaluate the Coarse Graining technique under vertical consolidation stresses. Fig. 5 shows the specifications of the simulation setup, including reference box dimensions. Particles are generated using a moving particle factory plate, which fills the box from bottom to top. This avoids compaction during the particle generation step, which might be caused by the kinetic energy of the particles. Next, particles are allowed to settle and to reach a static condition, where the ratio of the kinetic energy to the potential energy is less than 10 −6 , as defined in [21]. Afterwards, the loading (consolidating) stage starts by moving the lid plate downward at a velocity of V Lid . After the bulk material has been consolidated, the unloading stage starts by moving the lid plate upward at the same velocity. A lid velocity of 4 mm/s is used at S Box = 1, which is equivalent to an axial strain rate of 0.02 s −1 . Based on (17, by using the axial strain rate, an Inertial Number, I, smaller than 0.01 during the consolidation (up to 200 kPa) is created, and therefore, the simulation procedure creates a quasi-static regime.

Tests
Two different tests are conducted using the uniaxial simulation setup. The first test, 1.1, uses the non-adhesive non-linear Elastic mode of the contact model, which is basically similar to the Hertz-Mindlin contact model. In other words, λ P , Δγ and f 0 are set to zero in the first test with n = 1.5. Fig. 6 shows the experimental plan for Test 1.1 with the uniaxial consolidation simulation. Test 1.1 includes both the upscaling and downscaling of particle size and the upscaling of box size. The horizontal axis indicates the particle scaling factor S P , and the vertical axis indicates the box scaling factor S Box . The box dimensions and top lid velocity V Lid are upscaled linearly, by scaling factors of 2.5 and 5, relative to the reference simulation setup. For each level of S Box , five different levels of S P are investigated. This results in smallest and largest particle radii of 2.2 and 33 mm respectively in Test 1.1.
In the second test, the contact plasticity is activated by setting λ P = 0.75. Test 1.2 with uniaxial consolidation simulation uses the non-linear Elasto-Plastic mode of the EEPA contact model. S P is varied from 2 to 6 while S Box is kept constant at 5.

Objectives
The tests with the uniaxial consolidation setup determine two main responses. The first is the average porosity of bulk material, as formulated in [49]. This parameter compares the packing densities after filling the box in different simulations. The second response is the pressure on the moving lid σ Lid during its displacement in the z direction. This enables the stiffness of bulk material during loading and unloading to be compared at different levels of S P and S Box .

Simulation setup
The angle of repose (α M ) is an important characteristic in bulk handling processes; according to [50], angle of repose results are useful to categorise flow properties. Multiple test procedures are available in literature to determine α M , some example are described in [49]. Using different test procedures, different values of angle of repose (α M ) for a same bulk material can be expected [1]. A ledge method setup or shear box is used to simulate the static angle of repose. Fig. 7 shows the reference test box dimensions. The container is 200 mm high, 200 mm long and 80 mm wide. The fixed parts are coloured black. The red part, which is the flap opening, starts moving at a velocity of V = 1 m/s to initiate particle flow by removing the lateral support. This leads to a shear strain rate of up to 0.2 s −1 , and I larger than 0.01 during the flow under the gravity force. Based on (17, the simulation procedure creates a dynamic regime that ends with a static condition once an angle of repose is formed.

Tests
Three tests are conducted using the ledge angle of repose setup. To be consistent with the previous uniaxial consolidation test experiments, a reference particle radius of 5.5 mm, equivalent to S p = 1, is used in the current experiments. Fig. 8 shows the experimental plan for Test 2.1, which verifies the developed coarse graining technique by applying a hybrid particle-geometric scaling approach, in which both particle size and box size are varied in a sequential manner.
In the second test, 2.2, particle size and level of cohesion are varied using a full factorial experiment to verify the coarse graining technique for different levels of cohesion. In Table 3, we define a relative cohesion term to distinguish between the expected levels of bulk cohesion, from a relative low angle of repose to high values. For example, at a medium level of relative cohesion for S P = 2, Δγ and f 0 are equal to 20 J/m 2 and − 0.32 N respectively. To create relatively high bulk cohesion for S P = 2, Δγ and f 0 are increased by 50% respectively, compared to the medium level. A decrease of 50% is also applied to create relatively low cohesion.
All the simulations in Test 2.1 are conducted at medium relative cohesion. Thus, Δγ and f 0 are varied with respect to the particle scaling factor S P according to Eqs. (13) and (11) respectively. The test box dimensions are also scaled by geometric scaling factors S Box of 0.5, 2.5, 5 and 10. Five levels of S P , the particle scaling factor, are tested at each level of S Box . This results in a total of 25 simulations in Test 2.1. First, in the smallest box size, S Box = 0.5, five different particle scales are simulated, S p equal to 0.2, 0.3, 0.4, 0.5 and 0.6. Next, five particle sizes are simulated in S Box = 1.0, that are equal to S p = 0.4, 0.6, 0.8, 1.0 and 1.2. S p equal to 0.4 and 0.6 are simulated in both S Box = 0.5 and 1.0. It is expected that by maintaining the particle size constant, the angle of repose can be compared under the effect of varying box dimensions. In other words, the link between S Box = 0.5 and 1.0 in the experimental plan is S p = 0.4 and 0.6, which allows to verify the adequacy of hybrid particle-geometric scaling using the ledge angle of repose simulations setup. The hybrid scaling is continued by using S Box = 2.5 and S p = 1, 1.5, 2, 2.5, 3. In S Box = 5, particle scales of 2, 3, 4, 5, and 6 are simulated. Using the largest geometry dimension, S Box = 10, particle scales of 4, 6, 8, 10 and 12 corresponding to particle radii of 22, 33, 44, 55 and 66 mm are simulated. Therefore, in Test 2.1 with the ledge angle of repose setup, adequacy of the following hypotheses are checked: a) The particle scaling rules that are established through Eqs. 6 to 15 is applicable for various particle sizes. b) The exact scaling, where both particle and geometry are scaled with a same scaling factor, is not applicable on cohesive materials when   an elasto-plastic adhesive contact model is used. In other words, the effect of particle scaling should be decoupled from the effect of geometry scaling, which can be done using a hybrid particlegeometric scaling approach.
In Test 2.2, the dimensions of the text box are kept constant at S Box = 5. The particle size is varied from S P = 1 to 6, and three levels of relative cohesion are investigated. This results in running a total of 15 simulations.
In Test 2.3, the effect of coarse graining is also evaluated for the case of enabling rotation of particles. The following variables are included in Test 2.3: particle scaling factor, coefficient of static friction, and coefficient of rolling friction. The rolling friction model follows the recommendation by Ai et al. [47] to use model C for quasi-static conditions. Following the suggestion of [51], the rolling stiffness of [52] is used and the viscous rolling damping torque is disabled. By enabling rotation of particles, their mobility increases, so higher restrictive forces (e.g. cohesive and friction) needs to be used, compared to the case when the rotation restricted option is used. The reference values of coefficients of static friction and rolling friction are chosen according to the calibrated model of wet sand in [53], which are equal to 0.7 and 0.8 respectively. For S P = 2, Δγ and f 0 equal to 100 J/m 2 and − 1.32 N are used respectively. This resulted in α M = 53°, similar to the angle of repose measured for S P = 2 in S box = 5 with the medium relative cohesion level when the restricted rotation option is used. Δγ and f 0 are varied with respect to the particle scaling factor S P according to eqs. (13) and (11) respectively. In this test, the ledge angle of repose simulation is done using S box = 5, and S p = 2, 3, 4, 5, and 6. Using a one-variable-at-a-time approach, level of μ s,p-p is varied between 0.3 and 0.9, and level of μ r,p-p is varied between 0.4 and 1. This allows to confirm the particle scaling rules are independent of levels of coefficient of static and rolling friction.

Objectives
The equilibrium of forces and stresses can be drawn for a critical failure angle α M,critical , as illustrated in Fig. 9. In an arbitrary cutting plane A s , different normal and shear stresses will act, depending on α M,critical . All normal and shear stresses at the free surface are equal to zero. Failure will occur once τ α exceeds the shear strength of the bulk material. According to the Mohr-Coulomb equation, the shear strength of bulk material τ s is often approximated by Eq. (17) [54]: where tan(φ) indicates the angle of internal friction of the bulk material. c denotes the cohesion of the bulk material: in other words, c is the shear strength of the bulk material if σ α = 0.If the box dimensions in the Test 2.1 are scaled up, the vertical stress acting on A s increases, due to the greater weight of bulk material. Although the exact location of A s is unknown, Eq. (17) suggests that increasing the normal stress σ α, decreases the contribution of c to the shear strength. A negative correlation between S Box and α M can therefore be expected. To enable the effect of box scaling on the angle of repose of cohesive materials to be evaluated, we need to ensure that the vertical pressure in the z direction (and consequently σ α ) is scaled correctly in DEM simulations. According to [55], normal stresses in vertical sections can be calculated using Eq. (18), in which the constant vertical stress σ v is assumed to act across the cross-sectional area A.
where g, ρ b , K, φ x , which denote standard gravity, bulk density, lateral stress ratio and wall friction angle respectively, are assumed to be constant [49]. U is the cross-sectional perimeter, and z is the bulk material height above the cross-section.
where the parameters with subscript 1 refer to the original box dimensions with S box = 1. σ v , bottom increases linearly when S box is increased, where the constant α is the rate of change. To check whether a unique α can be obtained during hybrid particle-geometric scaling, and also to ensure that the simulation results match the analytical solution, σ v,bottom is measured for all the simulations. The coordination number, or the average contacts per particle, is another important factor in quantifying particle packing [56]. This indicator is used here to evaluate the effect of both Coarse Graining and geometric scaling on particle packing. The average coordination number over all the particles is measured after the bulk material has reached a static condition. Since the number of contacts between geometry and particles is basically higher for simulations with smaller particles, only contacts between particles are included when calculating the average coordination number.
In addition, the angle of repose α M is measured using a computer image-analysis technique as proposed in [57]. Once a static angle of repose is created, α M is determined from the images by taking the coordinates of ten equally spaced points on the slope of the material. Then, the linear regression technique is used to fit a straight line to the data points and the angle of the line with the horizontal represents the angle of repose. Using the images, the bulk surface profiles after the creation of the angle of repose are also compared between simulations.
The four parameters mentioned, σ v,bottom , average coordination number, α M and bulk surface profile, are used to evaluate the effects of particle and box scaling on the angle of repose simulation.

Simulation setup
A DEM simulation of the ring shear test (RST) is set up based on [15,58]. The RST is commonly used to characterize the flowability of cohesive and free-flowing materials. A schematic cross-sectional view of the RST cell filled with particles is shown in Fig. 10.
The simulation setup and test procedure is similar to [45], in which the calibrated DEM parameters of cohesive material using the elastoplastic adhesive contact model is developed. In the simulation, the shear cell is first filled with particles. After they reach the static condition, the pre-shearing stage starts. In this stage, a uniform normal pressure of 20 kPa is first applied to the bulk material surface using the lid plate. Next, the shear cell starts to rotate over the centre line (CL) with a rotational velocity of ω cell . A range of rotational velocity between 5 and 20 degree/s was used. It was found that this range of velocity is enough to create a steady-state shear flow in the bulk material. Using a lower rotational velocity a longer time is required to to create the steady-state shear stress in simulation. For that reason, to ensure that the steadystate is always reached within the simulation time, a ω cell = 15 deg./s is used in all simulations with ring shear test. This leads to a shear strain rate of 0.49 s −1 , and a dynamic regime. In the second stage of the test, the shearing stage, the normal pressure is reduced to 2 kPa to measure the shear stress of the pre-consolidated bulk material. Fig. 11 shows the experimental plan designed to scale RST simulations. Three different shear cell sizes are created by scaling the reference cell by scaling factors of 2 and 5, while ω cell is kept constant during this test. The particle radius is also varied from 2.2 to 27.5 mm.

Objectives
The shear stress values are measured at both the pre-shearing and shearing stages, referred to in this paper as τ pre and τ shear respectively. The measured values of τ pre and τ shear for different combinations of S Cell and S P are used to evaluate the hybrid particle-geometric scaling in RST simulations. Fig. 12 shows the porosity values measured for different levels of S Box . The average porosity is measured after the particles reach a static condition, and before starting the loading stage. Upscaling the box dimensions decreases the porosity, due to a greater total mass of particles, which results in denser packing in general. These results show that S Box and average porosity are strongly correlated, with a correlation coefficient of −0.985. The maximum porosity at each level of S Box is measured for the smallest value of S P /S Box , thus the minimum porosity is measured for the highest value of S P /S Box respectively. Hence the minimum porosity occurs when S P = 6 at S Box = 5 and is equal to 0.385. This value is higher than the theoretical limit of minimum porosity for rigid spheres. According to [59], n min (also known as n Keppler ) for rigid spherical particles = approximately 0.22. The maximum variation in porosity due to the scaling of particle radius is 0.5%. This variation is probably caused by the particle generation method, where the new particle is placed in the simulation domain without any contact with neighbouring particles. The initial conditions are therefore adequately comparable at each level of S Box . Fig. 13 shows the outcome of Test 1.1, which illustrates the influence of both S Box and S P on the loading and unloading behaviour of bulk material. The vertical axis represents σ Lid and the horizontal axis indicates the vertical location of the lid plate Z Lid divided by the height of the box H Box . The σ Lid variation due to particle upscaling is determined by calculating the maximum difference between σ Lid in the simulation with the smallest S P and other simulations.

Uniaxial consolidation
As shown in Fig. 13a, for S Box = 1, loading is continued until Z Lid /H Box = 0.15, and a maximum σ Lid of 184 ± 6 kPa is measured for S P from 0.4 to 1.2. The minimum σ Lid , 0 kPa, during the unloading stage for S Box = 1 occurs at Z Lid /H Box = 0.0302 ± 0.0020. This value represents the residual deformation due to one complete cycle of loading (consolidating) and unloading the bulk material.
As shown in Fig. 13b, for S Box = 2.5, loading is continued until Z Lid / H Box = 0.12, and a maximum σ Lid of 165 ± 4 kPa is measured. A residual deformation of 0.0243 ± 0.0000 is measured. As shown in Fig. 13c, for S Box = 5, loading is continued until Z Lid /H Box = 0.12, and a maximum σ Lid of 151 ± 4 kPa is measured for S P from 2 to 6. A residual deformation of 0.0276 ± 0.0004 is measured.
On average, a standard deviation of 3% is measured for both maximum σ Lid and residual deformation at three levels of S Box . Test 1.1 therefore confirms that the uniaxial confined consolidation simulation using the non-linear non-Adhesive mode of the contact model is adequately insensitive to particle scaling. This mode of the contact model is equivalent to Hertz-Mindlin. Fig. 14 shows the outcome of Test 1.2, which uses the non-linear Elasto-Plastic mode of the contact model. Loading is continued until Z Lid /H Box = 0.095, and a maximum σ Lid of 120 ± 1 kPa is measured. This results in a standard deviation of less than 1% in maximum σ Lid with varying the particle radius from 11 to 33 mm. A residual deformation of 0.0749 ± 0.0006 is measured in Test 1.2, which is 2.7 times the residual deformation in the equivalent S Box in Test 1.1 when λ p = 0. The considerable difference in the value of residual deformation, as well as the difference in maximum σ Lid , is due to enabling the plasticity of the contact model, which has been discussed in [33] as well. A   Fig. 11. Experimental plan for hybrid particle-geometric scaling in Test 3 with RST. Fig. 12. Porosity in Test 1 (uniaxial consolidation test) before loading. standard deviation of less than 1% is measured for the residual deformation in Test 1.2. The above shows that, in both the Elastic and Elasto-Plastic modes of the contact model, the loading and unloading paths are adequately independent from S P . Coarse Graining as well hybrid particlegeometric scaling are therefore applicable to the uniaxial confined consolidation test. Fig. 15 shows the relationship between σ v,bottom , S Box and particle size. The vertical axis indicates σ v , bottom ; the horizontal axes in Fig. 15a and Fig. 15b show the box scaling factor S Box and S P /S Box respectively. According to the fitted linear regression in Fig. 15a, the simulation results match the analytical solution with α = 3.43 and a coefficient of determination R 2 of 0.9989. This confirms that the simulation is able to capture the effect of geometry scaling on the vertical pressure distribution properly. A negligible effect of S P /S Box on σ v,bottom is measured in Fig. 15b. The normal pressure on the cutting plane σ α is therefore scaled correctly in DEM simulations, and it is comparable at every level of S Box .

Angle of repose
As Fig. 16 shows, for the entire range of S P /S Box , the average coordination number follows a similar trend as S Box increases. Particle packing therefore depends on S Box . Although some deviations of the average coordination number, up to 10%, are captured for S Box = 0.5 and 1, the influence of S P /S Box on particle packing is significantly less than the effect of S Box . The increase in the average coordination number due to the scaling of box dimensions is caused by the strong correlation between S Box and vertical pressure in the bulk material, as demonstrated previously in Fig. 15. As illustrated in Fig. 17, the coordination number (averaged in horizontal directions) increases in the z direction. In other words, particle packing is in the loosest state near the bulk surface (z → 0), and the bulk material becomes denser as z/H Box increases. In general, increasing H Box in Test 2.1 creates denser packing at a higher level of S Box . Fig. 18 shows the angle of repose results of Test 2.1, in which S P and S Box are varied. With S Box = 0.5 and 1, an angle of repose of 90°is measured for all the particle sizes. If the box dimensions are scaled up with S Box = 2.5, the particles start to form an angle of repose smaller than 90°, which is equal to 63°in average. A standard deviation of only 1°is measured for S Box = 2.5. The average values of α M at S Box = 5 and 10 are 53°and 43°respectively; standard deviations of less than 1°are measured in these tests when varying the particle size. Considering the low standard deviation values, the particle-scaling thus successfully replicated similar angles of repose for the S P /S Box range investigated at each level of S Box . Therefore, hypothesis a (described in Section 3.2) on the adequacy of the particle scaling rules in the ledge angle of repose simulation is confirmed.
Furthermore, the performance of Exact Scaling, where both S P and S Box are scaled at the same time, is analysed in the angle of repose simulation. Fig. 19 shows the α M of the simulations where S P = S Box . Increasing the scaling factor, as Eq. (17) suggests, makes the cohesion term c less of a contributory factor to the shear strength compared to the angle of internal friction, tan(φ). This led to a negative non-linear relationship between the geometry scaling factor and α M . Exact Scaling is therefore an inadequate solution to scaling the angle of repose simulation for cohesive materials when an elasto-plastic adhesive contact model is used. This is consistent with findings of [12] in DEM modelling of a draw-down test. This confirms hypothesis b, which was described in Section 3.2.
In addition, in the hybrid experimental plan designed, at each level of S Box at least one S P intersects with a higher level of S Box . For example, S P = 1 is simulated in both S Box equal to 1 and 2.5, corresponding to S p / S Box equal to 1 and 0.4 respectively. According to Fig. 18, for S P = 1 at S Box = 1, an α M = 90°is captured, while at S Box = 2.5 using the same particle size and contact settings, α M is equal to 64°. The difference in the angle of repose is caused due to the differences in level of normal pressure and average coordination number, that were demonstrated earlier in this section. For that reason, following a hybrid particle-geometric scaling, first particles are scaled up from S P = 0.4 to 1.2. Second, geometry dimensions are scaled up from S Box = 1 to 2.5, by keeping the particle size and contact settings constant for S P = 1. Next, particles are scaled up from S P = 1 to higher scales i.e. S P = 3. Given the success of the particle scaling rules, α M is adequately equal between S P = 1 and 3 at S Box = 2.5. Additionally, the intersection between S Box = 1 and 2.5 is S P = 1. Therefore, up-scaling is done from S P = 0.4 to 3 in this case, by decoupling the geometry scaling from particle scaling. This is done by using an intersection point between S Box = 1 and 2.5, which is S P = 1. Applying a similar rationale, the effect of geometry scaling is decoupled from the particle scaling for the other levels of S Box. Using hybrid particle-geometric scaling, the ledge angle of repose simulation is therefore scaled up by increasing the particle size 60 times, from a particle radius of 1.1 mm to 66 mm.
Keeping S Box constant in the Test 2.2 enables the effect of the cohesion parameters Δγ and f 0 under Coarse Graining to be investigated. To analyse the results of Test 2.2, in addition to α M , the difference in bulk surface profiles is shown. Fig. 20 illustrates an example of  comparing bulk surfaces between two DEM simulation outputs. The left image shows the reference image, in which S P = 2 and relative cohesion is set at a low level. The right image compares the outcome of the simulation with S P = 5 and a similar level of relative cohesion with the reference image. The green area indicates the missing area in the second image compared with the reference image. Red lines, which represent the bulk surface offset by particle radius in the simulation with S P = 2 (2 R P,ref. = 11 mm), are used to evaluate the variation in bulk surface. The bulk surface in the coarse grained simulation therefore matches the reference simulation well. Fig. 21 compares the bulk surfaces between simulations with low (left) and high (right) relative cohesion levels. The magenta area indicates the difference due to the increase in cohesion. Both the angle of repose and the bulk surface irregularity increase as the relative cohesion increases. Fig. 22 shows comparative bulk surfaces at low, medium and high levels of relative cohesion. In general, the mismatch is zero in the upper part of piles, and is only present in the middle and lower parts of the pile. Overall, at all three levels of relative cohesion, the difference in bulk surface between coarse grained particles and the simulation with R P = 11 mm (S P = 2) is limited to 2 R P,ref = 11 mm. Fig. 23 compares α M in Test 2.2. The dashed line shows the average α M at all five levels of S P . The average α M with low relative cohesion is equal to 47°with a standard deviation of less than 1°. Increasing the relative cohesion to a medium level increases the average α M to 53°, and an average α M of 62°is measured at high relative cohesion with a standard deviation of 1°. This confirms that the Coarse Graining technique is applicable to different levels of relative cohesion.      Angle of repose results in Test 2.3, in which the rolling friction C is used, are presented in Table 4. By increasing coefficient of static friction from 0.3 to 0.9, the average α M is increased from 45°to 56°. By increasing coefficient of rolling friction from 0.4 to 1, and maintaining a constant coefficient of static friction at 0.7, average α M is increased from 47°to 54°. In all cases, by varying particle scaling factor, a standard deviation of 1°or less is captured for α M . This shows that the particle scaling rules were established in Section 2.1 are applicable when a rolling friction model (e.g. model C) is used. Additionally, scalability of particles in the ledge angle of repose simulation for elasto-plastic cohesive materials is confirmed independent of values of coefficients of static and rolling friction. Fig. 24 shows the outcome of the test with the RST simulation. For each simulation, two different shear stress values are plotted. The results of the pre-shearing stage τ pre and the shearing stage τ shear are plotted in the left and right graphs respectively. τ pre and τ shear of 20.5 kPa and 5.3 kPa respectively are measured for the smallest particle size, S P = 0.4. ±20% of the measured shear stress values for S P = 0.4 is used to evaluate the influence of particle and geometric scaling. In other words, a variation of ±20%, compared to the measured shear stress values for the smallest particle size, is considered acceptable during scaling of the ring shear test. As the left graph shows, τ pre increases by 20% when particle size is scaled up from S P = 0.4 to S P = 5 and the geometry of S Box is scaled up from 1 to 5. Similarly to the pre-shear stage, the shearing stage results show that τ shear increases by 20%, due to the scaling up of particles and geometry. In addition, a comparison of the results of the pre-shear and shearing stages shows that during the   Table 4 Effect of coefficients of static and rolling friction on the angle of repose for different particle scaling factors, S P = 2 to 6, at S Box = 5. shearing stage, a relatively higher variation in τ due to upscaling is measured. The measured values of τ shear for intersects of different levels of S Cell (i.e. S P = 1.2 and 1.8) vary considerably when geometric size is scaled up. The hybrid particle-geometric scaling approach therefore enables the upscaling of both particles and geometry in this test by using the intersect points that decouples the effect of particle and geometric scaling on shear stress .

On applying hybrid particle-geometric scaling
Using a hybrid particle-geometric scaling, the upscaling of particles and geometry can be used to develop large-scale DEM simulations of industrial granular processes, such as grabs, silo flow and transfer chutes, while a minimized computational time is maintained. To achieve this, steps shown in Table 5 are recommended to follow.
Step I is to conduct the laboratory tests to characterize complex behaviour of cohesive and elasto-plastic materials for various bulk responses, denoted by y.
Step II is to calibrate the DEM simulation replicating the laboratory tests at a scale of 1:1, which is a common calibration procedure.
Step III is to vary the geometry scale by maintaining constant particle size and contact settings. Values of bulk responses are expected to be affected by geometric scaling.
Step IV is to vary the particle scale and to compare bulk responses with outcome of step III. More scaling steps can be added to reach the desired trade-off between computational time and accuracy. Once a scaled up simulation with a reduced computational time is developed, validation should be done using in-situ experiments [1]. Validation experiments can be done in quantitative and qualitative ways; some examples can be found in [4,[10][11][12]46,[60][61][62][63][64][65].

Implementation for other contact models
The current study applies hybrid particle-geometric scaling to the Edinburgh Elasto-Plastic Adhesive (EEPA) contact model. In general, the hybrid scaling approach can be used for other cohesive DEM contact models simulating cohesive materials. For instance, the Hertz-Mindlin contact model combined with the Linear Cohesion model, as formulated in [25], can be used to model Elastic Adhesive (cohesive) bulk materials. Using the Linear Cohesion model adds an additional normal cohesive force to the Hertz-Mindlin model. The additional normal cohesive force is calculated according to Eq. (20): where k adh is the cohesion energy density (J/m 3 ) and A is the contact area (m 2 ). The particle scaling rules can be established similarly to the approach that was used earlier in Section 2.1. Using the superposition principle, the equivalent force of the original system f adh,eq , which consists of s 2 pairs of series springs, can be derived using Eq. (21). Using the fact that in series springs, δ series = s δ N , the contact area is scaled by s 2 . Maintaining k adh constant during scaling makes the additional normal cohesive force f adh scale-invariant.
The concept of the EEPA contact model is similar to the Adhesive Elasto-Plastic contact model that was earlier developed by S. Luding in [66]. The main difference between these two contact models is that in Luding's model, the contact stiffness during unloading and reloading also depends on the plastic overlap δ P . This allows the level of nonlinearity of plastic displacements during unloading and reloading to be adjusted. This difference should be taken into account if Coarse Graining is applied to the contact model described in [66].

Conclusions
This paper successfully develops a hybrid particle-geometric scaling approach that allows to scale DEM simulations by isolating the effects of varying particle size and geometric dimensions on bulk properties. Additionally, particle scaling rules were developed by extending the Coarse Graining technique proposed in [1] to incorporate two important aspects of bulk materials, their elasto-plastic behaviour and their cohesive forces.
Three different types of tests were used to confirm that the proposed particle scaling rules as well as hybrid particle-geometric scaling are applicable to quasi-static and dynamic regimes.
1. Uniaxial consolidation test at various vertical confining pressures, up to 190 kPa: the Coarse Graining technique is applicable to both the non-linear Elastic and non-linear Elasto-Plastic modes of the EEPA contact model. This was confirmed for a range of particle sizes. Table 5 Necessary steps to apply hybrid particle-geometric scaling in combination with DEM calibration.
I) Laboratory test II) Calibrated DEM simulation III) Simulating geometric scaling IV) Simulating particle scaling S Box 1 1 N1 (e.g. 5) N1 (e.g. 5) Sp 1 1 1 N1 (e.g. 5) Response 1 (e.g. bulk stiffness) y 1 y 1 y 1 ' y 1 ' Response 2 (e.g. angle of repose) y 2 y 2 y 2 ' y 2 ' Response 3 (e.g. shear strength) y 3 y 3 y 3 ' y 3 ' 2. Ledge angle of repose that investigates the shear flow of cohesive materials under gravity force: firstly, the hybrid scaling approach was successfully applied to scale up the particle size as well the geometry size. In other words, the particle size was scaled up to 60 times by isolating the effect of particle scaling from geometric scaling. Furthermore, the Exact Scaling technique, where both particle size and geometry are scaled using the same scaling factor, is inadequate for the ledge angle of repose test for cohesive (elasto-plastic) materials. Comparable initial conditions (e.g. average coordination number) cannot be created using Exact Scaling. Secondly, our particle scaling rules have been successfully applied to different levels of cohesion parameters f 0 and Δγ. There was a positive correlation between cohesion and ledge angle of repose. 3. Shear stress values in ring shear test for both the pre-shear (τ pre ) and shearing (τ shear ) stages: the variations were limited to 20% by applying the hybrid particle-geometric scaling.
We have demonstrated that the constant pull-off force f 0 and the surface energy Δγ should be scaled by the factors s 2 and s respectively during particle scaling. Furthermore, in hybrid particle-geometric scaling, only particle properties (e.g. particle size, particle interaction parameters) or geometric properties (e.g. dimensions) are varied at a time, which is the main novelty over Exact Scaling or Coarse Graining. Using a hybrid scaling, the upscaling of particles and geometry can be used to develop large-scale DEM simulations of cohesive (elasto-plastic) bulk solids with a minimized computational time.
Future research will focus on validating a scaled up DEM simulation in an industrial application where the material model is calibrated using the EEPA contact model in the three laboratory test setups that were investigated in the current study. Additionally, using the hybrid particle-geometric scaling, the bulk calibration approach can be applied in combination with the in-situ calibration approach. This ultimately enables the possibility of calibrating a DEM parameter set independent of design and type of process.

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.