Influence of various DEM shape representation methods on packing and shearing of granular assemblies

Realistic yet efficient representation of particle shape is a major challenge for the Discrete Element Method. This paper uses angle-of-repose and direct-shear test simulations to describe the influence of several shape representation methods, and their parameters, on the bulk response of granular assemblies. Three rolling resistance models, with varying coefficient of rolling friction, are considered for spherical particles. For non-spherical particles, superquadrics with varying blockiness and multi-spheres with varying bumpiness are used to model cuboids and cylinders of several aspect ratios. We present extensive quantitative results showing how the various ways used to represent shape affect the bulk response, allowing comparisons between different approaches. Simulations of angle-of-repose tests show that all three rolling friction models can model the avalanching characteristics of cube/cuboid and cylindrical particles. Simulations of direct-shear tests suggest that both the shear strength and the dilative response of the considered non-spherical particles (but not their porosity) can only be predicted by the elasto-plastic rolling resistance model. The quantitative nature of the results allows identifying values of the shape-description parameters that can be used to obtain similar results when using alternative shape representation methods.


Introduction
For granular materials, particle shape is an important factor that can highly influence the contact-level behaviour of single particles and, accordingly, affect the bulk-scale response. Non-spherical shapes lead to geometric interlocking between particles, which hinders the motion of the particles and changes the rheology of the granular assemblies.
Achieving a realistic, yet efficient, representation of particle shape is currently one of the most important challenges for the Discrete Element Method (DEM). In DEM, spheres are usually employed to describe particle shape, due to their simplicity and low computational cost. Restraining rotations through rolling resistance models has been suggested to overcome the inherent inability of the spheres to provide geometric interlocking [14,37], and thus represent indirectly the effect of particle shape.
A more realistic, but computationally more expensive, alternative is to simulate the particle shape either through rigid clumps of spheres (multi-spheres [12]) or through explicit descriptions of particle shapes (e.g. using polyhedrals, superquadrics or other shapes). While there is extensive literature on particle shape in DEM, we only consider below several studies that have compared the bulk response of various shaped particles, namely spheres with rolling resistance, multispheres, superquadrics, and polyhedrals.
Zhou et al. [36] proposed a methodology to feed DEM simulations with realistic-shaped particles of natural sand. This article is part of the Topical Collection: Multiscale analysis of particulate micro and macro processes.
They utilized X-ray micro-computed tomography to scan the sand particles. The obtained images are then passed through several image processing techniques, which are followed by the reconstruction of particle surfaces (based on spherical harmonic analysis). The real particle shapes are approximated by clumps of overlapping spheres and fed into DEM simulations. A set of direct shear tests are then performed with the generated complex shape particles and results are compared to spherical assemblies. It is shown that non-spherical particles, due to increased interlocking, yield a higher initial void ratio, peak, and residual strength.
Seelen et al. [27] developed a DEM code to simulate particles of any convex shape. Their study consists of assessing positional and orientational ordering of the generated random granular packings in a cylindrical container. It is observed that the solid fraction is highest in the core of all samples (loose packings of spheres, cuboids, cylinders, and ellipsoids), while it is gradually decreasing towards the boundaries (walls). For cuboids and cylinders, it is shown that the solid fraction smoothly decreases as the aspect ratio AR increases (for AR > 1 ), while spheres provide the lowest solid fraction. In the case of ellipsoids with aspect ratios higher than one, however, the solid fraction of the samples increases.
Zhao et al. [33] compared the response of spherical particles with rolling resistance to superquadric shape particles (i.e. superballs and ellipsoids) in triaxial compression tests. They showed that increasing shape deviation from the sphere affects the shear strength and dilative response of the granular assembly. Additionally, adjusting rolling resistance can provide a similar influence as the shape parameter on the shear strength, however, it greatly overpredicts the volume change during shearing. They also suggested that the application of rolling resistance, to mimic the behaviour of non-spherical particles, is limited to particle shapes without elongation.
Zhou et al. [34] employed the rolling resistance model proposed for disks by Iwashita and Oda [14], and compared the results with the response of non-spherical particles (triangle and square clumps). They reported that the two approaches present different localisation modes of particle rotation and shear strain at the peak state. The sample with disks shows a clear localization band, while the sample with clump particles exhibits a more uniform localization pattern. They stated that, although the consideration of rolling resistance for disks affects the macroscopic strength and dilatancy of the samples, it is not able to reproduce the rotational behaviour as exhibited by non-spherical particles.
Wensrich and Katterfeld [32] proposed relating the coefficient of rolling friction to a normalized average eccentricity of contact. They evaluated the proposed solution through conducting angle-of-repose tests and compared the results for spheres with rolling resistance, to simple clumped particles. It was seen that similar angle of repose with nonspherical particles can be reached for spheres by considering approximately half of the estimated rolling friction. They suggested that the discrepancy can be justified by considering that the artificial rolling resistance, applied through the proposed methodology, is always acting in the reverse direction of rotation, while real interlocking among particles sometimes can lead to acceleration of the rotation.
Estrada et al. [10] compared the response of discs with rolling resistance, with regular polygonal particles (in 2D).
The tests were conducted in simple shear conditions, assessing various aspects of the bulk response. They observed that the rolling friction coefficient can be successfully related to the number of sides in polygonal particles. They also reported that there exists a good agreement between the two types of particle packings in terms of shear strength, solid fraction, force and fabric anisotropy, and probability distribution of contact forces.
In a previous study [28,29], we have compared the influence of multi-sphere bumpiness and superquadric blockiness on the response of assemblies of cubical particles in several micro-and macro-scale tests. We have observed that singleparticle behaviour is a direct function of both surface and edge properties, and that shape complexity only significantly affects the shear strength, porosity and mode of motion once the packing is dense.
A comprehensive comparison between spheres, with restrained rotation, and non-spherical particles is yet missing, with the above-mentioned studies only partially clarifying this issue. We address this by testing various shape complexity factors under various stress states, namely free flowing and dense shearing regimes.
We significantly extend our previous work [29], by considering cylinders in addition to cubes, and varying the particle aspect ratio. Additionally, the present work not only discusses the influence of shape complexity factors, but also specifically focuses on the capability of spheres with rolling friction to simulate characteristics of complex-shaped particles.
The main novelty of the current study is therefore to compare shape representations across multiple aspects, considering different rolling friction models and comparing with different particle shapes and aspect ratios, shape-representation methods, and shape parameters, considering also different testing scenarios and initial densities.
To do so, we simulate angle-of-repose and direct-shear tests to assess the macro-scale influence of several shape complexity factors: rolling friction for spheres (using three models) and shape representation (multi-sphere or superquadric), shape parameter (bumpiness or blockiness), and aspect ratio for cuboids and cylinders.

Methodology
This section provides information on the material properties and the testing procedures that have been followed. The considered particles are spheres with 10 % of polydispersity and mean radius of 1 mm (diameter ranging from 1.9mm to 2.1mm and following a normal distribution with a mean diameter of 2 mm and a standard deviation of 0.1). Both cubes and cylinders are simulated with various complexities at aspect ratios (AR) between 1 and 3. Additionally, nonspherical particles are monodisperse. The diameter of the spheres is equal to the edge length of non-spheres at AR = 1.
Cubes/cuboids and cylinders are approximated by two approaches: superquadric (SQ) particles in LIGGGHTS [21,26] and multi-sphere (MS) particles in EDEM [9]. Regarding the use of two distinct DEM codes, it should be noted that the contact detection algorithm and force calculation methodology are different for multi-sphere and superquadric particles. Furthermore, for modelling spherical particles, both codes are used to consider various rolling resistance models, see Sect. 2.3. The considered rolling resistance models are already implemented in the employed codes and are widely used in the literature.
The parameters for the considered material are chosen so that the computational cost is reasonable for non-spherical particles. Table 1 shows the material properties used in the simulations. The Hertz model with viscous damping [7] is used for normal forces, the Mindlin-Deresiewicz [23] model is used for tangential forces, and no rolling friction model is used for MS and SQ particles.
Anand et al. [3] used DEM to model the discharge of bulk granular materials from a quasi-3D, wedge-shaped hopper. The coefficient of restitution (CoR) was shown to have negligible influence on the discharge rate. Kasina [19] considered four different values of CoR ( 10 −4 , 0.1, 0.6 and 1) and investigated its influence on the bulk internal friction. It is observed that the predicted shear response is insensitive to CoR, since collisions are absent during shearing. Accordingly, as suggested by Kasina [19], a low CoR value of 10 −4 is considered to enhance the settlement of the particle assembly (i.e. to accelerate the process of reaching the steady-state during the filling).
The current study focuses on the impact of various shape complexity factors (i.e. three rolling friction models, bumpiness, blockiness, and aspect ratio). Polydispersity is only considered, to avoid crystallization, for spherical particles. The size distribution is rarely discussed for non-spherical particles in the literature and is not considered in this study. This is to limit the number of parameters for non-spherical particles in the current study. Additionally, the computational cost increases for polydisperse systems, and this is more critical for multispheres, where the time-step is calculated based on diametre of the smallest sub-sphere. For a 2D set of particles, Nguyen et al. [25] presented a detailed numerical investigation of the combined effects of size and shape polydispersity, reporting that triangles and hexagons spontaneously assemble into crystalline structures. However, for square particles with mostly side-side contacts, packings tend to form column-like structures.

Superquadrics
The equation of the shape of a SQ particle in its local coordinate system is [6]: where a, b, c are the half-lengths of the particle along its principal axes, and n 1 and n 2 are parameters that control the shape. For instance, cubical ( a = b = c = d∕2 ) and prolate cuboidal ( a = b = d∕2, c = AR ⋅ a ) particles can be modelled by superquadrics with n 1 = n 2 = N > 2 , while prolate cylindrical particles ( a = b = d∕2, c = AR ⋅ a ) can be modelled by taking n 1 = N > 2, n 2 = 2 , where N controls particle blockiness.

Multi-sphere approach
Multi-spheres, which approximate the shape of particles by overlapping or touching spheres, are used as an approxima- any desired shape. The contact force between neighbouring particles is calculated from their element spheres, using sphere-sphere contact detection. Details of the algorithm and mechanical calculations can be found in [12]. Cubes, as multi-spheres, were modelled using the EDEM software with equal-radius (d/4) overlapping sub-spheres. The number of sub-spheres in each edge of the cubes varies between 2 and 4, resulting in 8, 27 and 64 total sub-spheres per particle, denoted as MS (8), MS (27) and MS(64) correspondingly (see Fig. 1). We have also considered an extension in the number of sub-spheres around a single axis to investigate the effect of the aspect ratio (AR), another shape complexity parameter. Due to an excessive increase in the number of sub-spheres, the extension of MS(64) is not considered (remains with AR = 1 ). For the angle-of-repose test, AR is increased up to 3 (for both MS (8) and MS (27) particles), while for the direct shear test, the particle extension is considered up to AR = 2 (due to the high computational time required).
Furthermore, the representation of the cylindrical particles is addressed through considering overlapping spheres, which are aligned along a single axis. The particles consist of 2, 3, 4 and 30 sub-spheres, which leads to variation of the bumpiness. For MS(2), MS(3) and MS(4), the bumpiness is changed along the circumference of the particles, while for MS(30) a top and bottom bumpiness are also added, see Fig. 1. Note that all the non-spherical particles have equal edge length (2 mm), which is equal to diametre of the simulated spherical particles.

Rolling friction models
Spherical particles have no geometrical interlocking, thus failing for example to capture the formation of stable angles of repose [1]. Consequently, additional rolling friction models (RF) are often used to represent the effects of particle shape. Therefore, when a bulk solid is represented by perfect spheres, the coefficient of rolling friction r is considered as a DEM shape parameter.
Ai et al. [1] reviewed the commonly used rolling friction models proposed in the literature and assessed them in a series of 2D small-scale tests.
The current study uses three rolling resistance models. The first model, called here RF RVD , is the Relative Velocity Dependent (RVD) rolling friction model in EDEM. It considers equivalent radius R * and unit average rotational velocity vector RVD of the particles. The RF RVD model applies to each particle i in contact with a particle j a torque where and n ij is a unit vector pointing from particle i to the point of contact.
The second model, called here RF CI , is the default rolling resistance model in EDEM, which incorporates the radius R i and unit rotational velocity vector ̂ i of an individual particle into the torque calculation According to the definition by [1], this model is a contactindependent model (type D).
The third model, called here RF EPSD , is the EPSD2 model in LIGGGHTS [21], which is a modified version of the elastic-plastic spring-dashpot model [14,16,17] (type C in [1]). In this modification, the viscous damping term is switched Non-spherical particles (cubes and cylinders), simulated by means of superquadric (top 2 rows, with increasing blockiness parameter from left to right) and multi-sphere (bottom 2 rows, with increasing bumpiness from left to right) shape representation techniques off, since it is only effective during dynamic events [32] and to reduce the number of unknown parameters (i.e. rolling viscous damping parameter, which needs to be calibrated [20]). The torque is therefore computed incrementally using only the spring term, as where r is the incremental relative rotation angle between two particles and k r is the rolling stiffness. The magnitude of the spring term is bounded by: For determining the k r value, Ai et al. [1] suggested the approach proposed by Bardet and Huang [5]: where J n is a dimensionless parameter varying between 0.25 and 0.5 [32]; they showed however that, when considering only the spring term, the model is susceptible to oscillations. To overcome this deficiency, Wensrich and Katterfeld [32], recommended the use of the following equation, initially proposed by Iwashita and Oda [14]) where, k t is the tangential stiffness. Wensrich and Katterfeld [32] pointed out that the rolling stiffness k r calculated from Eq. (9) results in a larger value than from Eq. (8). It is shown that higher k r value leads to mobilization of shear dissipation in contacts, which further hinders the rotation of particles (nullifying the need for a damping term).
As suggested by Jiang et al. [16] the relative rotation of two 3D particles in contact, can be decomposed into two components: a) about the contact normal direction, which renders a twisting resistance torque b) on the tangential contact plane, which leads to rolling resistance moment. It is shown that both components are influencing the relative motion of the particles, which results in increased dilation (in sample volume) at both triaxial and plane-strain compression tests. However, it is shown that once rolling resistance is excluded from contacts, the strength and dilatancy reduce more (compared to the case where twisting resistance is off), which suggests the importance of rolling resistance. Consequently, in this study, to minimize the number of varying parameters and to highlight the influence of the rolling resistance component, the twisting resistance is not activated. Table 1 shows the properties of the simulated particles. Figure 2 shows the simulation set-up and the dimensions used. The orifice of the funnel has a diameter of 24mm , which is almost 12 to 4 times the equivalent particle diameter. This ensures free flow of all particles.

Test set-up and methodology
Miura et al. [24] observed a decrease in angle of repose (AoR) with increasing number of particles. Wensrich and Katterfeld [32] simulated the AoR test using both clumped particles and spherical particles with rolling friction (7000 particles). Similarly, Chen et al. [8] simulated the AoR test for super-ellipsoids, with 8000 particles. We use here 6000 particles, to avoid excessive computational time.
The particles are distributed randomly in the cone and allowed to settle under gravity for 1 s simulation time. After the packing is formed, the orifice is opened, and the discharge commences. The simulation continues until a heap is formed, and then the angle of repose is estimated using the algorithm described in [29].
For the spherical assemblies with RF RVD , special care must be taken for the simulation period after formation of the heap, as particles with RF RVD are susceptible to oscillation in static packings [1]. In the angle-of-repose test, after particles are discharged and an initial heap is formed, the rotational and translational velocities tend to reduce, however the constantly generated torque prevents static equilibrium. Applying the artificial torque to the single particles, which are at rest, increases the kinetic energy of the system and this can lead to instability. RF CI has a similar mechanism as RF RVD for applying the torque, thus facing the same problem of ignoring the particle being static or in motion. Based on our observations and those in [1], this deficiency is more problematic for higher values of r and larger timesteps. Consequently, a small time-step is assigned (1% of Fig. 2 Angle-of-repose simulation set-up the Rayleigh time) for these simulations to minimize the effect of the instabilities. Furthermore, the angle of repose (AoR) is measured at several instances (after the initial heap is formed), to detect any time-dependent changes in the heap geometry. The error bars in Fig. 3 show that for RF RVD the difference in AoR is less than 1 • for r ≤ 0.2 , while this difference is less than 2 • for RF CI (between 1 s and 3 s of simulation time). For larger values of r the overall difference is less than 3 • for both models. Thus, to have consistency between the measured data, all simulations are stopped at exactly 3 s from the opening of the orifice. Figure 3 shows that the rolling friction models highly affect the heap formation characteristics of the spherical particles. For r ≤ 0.05 comparable AoR values are obtained for the three models, however by increasing r , RF CI leads to larger heap angles compared to those of RF RVD and RF EPSD (at saturation, the difference is 25%). The difference in AoR between RF RVD and RF EPSD is relatively small (especially for r ≤ 0.2 , which is the most frequently used range). The very high AoR for the packings with RF CI can be explained by the fact that this model applies the torque to every single particle always in the opposite direction of the rotation (ignoring the relative rotational velocity of the contacting particles).

Angle or repose results
Another important objective of this study is to evaluate the capability of the rolling resistance models in describing the bulk response of non-spherical particles. Accordingly, we have considered the simulation of non-spherical particles using multi-sphere and superquadric approaches for cubes/ cuboids and cylindrical particle shapes (varying bumpiness and blockiness). Figure 4 shows the AoR for cube/cuboid shape particles. For AR = 1 , as already reported in [29], increasing blockiness from 4 to 8 contributes to an increase of the AoR of 7 • . The increase of bumpiness from MS(64) to MS (27) has an insignificant effect, whereas MS(8), with higher interlocking, provides a larger AoR value. Additionally, varying the aspect ratio from 1 to 3 increases the AoR for all MS and SQ particle types (note that results for MS(64) are only computed for AR = 1 due to the requirement of high computational time). The results also suggest that, compared to the sharp increase for SQ(N4), SQ(N6) and MS (27), the AoR for particles with the highest bumpiness and blockiness are less affected by the increase of AR. Another important observation in Fig. 4 is that increasing AR (specially for AR > 2 ) the difference for the measured AoR between different particles decreases, compared to the difference at AR = 1.
To further investigate the avalanching characteristics of the non-spherical particles, the AoR test has been also carried out considering cylindrical particles with different aspect ratios (Fig. 5). The first set of tests has been conducted considering SQ particles with different blockiness levels. The results show that increasing the N parameter from 4 to 6 contributes to formation of heaps with higher AoR values (at all aspect ratios). However, cylindrical particles with further blockiness (i.e. N = 8 and N = 10 ) do not result in any further significant effect on AoR value (this holds true for all ranges of AR).
The second set of simulations, for cylindrical particles, has been done using MS particles with varying number of sub-spheres. These tests are only conducted at AR = 1.5 with particles consisting of 2, 3, 4 and 30 sub-spheres, where the former three types have only bumpiness at the circumference, whereas MS(30) has bumpiness both on its edges (top Fig. 3 Angle of repose for spherical particles with different r values Fig. 4 Angle of repose for cube/cuboid particles and bottom) and on its circumference, see Fig. 1. Results for MS particles suggest that increasing the bumpiness at the circumference (from MS(4) to MS (2)) has little influence on the AoR. Increased bumpiness at the edges (i.e. MS (30)) is also required to achieve higher AoR angles. Particles like MS(30) have significantly increased computational cost (both due to increased number of sub-spheres and the need for sub-spheres with smaller radii), whereas changing N has no effect on the computational performance of SQ particles, see [29]. Consequently, the representation of cylinders with MS approach is not extended to other aspect ratios.

Simulation of direct shear test
This section evaluates the response of granular systems in a shear cell and provides detailed information regarding the influence of different rolling resistance models on initial porosity, shear strength and volume change characteristics of spherical particles. These results are further compared with the response of MS and SQ particles with various shape complexities.
The direct shear test (commonly referred as Jenike shear in bulk solid handling problems) is widely used for measuring flow properties of particulate solids, [13,30,36]. In this test, the granular material is placed in a split cylindrical box. Then, the material is compressed by applying a constant vertical load v ( 10kPa ) to the lid section (consolidation state). Later, the top half of the cylinder (ring) is sheared at a constant translational velocity ( 2mm∕s ), see Fig. 6. The measured quantity is the force required for this movement that can be converted to average shear stress . The velocity of the lid in LIGGGHTS is controlled by a standard PID controller [2] that compares the current acting force with the predefined target value (in the vertical direction). In EDEM, the position and velocity of the lid are controlled by means of a body-force model written in the API interface of EDEM software.
The DEM time-step was chosen as Δt = 2 × 10 −6 s (5% of Rayleigh time) in all simulations.
Two types of packing were generated to assess the dependence of the results on the density state of the initial packing. In the dense packing, the particle-particle sliding friction s and rolling friction r coefficients were set to zero at the filling stage (thus all simulations were carried out based on a single initial packing). Next, before applying v , s is changed back to 0.56 and r is systematically increased for individual packing ( r = 0.01, 0.1, 0.2, 0.3, 0.5, 0.8 and 1). For these packings, the system is allowed to rest for 0.5 s of simulation time, prior to consolidation. Accordingly, this set of simulations are referred to as equivalent dense initial packing (EDIP), or simply as dense packings.
The sample preparation is carried out by air pluviation. Particles are discharged from a funnel from a certain height to achieve the desired void ratio in the system. To have identical density all over the sample, the factory which generates the particles moves upward with a constant speed. It must be noted that the friction coefficient equals zero for nonspherical particles during the filling.
A second type of packing was established, with s = 0.56 at all stages of the test. Furthermore, r is applied to the granular assemblies both during filling and shearing (approximating the inclusion of particle shape during the packing). This is done to clarify the importance of the shape parameter (i.e. rolling resistance) in imposing the arrangement of the particles and further evaluate the shearing response of the looser packings (depending on the incorporated value of r ).
For all simulations, material properties were kept identical and packings of equivalent bulk volume were generated for all types of particles (for detailed information regarding the packing density of MS and SQ cubes refer to [29]). Furthermore, Soltanbeigi et al. [29] reported that the shape factor is more influential (on bulk response) in densely packed assemblies, accordingly the cylindrical and cube/ cuboid particles are only simulated in the dense state (i.e. by having zero sliding friction during the filling stage). In the following section (4.1), the effect of including an artificial particle shape parameter on the packing density is evaluated. Successively, the shearing response of the particles is assessed in Sect. 4.3.

Initial porosity
It is known that applying any type of rolling resistance model affects particle arrangement, and thus the porosity of the granular packing (this is an inherent feature since all models reduce the movement of single particles). Accordingly, the initial porosity ( init ) of the samples was measured before consolidation, by dividing the total volume of the voids over the volume of the shear tester, as an indication of how densely the particles are packed in the system. As mentioned before, the packings are prepared in dense and loose states (by turning on or off s and r ).
There are many studies in the literature for determining the packing characteristics of spherical particles with various properties. For the spherical particles (using DEM), Jerier et al. [15] varied the particle radius (ratio of the largest radius to the smallest one is changed between 2 and 7) and determined a solid fraction within 0.58 to 0.75 (corresponding to 0.25 < init < 0.42 ). Aste [4] reported that pouring mono-sized spheres into a container, typically a solid fraction between 0.61 and 0.62 ( init = 0.38 and 0.39) is obtained (depending on the geometry of the container, speed and height of filling). It is mentioned that slightly denser packings (i.e. solid fraction of 0.63) can be obtained by gently tapping the container. Additionally, it is reported that to achieve the random close packing limit (i.e. solid fraction of 0.64), the sample must be simultaneously tapped and compressed. Figure 7 shows porosity for different values of the rolling friction coefficient, for each of the three models. The sliding friction is only changed once, from 0 to 0.56 and then it is kept constant. For dense packings, init is shown as a single point in Fig. 7 ( s = r = 0 ), which has the lowest value as init = 0.377 (this is in line with results from [4]). Additionally, it is observed that enabling the sliding friction (i.e. s = 0.56 ) affects the init (loosens the packing to init = 0.4 ). Moreover, considering higher r values, further increases the porosity in the granular system (compared to the case with no rolling resistance, the maximum difference is about 20 %). For RF RVD , init increases linearly, while its values are smaller than for the other two models (for the full range of r ). Additionally, for r < 0.2 the RF EPSD model provides the highest init values (i.e. looser packings). For r ≥ 0.2 , results of RF CI and RF EPSD are almost comparable (i.e. both models have similar influence on packing, despite their different rotation retardation mechanisms).
Overall, all three models are seen to be effective in governing the particle motion and thus influence the particle arrangement. Additionally, it is clear from the comparison of Figs. 7 and 8 that spheres always reveal larger porosities compared to non-spherical particles. Ma et al. [22] simulated a monosize pack of spheres and varied rolling friction between 0 and 0.2. It is shown that the initial porosity is increased as rolling friction approximates the particle shape effect and enhances inter-particle interlocking.
Values of init for particles with various shapes are also computed, and summarized in Fig. 8 (results of cubes with AR = 1 have already been discussed in [29]). For SQ cuboids, it is seen that increasing aspect ratio leads to slightly increased init for SQ(N6) and SQ(N8) with AR = 1.25 , while SQ(N4) is independent of the AR value,  Fig. 8a. For MS particles, the sharpest increase in init is seen for the MS(8) particle. However, MS (27) shows almost a constant value for init (a slight increase for init is seen while moving from AR = 1 to AR = 1.25 ). Note that, similar to Sect. 3, MS(64) is not simulated for AR > 1.
Furthermore, the results for cylindrical particles are plotted in Fig. 8b. A cylindrical particle with N = 2 would represent a sphere if AR = 1 , while ellipsoids are formed for AR > 1 . It is clearly seen that the initial porosity is affected once the spheres are converted to ellipsoids. In general, it can be suggested that cylinders, at any level of shape complexity, have lower porosity than spheres. The difference in porosity for spheres (with r = 0 and s = 0.56 ) in Figs. 7 and 8b can be justified by the fact that the spheres generated by the superquadric approach are mono-sized (thus have larger init ). Moreover, for SQ(N4), SQ(N6) and SQ(N8) the aspect ratio parameter is not affecting the packing porosity (at AR > 1.5 , identical init values are obtained for all N values). Information regarding particle number and porosity is summarized in Table 2.
Finally, the packing characteristics of spherical particles (in the presence/absence of rolling resistance models), are compared to those of cuboids and cylinders through assessing Figs. 7 and 8. It can be pointed out that except for MS (8), which has an extra porosity between adjacent sub-spheres, all cubes/cuboids and cylinders (in dense state) have lower porosities than spheres.

Volume change during shearing
It is also useful to monitor the vertical displacement of the lid, during the shearing of the granular assemblies, to determine the mode of volume change (i.e. contractive or dilative).
For the dense packings, as observed in Fig. 9, for all three rolling friction models contraction of the assemblies is insignificant and dilation is the dominant mode of volume change. The packings with RF RVD experience similar lid displacement and are thus almost independent of the r value (at the shearing stage).
For a certain range of rolling friction coefficients (i.e. 0.01 < r < 0.3 ) the packings with RF CI follow a similar trend compared to RF RVD . However, for r > 0.3 , the vertical displacement D L decreases slightly (i.e. the system starts to dilate to a lesser extent). It is already seen in Sect. 3 that RF CI can provide comparable AoR to those of non-spherical particles by assigning r < 0.18 . Accordingly, it can be suggested that in the case of RF CI , the value of r has to be less than 0.2, to get reasonable results.
For RF CI with r > 0.3 sliding starts to happen among most of contacts, and accordingly the packing experiences less dilation. A high value of r or s nullifies the respective degree of freedom, and therefore contributes to occurrence of alternative mode of relative motion. A similar observation is reported by Estrada et al. [11], in which they simulated shearing of a polydisperse system under simple shear conditions and systematically varied sliding and rolling friction coefficients. Their results suggested that when sliding friction is high enough, while r is constant, the dominant mode of motion is contact rolling (the reverse case is also valid, where sliding occurs for major portion of the mobilized contacts).
For packings with RF EPSD , it can be seen that increasing rolling resistance has a more significant impact on the increase of the lid's vertical movement (i.e. D L , in this case, can be 2 times larger than that of for RF RVD and RF CI ).
For relatively loose packings (with enabled s and r for both filling and shearing), RF RVD shows a contractive Fig. 8 Initial porosity of the packings with a cube/cuboid b cylindrical particles behaviour (i.e. the samples follow the characteristics of a loose assembly) for almost all applied r , see Fig. 10a. It is only the packing with a very small value of rolling friction coefficient (i.e. r = 0.01 ) that undergoes a volume reduction  . 9 Vertical displacement of the lid during shearing for dense packings a RF RVD b RF CI c RF EPSD ▸ for the first part of shearing and then starts to present a volume expansion (i.e. behaves approximately similar to a medium dense sample). Consequently, it can be said that although RF RVD contributes to the formation of looser packings, it is unable to provide particles with dilative behaviour during the shearing process once r > 0.01. Moreover, it can be seen that RF CI is approximately following RF RVD (except for higher volume expansion when 0.01 ≤ r ≤ 0.1 ). Meanwhile, packings with RF EPSD provide a different response compared to those of RF RVD and RF CI ; regardless of the initial porosity, all the samples (except r = 0.01 ) present a half-half contractive/dilative behaviour (i.e. until the half of the total shearing distance, samples are compressed and once a threshold for is reached, dilation is initiated), which represents characteristics of a mediumdense sample.

Shear strength under direct shear
The curves of shear stress as a function of shear displacement (D), for spherical particles, are shown in Figs. 11 and 12. It can be seen that for dense packings, RF RVD has the following characteristics: a) there is a slight increase in maximum shear strength once rolling resistance is considered, but it saturates at r = 0.1 (overall, the maximum shear strength is increased by 15% compared to r = 0 ); b) for the initial stiffness and the residual strength, the effect of rolling resistance is indeed insignificant.
Additionally, for r < 0.3 , packings with RF CI have a similar shearing response to that of r = 0 . However, for larger values of r a lower shear strength is observed. As mentioned earlier, for RF CI with r > 0.3 sliding starts to happen among most of contacts, and accordingly the packing experiences a strength loss.
Results for RF EPSD are markedly different to those for RF CI and RF RVD ; increase of r results in significantly higher shear strength for the dense samples and its effect saturates in the range r = 0.5 … 0.8 . Additionally, the largest peak shear strength is almost two times the maximum shear strength for RF CI and RF RVD . For dense packings, therefore, only RF EPSD has a significant influence on the shear strength of the granular assembly, which increases with increasing rolling friction value. Figure 7 shows that rolling resistance contributes to the formation of loose packings. However, it is not clear whether applying rolling resistance for loose packings can influence the shearing response of the granular assemblies. Accordingly, the shearing response of the relatively loose packings is summarized in Fig. 12.   Fig. 10 Vertical displacement of the lid during the shearing for loose packings a RF RVD b RF CI c RF EPSD ▸ For both RF RVD and RF CI , increasing r decreases to a large extent both the initial stiffness and the maximum shear stress for the packings; nonetheless, the shear strength at critical state is approximately comparable. This reduction in shear strength can be partially related to lower packing density, however, it seems that these two models fail to calculate appropriate opposing torque once particles are compressed. This hypothesis arises from the dependence of the torque calculation on the normal force magnitude; here, due to the presence of higher normal forces in a shear cell, the calculated torque is over-predicted for these models, which forces the particle to mobilize by sliding.
Some kinks are observed for RF CI and RF RVD models, at intermediate shear displacement in Fig. 12. This can be justified by the fact that the volume contraction, happening during the shearing in loose samples, can slightly contribute to an increase in shear strength. A similar behaviour is seen in results provided by Wang and Gutierrez [31], and similarly reported by Jo et al. [18].
For packings with RF EPSD it is clear that, despite the initially loose packings, the shear strength is increased by increasing the rolling friction coefficient (up to r = 0.5 ). This response is completely different from the observations for RF RVD and RF CI , which makes RF EPSD the only appropriate rolling resistance model for obtaining reasonable direct shear results (for both loose and dense packings).
A more detailed comparison is provided in Fig. 13 for the maximum friction angle of all the spherical packings. For dense packing, results from the three models are shown, however, the loose packings are only presented for RF EPSD (since both RF CI and RF RVD fail to influence the shearing response reasonably). It is seen that only RF EPSD , with both loose and dense packing, increases its shear strength with r variation.
Results provided in Sects. 4.1, 4.2 and 4.3 give an overall understanding regarding the applicability of the three rolling resistance models to simulate packing, shearing and dilative behaviour of non-spherical particles. It is observed that all three models can affect the density of packing. However, it must be noted that RF CI and RF RVD models are not suitable for simulation of the dense shearing regimes (i.e. direct shear test). This can be justified by the fact that calculated torque in RF CI and RF RVD is a direct function of normal contact force, and no limitation is applied to gauge a maximum torque value. Since in direct shear test, an external vertical force is acting through the lid, a higher normal contact force, and thus a high torque value is foreseen. A high value of r cancels the rotational degree of freedom. Therefore, for both Fig. 11 Shear stress-displacement curves considering various values of r , for dense packings a RF RVD b RF CI c RF EPSD ▸ rolling resistance models sliding starts to happen among most contacts and accordingly the packing experiences a non-physical change in strength and dilation (this is not in agreement with mechanical expectations from shape factor). Conversely, RF EPSD can be applied for estimating shear strength and dilative response of the non-spherical granular assembly.
Non-spherical particles Similar to Sect. 3, the effect of aspect ratio on cuboids and cylindrical particles (with various surface and edge properties) is investigated within the shear cell, see Fig. 14. The peak friction angle ( Φ p ) of the packings, which is obtained through dividing the maximum value of normal and shear stresses, is plotted with respect to AR. For SQ cuboids, Φ p increases sharply for SQ(N4) up to AR = 1.5 , whereas for SQ(N6) and SQ(N8) the change in Φ p is smaller and limited by AR = 1.25 , see Fig. 14a. However, MS cuboids show less dependency on the AR value (only MS (27) shows an increase for AR = 1.25).
On the other hand, cylindrical particle with N = 2 and AR = 1 (i.e. spheres) have the lowest Φ p , which is significantly increased at AR = 1.25 (i.e. converting to ellipsoid), see Fig. 14b. For AR > 1.25 results for SQ(N2) and SQ(N4) overlap. Other cylindrical particles with sharper edges only show a small change in Φ p for AR = 1.25 (no changes are seen for AR > 1.25).

Discussion
In this section results from angle-of-repose and direct-shear tests are compared for spherical and non-spherical particles. Considering AoR results for cuboidal and cylindrical particles with various AR and comparing them with those of spherical particles (accompanied by three rolling resistance models), the following can be pointed out: a) RF RVD and RF EPSD can provide a similar AoR to that of AR = 1 for cubes with 0.2 < r < 0.4 and for cylinders with 0.1 < r < 0.2 (see Fig. 15), b) for AR > 1 , to compensate a similar AoR (for both cuboid and cylinders), spherical particles with RF RVD and RF EPSD need to have a maximum value of r = s = 0.56 (it is usually suggested that rolling resistance coefficient must be smaller/equal to sliding friction coefficient) c) for RF CI , it can be seen that assigning 0.08 < r < 0.18 can provide AoR values comparable to all range of AR (for both cubes and cylinders). Consequently, it is clarified in this section that for the angle-of-repose test, the incorporation of rolling resistance models can well mimic the avalanching response of non-spherical particles.
Furthermore, the shear strength for cubes/cuboids and cylinders is compared with spherical particles. As already discussed, only RF EPSD provides an approximately monotonic increase for the peak friction angle of spherical assemblies (i.e. with increase of r ). Accordingly, the peak friction angle of the spherical packings with RF EPSD is summarized in Fig. 16. Results in Fig. 16 suggest that dense packings with RF EPSD can be calibrated (by varying r ) to obtain Φ p values similar to that obtained for non-spherical particles.

Conclusions
The current work compares the avalanching, packing, shearing and dilative response for spherical particles (in presence of three different rolling resistance models) with common non-spherical particles with aspect ratios of one and higher (cylinders and cubes/cuboids). The non-spherical shapes are simulated by multi-sphere (MS) and superquadric (SQ) shape representation methods, each of which has unique shape complexity properties (MS particles have an inherent surface pseudo-bumpiness, whereas SQ approach can well-produce sharp edges). Accordingly, the ability of the artificial shape representation techniques (i.e. rolling resistance models) to simulate a wide range of shape complexity parameters is evaluated.
For the AoR test, results suggest that tuning the rolling friction parameter r , as a shape parameter for spheres (considering all three models), can lead to the desired heap characteristics (which is comparable to the results of nonspherical particles with various bumpiness, blockiness and aspect ratios).
For the direct shear test, it is seen that including rolling resistance can provide a wide range of initial porosities init . However, it is seen that spherical packings provide always a higher init than cubes and cylinders.
Additionally, the capability of the considered rolling resistance models in providing dilative response, as it happens for non-spherical particles during shearing, is evaluated. Results suggest that RF RVD and RF CI are addressing this phenomenon only for dense packing. However, RF EPSD is well capable of providing dilation for the granular assemblies almost for the entire range of packing densities (which is in agreement with results reported in [29], where it is shown that even loose samples of non-spherical particles tend to dilate).
For the shear strength, once r is varied for packings with various init , only RF EPSD provides Φ p that is well matching with those obtained for cubical/cuboidal and cylindrical shape particles. The issue associated with RF CDT and RF CI , in case of dilative response and shear strength, has been explored further. It is proposed that since the calculated torque is dependent on normal force magnitude, which is high in case of shear test, the calculated torque is over-predicted. Restricting rotational freedom of the spheres results in pure sliding of particles (this deficiency is less prominent for dense packings, where due to dense arrangement of particles sliding is happening the least).
In summary, the results for RF EPSD suggested that it is possible to calibrate the rolling friction r to obtain the desired individual macro-scale response of non-spherical particles (e.g. angle of repose and shear strength).
Meanwhile, results clearly showed the defects of both RF CI for RF RVD models in a dense shearing regime.
A future study should consider the influence of polydispersity on bulk response of non-spherical particles, as well as consider a loose initial packing with different rolling friction only during the shearing (similar to that applied to dense packings).
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://creativecommons.org/licenses/by/4.0/.

Fig. 15
Graphical description of equivalent r for obtaining similar AoR values with a cubes and cylinders with AR=1 b cuboids and cylinders with AR> 1 Fig. 16 Graphical description of equivalent r for obtaining similar peak friction angle to cuboid and cylinder particles