Comparing the effects of interparticle friction coefficient and intermediate stress ratio on critical-state DEM simulations using Delaunay triangulations

Strong force chains form when any granular material is subjected to load. A prior study used Delaunay triangulations to investigate the role of interparticle friction coefficient, μ, in stabilising the strong force chains [1]. In this follow-on paper, the effects of μ and the intermediate stress ratio, b, are compared. Numerical samples were sheared triaxially until critical state was attained. The contact networks and Delaunay triangulations of the particle centroids were both obtained at the end of each simulation. As μ is increased, the numbers of contacts in the contact networks decrease consistently. The numbers of edges, faces or tetrahedra in the Delaunay triangulations all increase with increasing μ up to 0.25 and become approximately constant thereafter. Changing b has no significant effect. The percentage of faces in the triangulation comprising three contacts shows a linearly-decreasing trend with increasing angle of shearing resistance. This is because only orthogonal loads are applied. Triangular structures require larger lateral supporting forces to maintain their stability than columnar structures when subjected to an axial load; hence, is expected to be larger relative to when triangular motifs are more prevalent. An increased proportion of triangular structures therefore leads to a lower  .


Introduction
The stress distribution is highly heterogeneous within any granular material subjected to load.Strong force chains, consisting of contacting particles aligned with the major principal stress orientation, form upon loading to transmit the majority of the load and carry the deviatoric stress during shear deformation.These strong force chains are supported by a weak orthogonal contact network.A range of complementary investigation tools including DEM simulation [2], photoelasticity experiments [3] and analytical modelling [4,5] indicate that the dominant failure mechanism in a granular material is buckling failure of the strong force chains.The stability of the strong force chains increases with increasing interparticle friction coefficient, μ [6].A prior study used Delaunay triangulations to further investigate the role of μ in stabilising the strong force chains [1].In this paper, Delaunay triangulations are again used, on this occasion to compare the effects of μ and the intermediate stress ratio, b, on force chain stability.The intermediate stress ratio is defined as where ߪ ଵ ᇱ , ߪ ଶ ᇱ and ߪ ଷ ᇱ are the major, intermediate and minor principal effective stresses, respectively.

DEM simulations and analysis
A total of 20,164 polydisperse spherical particles were randomly placed within a cubical periodic cell, ensuring no overlaps between particles.Fig. 1 shows the particle size distribution of the numerical sample which is representative of Toyoura sand.A simplified Hertz-Mindlin contact model was used with a particle shear modulus of 29 GPa and Poisson's ratio of 0.12.The particle density was set at 2650 kg/m 3 .Gravity was inactive.The simulation code used was a version of LAMMPS [7] which included some modifications made by the authors.The starting point for all triaxial shearing simulations was the same: a dense, stable, frictionless isotropic assembly at a confining pressure of 100 kPa which was created by moving opposing pairs of periodic boundaries closer together under stress control.
Two different data sets were analysed: 'varying b' and 'varying µ'.The 'varying b' data set [8] consisted of eight individual simulations.Seven of these simulations used constant b values of 0.0, 0.2, 0.3, 0.4, 0.5, 0.6 or 0.8, while one plane strain simulation was run which had b varying between 0.3 and 0.4 (b = 0.3383 at critical state).Starting from the stable, isotropic sample described above, each of these 'varying b' triaxial shearing simulations was run by initially increasing µ from 0.0 to 0.25, and then using a stress-control algorithm to either impose plane strain conditions or else maintain b at a defined, constant value under a constant ߪ ଷ ᇱ condition.The 'varying µ' data set [6] consisted of six drained shearing simulations at constant µ values of 0.0, 0.1, 0.25, 0.5, 0.75 or 1.0.Taking the same stable, isotropic sample, each 'varying µ' simulation was run by changing µ from 0.0 to one of these six values before subjecting the sample to drained triaxial compression at b = 0.0.Comparing the 'varying b' and 'varying µ' data sets, one simulation with b = 0.0 and µ = 0.25 is common to both.Hence, thirteen simulations are analysed in this paper in total.Each shearing simulation continued to an axial strain of around 54%.Samples were compared only at critical state, i.e., the state at which volume and stress become constant, using data obtained at the end of each simulation.The complete macro-scale responses obtained in these simulations may be seen in [6] for five of the six samples from the 'varying µ' data set (excluding the frictionless simulation).The effect of b on the macro-scale behaviour under shearing is presented in detail in [8].
By joining the centroids of the contacting particles, the contact networks were obtained.The particle centroids were also used to compute Voronoi tessellations using Voro++ [9].Both the conventional unweighted tessellation and the radical/Laguerre tessellation which includes particle diameter in its calculation were obtained.These tessellations were converted to their Delaunay duals to enable a direct comparison between the edges of the contact networks and the edges of the associated Delaunay triangulations.As μ increases, the number of contacts decreases from around 57,000 at μ = 0.0 to 32,000 when μ = 1.0 whereas e shows the opposite trend (Fig. 2a).The number of contacts is independent of b (Fig. 2b).There is also no obvious relationship between e and b.Analysis of a more comprehensive data set has shown a slight sensitivity might be expected [8]; nonetheless, the effect of b on e is negligibly small compared to the effect of µ on e.
Results for the plane strain simulation, shown as light grey markers on Fig. 2b, are in line with those for the simulations at constant b values.

Results and discussion
The number of contacts decreases with increasing μ; however Fig. 3a shows that the number of edges in the Delaunay triangulations increases with μ up to μ = 0.25 and remains approximately constant thereafter.As for the number of contacts, the number of edges in the Delaunay triangulations is unaffected by b (Fig. 3b).Similar results are obtained for the numbers of faces or tetrahedra in the Delaunay triangulations (Fig. 4 and Fig. 5, respectively).The trends obtained for both types of triangulation are similar; the unweighted triangulations have a consistent offset from the radical equivalents which contain fewer edges (Fig. 3).These increases in the numbers of edges or faces with increasing μ (up to μ = 0.25) are implied in the literature, although usually for more idealised assemblies than were used for this paper.For mono-sized [10] or binary [11] assemblies of hard spheres, the average number of faces comprising a Delaunay tetrahedron decreases with increasing packing fraction (η).Redenbach [12] points out that this behaviour is expected as a dense packing of spheres becomes increasingly ordered as η increases.In this study, the total number of faces increases with increasing μ (Fig. 4a), with increasing e (Table 1) or equivalently with decreasing η. [13] show that the average number of edges per face is almost invariant for radical tessellations using varying ternary mixtures of spheres simulated using DEM.Hence, it may be inferred that the average number of edges per tetrahedron also decreases with increasing η.As η increases, so too does the coordination number.For an idealised 2D random tessellation, an approximate relationship is presented in [14] based on Euler's formula: increasing the coordination number, and hence η, reduces the average number of edges per tetrahedron.In 3D, a reduction in the average number of edges per tetrahedron also implies a reduction in the average number of faces per tetrahedron [14].The percentage of faces in the triangulation comprising three interparticle contacts (P fc ) decreases monotonically as μ is increased [1].Fig. 6 shows that increasing μ reduces P fc but increases ߶ ௩ ᇱ .The relationship between P fc and ߶ ௩ ᇱ is almost linear.The shear resistance increases when the proportion of triangular structures is reduced.A tentative explanation was proposed in [1], that the effect of μ on strength is much more significant than the proportion of triangular structures.It is likely that the observations in Fig. 6 are because only orthogonal loads are applied in these simulations.Triangular structures require larger lateral supporting forces to maintain their stability than columnar structures when subjected to an axial load; hence, ߪ ଷ ᇱ is expected to be larger relative to ߪ ଵ ᇱ when triangular motifs are more prevalent.Since the angle of shearing resistance is calculated as sin(߶ ௩ ᇱ ) = (ߪ ଵ ᇱ − ߪ ଷ ᇱ )/(ߪ ଵ ᇱ + ߪ ଷ ᇱ ), an increased proportion of triangular structures leads to a lower ߶ ௩ ᇱ .This contribution to strength due to microscale particle geometry reinforces the contribution to strength due to interparticle friction [6] as μ is increased.Six polydisperse numerical samples with different interparticle friction coefficients, µ, were sheared triaxially at b = 0.0 until critical state was attained for each.Eight similar samples were sheared at a constant µ of 0.25 but differing intermediate stress ratios from 0.0 to 0.8.At the end of each simulation, the particle centroids were used to obtain the contact networks and compute Delaunay tessellations.For all measures, the intermediate stress ratio had effectively no influence on behaviour, as shown in [8] for the coordination number.
As μ is increased, the numbers of edges in the contact networks decrease consistently whereas the numbers of edges, faces or tetrahedra in the Delaunay triangulations all increase with increasing μ up to μ = 0.25 and become approximately constant thereafter.The percentage of faces in the triangulation comprising three contacts (P fc ) decreases monotonically as μ is increased.The relationship between P fc and ߶ ௩ ᇱ is almost linear.The shear resistance increases when the proportion of triangular structures is reduced.A partial explanation was proposed in [1] based on the effect of μ on the rolling/sliding behaviour.Another reason is the contribution to strength due to micro-scale particle geometry.Columnar structures possess greater selfstability than triangular structures when subjected to orthogonal loads.Triangular structures require larger lateral supporting forces to maintain their stability than columnar structures when subjected to an axial load; hence, ߪ ଷ ᇱ is expected to be larger relative to ߪ ଵ ᇱ when triangular motifs are more prevalent.The definition of ߶ ௩ ᇱ thus indicates that an increased proportion of triangular structures leads to a lower ߶ ௩ ᇱ .This contribution to strength due to micro-scale particle geometry reinforces the contribution to strength due to interparticle friction [6] as μ is increased.

Fig. 1 .
Fig. 1.Particle size distribution of numerical samples compared with that of a real Toyoura sand (adapted from [1]).
Table1summarises the key macro-scale results at the critical state of each simulation.e is the void ratio (the volume of voids divided by the volume of solids), ߶ ௩ ᇱ is the angle of shearing resistance, p' = (ߪ ଵ ᇱ + ߪ ଶ ᇱ + ߪ ଷ ᇱ )/3 is the mean effective stress and ‫ݍ‬ = ඥ3 • ‫ܬ‬ ଶ where J 2 is the second invariant of the stress deviatoric tensor ‫ݏ‬ ᇱ = ߪ ᇱ − ‫‬ ᇱ • ߜ and δ ij represents the Kronecker parameter[8].For triaxial compression with b = 0.0, q simplifies to the conventional definition ‫ݍ‬ = ߪ ଵ ᇱ − ߪ ଷ ᇱ .

Fig. 2 .
Fig. 2. The variations in the number of contacts and the void ratio at critical state with (a) µ; (b) b.

Fig. 3 .
Fig. 3.The number of edges at critical state in the Delaunay triangulations against (a) µ; (b) b.

Fig. 4 .
Fig. 4. The number of faces at critical state in the Delaunay triangulations against (a) µ; (b) b.

Fig. 5 .
Fig. 5.The number of tetrahedra at critical state in the Delaunay triangulations against (a) µ; (b) b.

Table 1 .
Key macro-scale data at the end of all simulations after the attainment of critical state.