A study on the effect of grain morphology on shear strength in granular ma- terials

The Discrete Element Method (DEM) has been successfully used to further understand GM behaviour where experimental means are not possible or limited. However, the vast majority of DEM publications use simplified spheres with rolling friction to account for particle shape, with a few using clumped spheres and super quadratics to better capture grain geometric detail. In this study, we compare the shear strength of packed polyhedral assemblies to spheres with rolling resistance to account for shape. Spheres were found to have the highest shear resistance as the limited rolling friction model could not capture the geometric of rotation grains which caused reordering and dilation. This geometric arrangement causes polyhedra to align faces in the shear direction, reducing the resistance to motion. Conversely, geometric interlocking can cause jamming resulting in a dramatic increase in shear resistance. Particle aspect ratio (elongation and fatness) was found to significantly lower shear resistance, while more uniform aspect ratio’s increased shear resistance with shape non-convexity showing extremes of massive slip or jamming. Thus, while spheres with rolling friction may yield bulk shear strength similar to some polyhedra with a mild aspect ratio, the grain scale effect that leads to compaction and jamming from rotation and interlocking is missed. These results shed light on the complex impact that individual grain shape has on bulk behaviour and its importance.


Introduction
The granular materials used in geotechnical applications are either found in nature, such as sand, rocks, stones or created by industrial processing using quarrying, grinding and crushing. The former often have more considerable variability in shape and size, while the latter aims to create a uniform size and shape. While the uniformity of processed granular materials is a desired aspect, it is an energy-intensive task resulting in increased financial and environmental costs. Therefore, an understanding of how particle shape effects the strength of granular material is critical to enable the use of natural/recycled materials as well virgin material in applications such as railroad ballast, building foundations, and understanding natural phenomena such as landslides. Characteristics such as the strength, density, and wear of the solid material making up each grain are often well known with extensive data available in various conditions. However, the observed macro behavior of granular material due to interactions between the individual grains ( micro-scale) is far less understood. This is because the ability of experimental devices to measure forces between particles is minimal and often intrusive. Thus numerical simulation is the most viable option to gain insight into the microscopic behaviour of granular material. The Discrete Element Method (DEM) [1] ini- * e-mail: govender.nicolin@gmail.com * * e-mail: patrick.pizette@imt-lille-douai.fr tially developed for geotechnical applications is the most widely used method for the simulation of granular material. DEM simulations provide detailed metrics on particle behaviour at the grain scale that allows for the macroscopic behaviour to be manipulated or the microscopic behaviour of a given macro response to be determined.

Particle shape
Detailed grain-scale modelling comes at a high computational cost, with DEM simulations often taking weeks to months to complete making it unfeasible for most applications. Thus, to alleviate the computational cost of DEM, grain shape is most often simplified as a sphere. The use of contact models with non-physical parameters such as rolling resistance with extensive tuning may correctly capture the kinematic behaviour of materials with mild deviations from the spherical shape. However, the majority of materials used in geotechnical applications are highly non-spherical in stress dominated states. This is the ability of grains to rotate and pack in different configurations that forms the corner-stone of the practical uses of granular material in numerous applications. Thus low fidelity shape approximations significantly diminish the real-world value of the information extracted from such simulations. Figure 1(a) depicts a typical gravel particle with Figure  1 (b-c) representing linear curve shapes with the common feature of single point contact between individual grains.
The clumped sphere approach Figure 1 (d) is the most common complex shape abstraction; however, the number of sub spheres required to achieve the at faces that are typical of gravel is very high. Thus the bumpy surface artifacts is a caveat of this approach [2]. Figure 1 (e) depicts a polyhedral shape which approximates real gravel shapes to a high level of accuracy due to faces and sharp edges. Finally grains can be laser scanned and digitized as triangular meshes which is the most accurate generalized shape representation which can also be considered as a polyhedral representation of multiple tetrahedrons. Particle shape studies in general for large scale discrete element simulations have received little attention, although the industrial importance is well known [2][3][4]. Large scale industrial discrete element simulations can often only afford abstract particle shapes like spheres or multi-spheres with few particles. In the past decade, the progress of computing driven by graphical processing units (GPUs) has seen computational methods suited to parallel implementations achieve significant performance increases over traditional CPU implementations. Recently, it has been demonstrated that convex polyhedral shape abstractions can be computed efficiently and for millions of particles on GPUs [5][6][7] compared to the tens of thousands in the same time-frame on CPU computing architectures [2,8,9]. In this paper, the GPU based code Blaze-DEM developed by the author is used. While only a modest number of particles (22,000) is required for this study, the simulations take only a few hours with cases of a few hundreds of thousands running in a matter of days, making it possible to study effects in finer materials such as sand and crystalline powders.

Simulation of direct shear box tests
DEM simulations were performed to model the direct shear tests to highlight the shape effect. The system size are similar to the experimental works of [10]. Figure 2 (a) shows typical shear-box modeled, the bottom is sheared at a velocity of 0.1 mm/s while the top plate is held at a constant confining pressure P con f . Figure 2 (b) shows the DEM simulation with particles in an initial unconsolidated state after filling from a square inlet in a "rain" fashion. Each particle has a mass of 7.36 × 10 −5 kg with a friction coefficient of 0.312 determined from a static repose angle test and a restitution coefficient of 0.40. A numerical stiffness of 2E4 (linear spring dashpot) was chosen based on the bulk behavior by matching the height of the gravel bed in a uni-axial compression test up to confining pressure. The tangential model is a linear spring with incremental history having half the normal stiffness and bound by the Coulomb limit. In the case of spheres a Type B rolling model which employs a restive torque that is proportional to both the contact normal force and the relative velocity at the contact point is used with rolling friction coefficient of 0.05.  There is a qualitative match with the shearing of the colored bands of simulation in agreement to that of typical experiment from [10] given for sand particles (average size arround 4 mm).  Figure 4 (a) shows the shear stress as function of displacement for experiment [10] on sand particles (average size arround 4mm) and DEM simulations. There is a good match between both curves with a rapid increase in stress during the first 25% of motion there after reaching a plateau. The simulation reaches a slightly higher plateau with less fluctuation as result of the perfect geometry compared to the imperfect grains. There is a slight over-prediction of stress in the numerical simulation leading to the peak occurring a bit earlier however with the random nature of granular material packing such differences are well within expected deviations as found by O Sullivan et al. [11]. It can therefore be concluded that the choice of numerical parameters as well as the contact method in the Blaze-DEM software is indeed sufficient. The five polyhe-dral particle shapes considered in subsequent simulations are shown in Figure 4 Table 3 lists the geometric shapes properties for the particles used in the subsequent numerical study. All particles have the same mass/volume to ensure a point of similarity between the them with the bounding radius differing to account for this.   kPa. In the case of spheres, a high rolling friction coefficient (0,1) was chosen in this section in order to reproduce repose angle of the polyhedra. The result shows that the peak shear stress is the highest despite having the lowest solid fraction, the peak is also reached the quickest compared to the polyhedral shapes with the solid fraction reaching a plateau. The convex polyhedra on the other hand all reach the peak stress slower and then remain at a plateau before smoothly decreasing with no major slip events. There is also no correlation with the solid fraction as the Snome shape has the highest shear stress but not the highest solid fraction with the Rec shape having a very similar solid fraction to the Snome shape but half the peak shear stress. The HexP and Rec shapes which are "plate like" and in aspect ratio have similar shear stresses and solid fraction behaviors. The striking difference here is with the non-convex CrossP shape jamming after 5.8 mm of displacement with an corresponding to an increase of the solid packing fraction resulting a steep rise in the shear resistance. The evolution of solid packing fraction in Figure 6 (b) shows that generally a decrease in solid fraction indicates a slip event while a plateau results in a stable shear behavior. Conversely an increase in solid fraction as in the case of the CrossP shape indicates a higher shear resistance. The solid fraction also increases throughout the motion for the CrossP shape illustrating that there is re-arrangement as a result of in-homogeneity.  Figure 7 depicts the unbound and final states for selected shapes. The blocky (cube), elongated (Rec) and non-convex (CrossP) have the sparest packing exhibiting in-homogeneity with areas of high packing density as well as noticeable voids. The CrossP shape has the most inhomogeneity with a number of preferred orientations resulting in regions where shapes interlock as a solid while also having areas with very sparse packing. Under load the cubes orientate to form the most dense packing, while the spheres have a similar packing to their unbound state. The Snome shape on the other hand has a more homogeneous packing. Post shearing the cubes under go a fairly uniform shearing with the three middle stripes having a similar shape with the end bands have little shear resistance due to the ordered packing. All polyhedra have a similar symmetric deformation from the middle band to the end while the spheres show more shear either side of the middle band. This coupled with the longer plateau of the shear stress and no significant stick/slip indicates that the ability of polyhedra to rotate changing the geometric packing and interlocking allows for force to be dissipated in a more homogeneous manner. In order to understand the micro-structure the contact force network is plotted in Figure 8. There is a stark contrast between spheres (a) and the polyhedra with the spheres having the least dense network with shorter high load bearing chains. These shorter chains result in a transmission of force that is very short ranged occurring mainly at the moving boundary on the bottom right which is in contrast to Figure 8   The Snome shape shown in Figure 8 (c) has the same long chains like cubes but with a symmetric force dis-tribution that covers the width of the shear zone. This combination of long force chains and uniform contact gives the Snome shape the highest shear resistance. Conversely the non-convex CrossP shape Figure 8 (d) has short force chains that do not extend past the immediate contact neighbors resulting in a diminished shear response due to localized force transmission. The Rec and HexP shapes are both similar to the Snome shape with the HexP having a more ordered and uniform contact network. This difference in packing homogeneity explains the noisier force curve for the HexP and Rec shapes compared to the Snome due to micro-slip and in the case of cubes and spheres full slip events. Conversely the extremely noisy stress curve for the CrossP particle is a result of its very weak and chaotic contact network. The in-homogeneity of the contact network further asserts the notion that the solidfraction which is a single metric for the bulk does not provide a meaningful measure between relative systems with the same solid fraction.

Conclusion
The effect of physically considering grain geometry has been demonstrated to be important when assembly shear strength is of interest. Grain irregularity had a significant effect with shapes that packed such that the resulting contact angles were isotropic having the highest shear strength. Shape non-convexity resulted in a system with islands of high and low packing density resulting in a very low shear resistance with slip stick behavior due to grain reordering with jamming occurring due to shape interlocking. These results shed light on the complex effect that individual grain shape has on bulk behavior and the limitations of rolling friction with spheres.