Structural characteristics of ordered clusters in packs of ellipses

Structural characteristics of two-dimensional elliptic granular packs with various aspect ratios and intergranular friction coefficients were studied using the Discrete Element Method (DEM). Isotropic compaction from random unjammed state leads to a jammed state with polycrystals of orientationally ordered clusters (OOC). The OOCs were identified using a cluster labelling algorithm, based on the relative angle ∆θ between the major axes of two contacting particles. The threshold value of ∆θ was optimised to give the strongest correlation between OOCs and the force chain network. We found that the resulting OOC size distribution decays algebraically with an exponent of −2, independently of grain aspect ratio and material properties.


Introduction
The microstructural characteristics of granular materials play an important role in determining their macroscopic properties and behaviours. For simplicity, particles in granular materials are often treated as disk or sphere in two-and three-dimension, leading to considerable progress, both theoretically [1][2][3], experimentally [4,5] and numerically [4,[6][7][8]. However, this simplification makes it impossible to study a key issue: that grain shapes play a central role in, and correlates strongly with, the macroscopic structural characteristics [9][10][11][12][13]. The obvious next level in modelling grain shapes is by ellipses and ellipsoids and much of the work on such systems has focused on macroscopic nematic ordering. Yet, the details of grain-scale orientational ordering and cluster determination are also significant, as these affect strongly macroscopic features.
The aim of this paper is to investigate orientationally ordered clusters (OOCs) of two-dimensional elliptic particles on the grain-scale. We define clusters on the basis of the relative orientation of contacting grain pairs, with a threshold chosen to give the strongest correlation between OOCs and the force chain network. Then we investigated the characteristics of the distribution of the resulting OOCs.

Numerical simulations
Our numerical simulations were carried out using the discrete element method (DEM) [14]. The packing procedure of our systems is the same as in [7,8] for discs. Defining a particle's effective sizes as D = 2 √ ab, with a and b the major and minor axes of the ellipses, respectively, we generated N total = 10, 000 particles within a double * e-mail: tmasu@kz.tsukuba.ac.jp periodic domain. The particle sizes were distributed uniformly around a meanD = 1.02 and ranged between 0.663 and 1.38. The system sizes, L x = L y , ranged from 101D to 103D. This small dispersity stems from the slight total particle volume differences for each case. The initial packing fraction was set to 0.84 -just below jamming. To avoid overlaps between particles in the generation of the initial packing , we let particles move freely with a slight background viscosity until an unjammed static state was obtained. From that initial state, the random pack was compressed slowly and isotropically by changing the periodic length in both directions. The compression continued until the stress reached a desired value and the kinetic energy fell below a very small level. For more details about this procedure, see [7,8]. We studied systems of aspect ra- tios: α = a/b = 1.0 (discs), 1.5, 2.0 and 3.0, with the first used as a check against existing disc simulations. For each aspect ratio, we simulated five intergranular friction coefficients: µ = 0.01 (almost frictionless), 0.1, 0.2, 0.5 and 10.0 representing infinite friction) -altogether 20 different systems.
Following the isotropic compression, we computed for each system the probability density function of contact pairs P( f /f , ∆θ), where f /f is the contact force magnitude normalised by mean contact force, and ∆θ the relative angle between the major axes of the two interacting particles (Fig. 1). We observed that, in systems with large aspect ratios and low intergranular friction, contacts between pairs with small relative angle ∆θ are more likely to carry large contact forces and vice versa. This correlation between f and ∆θ is observed in all our systems, but it becomes weaker with decreasing α and increasing µ. Similar relation between contact forces and local grain shapes were observed numerically [15,16] and experimentally [17,18].

Cluster labelling algorithm
The cluster labelling is carried out by going over every contact and checking ∆θ. If it is smaller than a set threshold ∆θ 0 we either declare the two particles to be in the same cluster or combine the two clusters to which they belong. Expecting cluster structures to depend mainly on ∆θ 0 , we next introduce a method to evaluate the correlation between the cluster configurations and force chains. We show that this method allows us to set an optimal value for ∆θ 0 .

Correlation between clusters and force networks
Once clusters have been identified, we label contacts by Y = 1 or 0, depending on whether the pair in contact is fully inside a cluster or not, respectively. In particular, a contact between particles belonging to different clusters is labelled Y = 0. The correlation coefficient between f and Y, ρ f,Y , allows us to investigate the relation between clusters and force networks as a function of α and µ. Fig. 2 shows ρ f,Y as a function of ∆θ 0 for all the 20 systems. In the limits ∆θ 0 → 0, π/2, ρ f,Y → 0, corresponding to the extreme cases in which all contacts are either outside or inside clusters, respectively. In between, a peak appears, ρ max f,Y at some ∆θ c , corresponding to a maximum in the number of strong force chains inside the clusters. We regard ∆θ c as the optimal criterion for ∆θ because (i) the number of single particles belonging to no cluster is minimal; (ii) the fraction of strong force chains propagating within clusters is highest. The latter is expected because it is easier to transmit forces through a contact along the flatter surfaces of two ellipses and clusters contain the largest number of particles aligned in such a way.
In Fig. 3 we plot ρ max f,Y as a function of the aspect ratio α for different intergranular friction µ. While ρ max f,Y always increases with µ, the rate of this increase decreases with µ to a very weak dependence on α when µ → 10.
In Fig. 4 we show the value of ∆θ c as a function of α for different friction coefficients µ. While ∆θ c generally decreases with α, the dependence of the rate of this decrease on µ does not appear to be monotonic. This may be because the peak in Figure. 2 is not very sharp, which gives rise to an error in the measured value of ∆θ c . Nevertheless, the trend in this plot can be understood as follows: as α increases, ellipses shape approaches the long needle limit, in which case a very small value of ∆θ c is required to obtain ρ max f,Y .

Cluster size distribution
Defining the OOC size, n cl as the number of particles in the cluster, we plot in Fig. 5 the PDF of n cl , including single-particle clusters (SPC, n cl = 1). Although the range of sizes is relatively small, this PDF appears to de- creases roughly as a power law that depends weakly on α and µ. The noise at high cluster sizes is due to the reduced statistics for large clusters. In Fig. 6 we plot the ratio of single particles not belonging to any cluster, Φ = N S PC /N OOC . This ratio decays slowly but monotonically on α and clearly increases µ. This ratio is largest in systems with small α and large µ, which means that in such systems fewer particles are grouped into clusters. In Fig. 7 we plot the PDFs of the OOC sizes for all systems, excluding the single particle. Significantly, all the plots collapse onto a master curve for 2 ≤ n cl ≤ 16, which can be fitted well by a power law, P (n cl ) ∼ n A cl , with A ≈ −2. The exact fitting parameters are shown in the figure. The interpretation of this power remains to be understood. Fig. 8 shows examples of packing structures of various systems (top) and their mesoscale interpretation as a multi-phase continuum (bottom), in which the OOC clusters and the sin-  tom figure. The directions of the hatched lines represent the mean long axis orientation of the particles in the cluster, which makes these regions easy to shear in those directions. The mesoscopically-defined orientation makes this method useful as a basis for using crystal plasticity models of polycrystalline solids [19,20] to predict macroscopic granular material properties and behaviours with different grain shapes and intergranular friction coefficients.

Conclusions
To conclude, we have studied the dependence of oriented cluster structures in 2D elliptic granular assemblies under isotropic compression. We investigated the correlation, ρ f,Y , between the force chain network and clusters by varying the threshold value ∆θ 0 in our cluster labelling algorithm, and found that this correlation peaks at a certain value ∆θ c , implying that these clusters contain the highest fraction of the large force chains. Expecting such a correlation in real systems, we defined ∆θ c as the optimal   threshold for cluster identification. The value of the largest correlation coefficient, ρ max f,Y , appears to increase with α and decrease with µ. In contrast, our data suggests that ∆θ c only decreases with α. The fraction of single particles belonging to no cluster and the PDFs of cluster sizes were analysed and were found to be sensitive to both α and µ. Intriguingly, we found that P(n cl ) collapses onto a powerlaw master curve for all our systems, P(n cl ) ∼ n −2 cl . This study focused on relatively dense initial pre-compression states, with packing fractions of 0.84, and we intend to extend it to include looser pre-compression states, as well as different preparation protocols and shearing procedures.