1 Introduction
How to best pack or stack objects or particles is an age-old problem and of great interest to a large variety of applications (Börzsönyi & Stannarius Reference Börzsönyi and Stannarius2013). Space-filling packings can be found for equally sized simple objects such as squares, rectangles, equilateral triangles, cubes, etc. (all with flat faces). Packings of equally sized discs or spheres will, however, always contain voids with a maximum packing (solid) ratio of
${\sim}0.9069$
for discs and
${\sim}0.7405$
for spheres. One of the simplest extensions to the disc and the sphere are the ellipse (two-dimensional) and the ellipsoid (three-dimensional). While the densest known packing of ellipses is identical to the disc one (Tóth Reference Tóth1964), ellipsoids have been found both by experiments and simulations to pack denser than spheres. The densest known crystalline packing of ellipsoids is
${\sim}0.7707$
(Donev et al.
Reference Donev, Stillinger, Chaikin and Torquato2004b
), which is considerably denser than the corresponding sphere packing of
${\sim}0.7405$
, and found for aspect ratios equal to or greater than
$\sqrt{3}$
(prolate ellipsoids) and equal to or less than
$1/\sqrt{3}$
(oblate ellipsoids).
With increasing polydispersity (in size), the crystalline order vanishes and one finds instead an amorphous packing. Such a packing will consist of a set of configurations that can (initially) resist mechanical deformations such as shear. These sets of configurations collectively represent the jamming ‘point’ and are usually represented by their average packing fraction. The jamming point generated by isotropic compaction is referred to as random close packing (RCP). If the particles do not have any significant mutual friction the jamming point occurs at around
${\sim}0.84$
in packing fraction for discs and around
${\sim}0.64$
for spheres, both below their respective crystalline packings. For frictionless ellipses and ellipsoids this value naturally depends on the aspect ratio(s). However, these packings usually display a maximum in packing fraction as a function of the aspect ratio(s). For frictionless ellipses this occurs at a packing fraction of around
${\sim}0.895$
at an aspect ratio
$\unicode[STIX]{x1D6FC}$
of 1.4 when applying an isotropic compaction (Delaney et al.
Reference Delaney, Weaire, Hutzler and Murphy2005). The maximum packing fraction for frictionless and polydisperse ellipsoids is, to the authors knowledge, yet unknown.
RCPs can also be found for monodisperse ellipses and ellipsoids when they are either poured into a container (under gravitational forces) or isotropically compacted. The maximum random close packing fraction for monodisperse ellipses (and determined by pouring) coincides with the polydisperse maximum packing fraction of
${\sim}0.895$
and occurs at a similar aspect ratio (Guises et al.
Reference Guises, Xiang, Latham and Munjiza2009). For ellipsoids one finds maximum random close packing values of around
${\sim}0.70{-}0.74$
(determined by isotropic compaction), where the value depends on whether the ellipsoids are prolate (
$\unicode[STIX]{x1D6FC}^{max}\simeq 1.5$
), oblate (
$\unicode[STIX]{x1D6FC}^{max}\simeq 0.6$
) or general spheroids (
$\unicode[STIX]{x1D6FC}_{1}^{max}\simeq 1.25$
,
$\unicode[STIX]{x1D6FC}_{2}\simeq \sqrt{\unicode[STIX]{x1D6FC}_{1}}$
,
$\unicode[STIX]{x1D6FC}_{3}\simeq 1/\unicode[STIX]{x1D6FC}_{2}$
) (Donev et al.
Reference Donev, Cisse, Sachs, Variano, Stillinger, Connelly, Torquato and Chaikin2004a
; Man et al.
Reference Man, Donev, Stillinger, Sullivan, Russel, Heeger, Inati, Torquato and Chaikin2005). These jamming transitions are naturally shifted to looser packings, and consequently denoted random loose packings, if instead frictional particles are considered (Silbert Reference Silbert2010). For frictional discs and spheres the jamming transitions now occur at the lower packing fractions of
${\sim}0.80$
and
${\sim}0.60$
, respectively. Similar trend is found for monodisperse frictional ellipses and ellipsoids but with the exact values dependent on the aspect ratio(s) (Guises et al.
Reference Guises, Xiang, Latham and Munjiza2009; Garcia Reference Garcia2010).
Most aspherical particles can in fact be packed denser than the spherical cases (discs or spheres) within a certain ‘asphericity’-parameter space. This has been demonstrated for both amorphous packing of frictionless and polydisperse polygons or rods (Tillemans & Herrmann Reference Tillemans and Herrmann1995; Alonso-Marroquín et al. Reference Alonso-Marroquín, Luding, Herrmann and Vardoulakis2005; Azéma & Radjaï Reference Azéma and Radjaï2010; Azéma, Estrada & Radjaï Reference Azéma, Estrada and Radjaï2012; Azéma & Radjaï Reference Azéma and Radjaï2014; Azéma, Radjaï & Roux Reference Azéma, Radjaï and Roux2015) and monodisperse (quasi)crystalline packings of frictionless tetrahedra (Haji-Akbari et al. Reference Haji-Akbari, Engel, Keys, Zheng, Petschek, Palffy-Muhoray and Glotzer2009) or Platonic solids (Torquato & Jiao Reference Torquato and Jiao2009). Having flat faces, it is maybe not unexpected that one finds denser packings as the space between two particles can quite trivially be reduced by having two such flat surfaces in contact. Ellipses, ellipsoids, super-ellipsoids (Delaney & Cleary Reference Delaney and Cleary2010) and dimers (Schreck, Xu & O’Hern Reference Schreck, Xu and O’Hern2010) lack flat surfaces but can still pack denser than the spherical cases and must thus reduce the voids by other means (e.g. by increasing the directional ordering).
When instead applying shear to the above systems, the assemblies risk shear jam. This typically occurs at packing fractions lower than the corresponding RCP for same system (Bi et al. Reference Bi, Zhang, Chakraborty and Behringer2011). As shearing introduces anisotropies to the system this leads to both anisotropic force and contact networks. These anisotropies increase the amount of force bearing grains compared to the RCP at the same packing ratio. Intriguingly, this effect is absent for frictionless spheres, discs and polygons (Peyneau & Roux Reference Peyneau and Roux2008; Azéma et al. Reference Azéma, Radjaï and Roux2015), i.e. RCP and shear jamming coincide, suggesting that this effect is closely related to the frictional contacts.
Despite the increasing number of studies of packing of aspherical particles only a few have been devoted to flows, i.e. rheology, of those particles (Azéma & Radjaï Reference Azéma and Radjaï2010; Börzsönyi & Stannarius Reference Börzsönyi and Stannarius2013; Nagy et al. Reference Nagy, Claudin, Börzsönyi and Somfai2017). Here we carry out extensive numerical simulations of amorphous (polydisperse) packings using discrete element methods (DEM) of both frictionless and frictional ellipses of aspect ratios up to 3. The benefit of ellipses is that they have, compared to rods or spherocylinders, a continuously varying curvature. We generate packings by shearing, where we slowly decrease the shear rate to zero at constant confining pressure and without any gravitational effects. In this well-defined region we find shear-jammed packings. We report on the packing fractions, number of contacts per particle, fraction of non-rattlers, contact angles, effective friction coefficients, pressure-scaled normal stress differences, preferred orientations and nematic order parameters of those packings as a function of both aspect ratio and interparticle friction. The latter is varied from 0 (frictionless) to 10 (in the limit of ‘infinite’ frictional) in order to cover various types of flows (frictionless, frictional sliding or frictional rolling) and their cross-overs (DeGiuli, McElwaine & Wyart Reference Degiuli, McElwaine and Wyart2016; Trulsson, DeGiuli & Wyart Reference Trulsson, DeGiuli and Wyart2017).
2 Model and methods
2.1 Numerical simulations
We run simulations utilising a novel ellipse model. The model is a generalisation of the disc and sphere models commonly used for simulation of discrete granular matter. The ellipses are described by their long and short axes with lengths
$a$
and
$b$
respectively. The vector
$\boldsymbol{e}$
keeps track on the direction of the long axis during the course of the simulations, see figure 1 for definitions. Ellipses have an area equal to
$A=\unicode[STIX]{x03C0}ab$
and a moment of inertia of
$I_{z}=\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}_{s}ab(a^{2}+b^{2})/4$
, where
$\unicode[STIX]{x1D70C}_{s}$
is the particle density (mass per area). The interactions are dissipative and the normal force
$\boldsymbol{f}_{ij}^{n}$
between particle
$i$
and
$j$
is given by
if an overlap exists between them and is otherwise zero. Here
$k_{n}$
is the normal spring constant,
$\unicode[STIX]{x1D702}^{n}$
is a viscous damping which can be related to a coefficient of restitution
$e_{n}$
,
$\unicode[STIX]{x1D739}_{ij}^{n}$
the overlap projected on the normal vector
$\boldsymbol{n}$
and
$\boldsymbol{v}_{ij}^{n}$
the velocity difference in the normal direction. As standard for quasi-hard particles we chose
$k_{n}/P>1000$
. Above this threshold one finds no detectable effects on the results, which we also verified. The amount of energy dissipated upon a single collision is controlled by the coefficient of restitution
$e_{n}$
and is put equal to
$e_{n}=0.1$
in all of the simulations. This constrain is fulfilled by having the viscous damping
$\unicode[STIX]{x1D702}^{n}$
determined from
$k_{n}$
and
$e_{n}$
according to (Luding Reference Luding2008)
where
$m_{ij}=m_{i}m_{j}/(m_{i}+m_{j})$
is the reduced mass.
Definitions of the angles, distance vectors, etc.

The demanding task in the simulations is to estimate the overlap of two ellipses. The overlap is calculated by a difference between the centre-to-centre separation
$\boldsymbol{r}_{ij}$
and the closest separation
$\boldsymbol{r}_{ij}^{closest}$
between the ellipses with fixed orientations and along the
$\boldsymbol{r}_{ij}$
-direction (i.e. keeping the directional angles
$\unicode[STIX]{x1D713}_{i}$
,
$\unicode[STIX]{x1D713}_{j}$
and
$\unicode[STIX]{x1D713}_{ij}$
fixed) according to
$\unicode[STIX]{x1D739}_{ij}=\boldsymbol{r}_{ij}^{closest}(\unicode[STIX]{x1D713}_{i},\unicode[STIX]{x1D713}_{j},\unicode[STIX]{x1D713}_{ij})-\boldsymbol{r}_{ij}$
. The closest separation is determined by the algorithm of Zheng & Palffy-Muhoray (Reference Zheng and Palffy-Muhoray2007). Using the same algorithm we can also extract the two contact vectors
$\boldsymbol{c}_{i}$
and
$\boldsymbol{c}_{j}$
of the two touching ellipses defined as the vectors going from the centre-of-masses to the contact point as well as the normal vector
$\boldsymbol{n}$
at the contact point of these touching ellipses. Velocity difference is calculated through
where
$\boldsymbol{c}_{i}^{\bot ,pos}$
is here the orthogonal vector (positive oriented) to the vector
$\boldsymbol{c}_{i}$
with equal magnitude and
$\unicode[STIX]{x1D714}_{i}$
is an angular velocity. The tangential force is calculated in the same way as (2.1), but with different spring constant and damping. In all the simulations we use
$k_{t}/k_{n}=0.5$
and
$e_{t}=0.1$
. Rather than using a simple tangential projection of the overlap vector
$\unicode[STIX]{x1D739}_{ij}$
the tangential overlap spring
$\unicode[STIX]{x1D739}_{ij}^{t}$
is given by
where
$t$
is the time and
$t_{contact}$
is the time at which the contact is formed. If the particles are not in contact
$\unicode[STIX]{x1D739}_{ij}$
,
$\unicode[STIX]{x1D6FF}_{ij}^{n}$
and
$\unicode[STIX]{x1D739}_{ij}^{t}$
are trivially zero. This definition is also consistent with what has been used previously for discs and spheres. The particles also have microscopic friction modelled by a Coulomb friction law which puts an upper bound on the tangential force as
where
$\unicode[STIX]{x1D707}_{p}$
is the microscopic friction coefficient. Here we assume that the friction coefficient is an intrinsic property of the material and that it does not vary as a function of curvature. This of course an approximation but keeping it constant makes it easier to rationalise the results. The torque acting on particle
$i$
is given (in two dimensions) by
and unlike for discs and spheres
$\boldsymbol{c}_{i}$
and
$\boldsymbol{f}_{ij}^{n}$
do not have to be parallel to each other leading to a torque even in frictionless (i.e.
$\unicode[STIX]{x1D707}_{p}=0$
) cases. The tangential spring is furthermore projected down on the tangential line after rotations, a scheme identical to the one described by Luding (Reference Luding2008).
The components of the stress tensor can either be calculated via the virial stress formula
or the Love version (Edwards & Grinev Reference Edwards and Grinev2001)
where
$V$
is the area over which the stresses are calculated. From figure 1 one does, however, easily see that the two approaches yield the same answer for hard particles (i.e. where overlaps can be assumed to be zero). This was also verified in the simulations. The nematic order parameter was calculated according to
$S_{2}=\sqrt{\boldsymbol{Q}_{xx}^{2}+\boldsymbol{ Q}_{xy}^{2}}$
with
$\boldsymbol{Q}_{lk}$
being the components of the director tensor
We shear
${\approx}900$
particles with a polydispersity of
$\pm 50\,\%$
in the diameter and none in aspect ratio between two rough walls. The walls are composed of the same particles as the sheared material and with randomised directions. Besides having normal contact forces the particles also experience mutual friction modelled with the standard Coulomb model (see above). The onset of sliding (from a rolling motion) of these tangential forces is controlled by the microscopic friction coefficient
$\unicode[STIX]{x1D707}_{p}$
. Smaller microscopic friction coefficients mean weaker tangential forces and a higher fraction of sliding contacts. As all observables saturate for friction coefficients above
$\unicode[STIX]{x1D707}_{p}\gtrsim 4$
, values above this threshold can be regarded as having reached the ‘infinite’ frictional particle limit, where no tangential sliding occurs. The simulations are carried out at fixed normal pressure
$P$
which allows the walls to move in the normal
$y$
-direction, i.e. the cell can dilate or compress. The shear rate
$\dot{\unicode[STIX]{x1D6FE}}$
is set by imposing a constant velocity (difference) to the two walls. A snapshot of the two-dimensional simulation cell is found in figure 2. All the data are measured in a steady-state flow over a time
$t$
such that
$\dot{\unicode[STIX]{x1D6FE}}t\geqslant 9$
. No detectable finite size effects were found when running simulations with approximately twice as many particles with twice the width between the walls (see appendix A). Simulations are first carried out at high inertial numbers and then gradually lowered in order to find well-defined packings at
$I\rightarrow 0$
. For ellipses we redefined the inertial number by defining the average diameter of the grains as the average geometric mean of the lengths of the two axes, i.e.
$d=\sqrt{ab}$
where
$a$
and
$b$
are the lengths of the long and short axes of the ellipses, respectively. The aspect ratio is defined simply as the ratio between the long and short axes
$\unicode[STIX]{x1D6FC}=a/b$
. The cell dimensions were typically
$45d\times 20d$
. All quantities were measured in the centre of the cell, excluding roughly five layers (in units of
$d$
) close to each wall, in order to reduce boundary effects due to possible layering close to and induced by the walls. The simulations are performed in the units of
$d$
,
$Pd$
and
$\dot{\unicode[STIX]{x1D6FE}}^{-1}$
for distance, force and time. Error bars, corresponding to one estimated standard deviation, were estimated from ten equally sized (in time) block averages for each simulation.
Simulation snapshot of the system (here with
$\unicode[STIX]{x1D6FC}=3$
and
$\unicode[STIX]{x1D707}_{p}=0.4$
). Wall particles are coloured red and flowing particles blue. A darker blue indicates a larger number of contacts.

2.2 Granular rheology
It is by now well established that granular flows can be well described by simple constitutive relations as a function of the dimensionless inertial number
$I=\dot{\unicode[STIX]{x1D6FE}}d\sqrt{\unicode[STIX]{x1D70C}_{s}/P}$
(Midi Reference Midi2004; da Cruz et al.
Reference da Cruz, Emam, Prochnow, Roux and Chevoir2005), where
$\dot{\unicode[STIX]{x1D6FE}}$
is the shear rate,
$d$
the average diameter of the grains,
$\unicode[STIX]{x1D70C}_{s}$
the grain density and
$P=\unicode[STIX]{x1D70E}_{yy}$
the normal pressure. Both the effective friction coefficient
$\unicode[STIX]{x1D707}=\unicode[STIX]{x1D70E}_{xy}/P$
, with
$\unicode[STIX]{x1D70E}_{xy}$
being the shear stress, and the packing fraction
$\unicode[STIX]{x1D719}$
can be described by the simple equations
where
$\unicode[STIX]{x1D719}_{c}$
and
$\unicode[STIX]{x1D707}_{c}$
are the critical yield parameters at
$I\rightarrow 0$
(i.e. the point the flowing ellipses solidifies into an elastic solid),
$c_{\unicode[STIX]{x1D719}/\unicode[STIX]{x1D707}}$
are constants and
$\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D719}/\unicode[STIX]{x1D707}}$
are exponents that all depend on the exact details of the composition of the assembly (i.e. shape, microscopic friction coefficient, resilience and polydispersity). Similar expressions can also be found for the average number of contacts (per particle)
$Z$
and the pressure-rescaled normal stress difference
$N_{1}/P=(\unicode[STIX]{x1D70E}_{xx}-\unicode[STIX]{x1D70E}_{yy})/\unicode[STIX]{x1D70E}_{yy}$
, fraction of non-rattlers
$\unicode[STIX]{x1D701}$
(defined as the fraction of particles having 2 or more contacts), the average contact angle in the frame of the ellipses
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\boldsymbol{c}}\rangle$
(with the absolute angle between the contact vector
$\boldsymbol{c}$
and closest major axis
$\boldsymbol{e}$
of an ellipse), the average orientation (director) angle
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\hat{\boldsymbol{y}}}\rangle$
(with the angle between the major axis
$\boldsymbol{e}$
in the upper half-space and the
$y$
-axis and with positive angles if pointing in the direction of the upper wall velocity direction) and the nematic order parameter
$S_{2}$
along the alignment direction (which does not necessarily need to be in the direction of the flow). As standard in granular rheology the number of contacts (per particle) is only averaged over particles that are non-rattlers. Fits were performed either over the whole sampled
$I$
-interval or the lowest branch of this interval according to (2.10) for
$\unicode[STIX]{x1D719}(I)$
and (2.11) for
$\unicode[STIX]{x1D707}(I)$
and similar expressions for the other quantities. All through the main paper, simulation data are shown as symbols and corresponding fits are shown as lines. Critical parameters, such as those obtained from the best fits, are also shown as symbols in the figures.
3 Results
3.1 Typical rheological curves
Figures 3–5 show typical rheological curves for
$\unicode[STIX]{x1D719}(I)$
,
$Z(I)$
,
$\unicode[STIX]{x1D701}(I)$
,
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\boldsymbol{c}}\rangle (I)$
,
$\unicode[STIX]{x1D707}(I)$
,
$N_{1}/P(I)$
,
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\hat{\boldsymbol{y}}}\rangle (I)$
and
$S_{2}(I)$
for various aspect ratios and the three microscopic friction coefficients
$\unicode[STIX]{x1D707}_{p}=[0,0.4,10]$
representing respectively frictionless (figure 3), frictional sliding (figure 4) or frictional rolling (figure 5) flow regimes. As found previously for discs (and spheres) the quantities all show power-law behaviours at low inertial numbers. The best fits all show good agreement with the simulated data, seen as lines in the figures. The rheological curves for the ellipses (
$\unicode[STIX]{x1D6FC}>1$
) follow quite well the shapes (i.e. similar exponents and constants) of the corresponding disc case curves (
$\unicode[STIX]{x1D6FC}=1$
), except for the quantities which are trivially constant or ill defined in the disc cases (i.e.
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\boldsymbol{c}}\rangle$
,
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\hat{\boldsymbol{y}}}\rangle$
and
$S_{2}$
). One finds, just as for the disc cases, that
$\unicode[STIX]{x1D719}(I)$
,
$Z(I)$
and
$\unicode[STIX]{x1D701}(I)$
all are decreasing functions of
$I$
whereas
$\unicode[STIX]{x1D707}(I)$
is increasing. Typical ranges for the various fitted exponents (covering all
$\unicode[STIX]{x1D6FC}$
values with the same
$\unicode[STIX]{x1D707}_{p}$
) are found in table 1 for the three typical flows characterised by
$\unicode[STIX]{x1D707}_{p}$
. Note that the exact values of these exponents are hard to obtain, as the exponents are highly correlated with their corresponding constants, and should be considered with some scepticism. The exponents for
$\unicode[STIX]{x1D707}$
, which are close to 0.5 for the frictionless cases and close to 1 for frictional particles, are in good agreement with what have already been found for isotropic particles (Midi Reference Midi2004; da Cruz et al.
Reference da Cruz, Emam, Prochnow, Roux and Chevoir2005), even though the exact values scatter with
$\unicode[STIX]{x1D6FC}$
. The orientational ordering, measured via the quantities
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\hat{\boldsymbol{y}}}\rangle (I)$
and
$S_{2}(I)$
, are found for the frictional cases to be decreasing with increasing
$I$
. Surprisingly, the orientational ordering shows a clear non-monotonicity for the frictionless cases, where it first increases with
$I$
and subsequently decreases. The average contact angle
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\boldsymbol{c}}\rangle$
shows a strong
$I$
-dependence (increasing with
$I$
in the low
$I$
-branch) in the frictionless cases, while it is (quasi)constant over a large
$I$
-interval for frictional particles. The normal stress difference
$N_{1}/P$
is found to be a decreasing function of
$I$
for the frictional cases, becoming negative at high inertial numbers. Again the frictionless cases turn out to be the odd birds, with non-monotonic behaviours, where the normal stress differences first increase (for the high aspect ratios) and subsequently decrease, much as for the orientational ordering.
Rheological curves as a function of inertial number
$I$
for (a) packing fraction
$\unicode[STIX]{x1D719}$
, (b) contacts per particle
$Z$
, (c) fraction of non-rattlers
$\unicode[STIX]{x1D701}$
, (d) contact angle in frame of the ellipses
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\boldsymbol{c}}\rangle$
, (e) effective friction coefficient
$\unicode[STIX]{x1D707}$
, (f) pressure-rescaled normal stress difference
$N_{1}/P$
, (g) orientation angle
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\hat{\boldsymbol{y}}}\rangle$
and (h) nematic order parameter
$S_{2}$
for various aspect ratios
$\unicode[STIX]{x1D6FC}$
and with
$\unicode[STIX]{x1D707}_{p}=0$
(i.e. frictionless flows). Dotted lines are best fits according to (2.10), (2.11) and corresponding equations for the other quantities in the low
$I$
-branch.

Same as for figure 3 but with
$\unicode[STIX]{x1D707}_{p}=0.4$
(corresponding to frictional sliding flows).

Same as for figure 3 but with
$\unicode[STIX]{x1D707}_{p}=10$
(corresponding to frictional rolling flows).

Typical exponents,
$\unicode[STIX]{x1D6FD}_{X}$
.

a
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{c}\boldsymbol{\cdot }\boldsymbol{e}}^{o}\rangle$
is constant over the whole
$I$
interval, hence the exponent is irrelevant.
For a few of the frictional cases the
$\unicode[STIX]{x1D707}(I)$
-curves start to show weak non-monotonic trends at higher aspect ratios (
$\unicode[STIX]{x1D6FC}\gtrsim 2.5$
) and at low to intermediate inertial numbers. This implies that these systems have shear localisation/banding, in space or in time. However, such localisations were not detected (see appendix A for packing ratio, average velocity, pressure and stress profiles as well as time series of same quantities). A non-monotonicity could possibly be concealed by the increased noise in
$\unicode[STIX]{x1D707}$
for the very frictional cases (i.e. at
$\unicode[STIX]{x1D707}_{p}\geqslant 1$
). Looking closer at the data for
$\unicode[STIX]{x1D6FC}\gtrsim 2.5$
and
$\unicode[STIX]{x1D707}_{p}=10$
(see figure 5
e), one sees that the best-fit power laws are within the error bars of their corresponding simulated data. This suggests that non-monotonicity is a (statistical) artefact, even if more detailed studies are needed to confirm this.
The ‘critical’ values, in the mechanical sense, in the
$I\rightarrow 0$
limit do not, unlike the exponents and constants, scatter. These values describe well the quasistatic flow limit close to shear jamming and have, in general, the most evident impact on the rheological curves. In the following sections we will only discuss the critical values as
$I\rightarrow 0$
, i.e. the values at shear jamming and how they vary with aspect ratio and friction coefficient. Figure 6 shows values of the critical key quantities as a function of the aspect ratio
$\unicode[STIX]{x1D6FC}$
at various microscopic friction coefficients. These critical key quantities are all labelled by a ‘
$c$
’ in sub- or superscript.
Critical (i.e.
$I\rightarrow 0$
) (a) packing fractions
$\unicode[STIX]{x1D719}_{c}$
, (b) contacts per particle
$Z_{c}$
, (c) fraction of non-rattlers
$\unicode[STIX]{x1D701}_{c}$
, (d) contact angle in frame of the ellipses
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\boldsymbol{c}}\rangle _{c}$
, (e) effective friction
$\unicode[STIX]{x1D707}_{c}$
, (f) pressure-rescaled normal stress differences
$N_{1}^{c}/P$
, (g) orientation angle
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\hat{\boldsymbol{y}}}\rangle _{c}$
and (h) nematic order parameter
$S_{2}^{c}$
at yielding conditions as a function of the aspect ratio
$\unicode[STIX]{x1D6FC}$
at different friction coefficients
$\unicode[STIX]{x1D707}_{p}$
. Red dashed lines are best fits accordingly to (3.1) and (3.2). Grey shaded regions show the values within the ‘realistic’ range of microscopic friction coefficients. Black dashed curve in (a) shows the maximum packing for monodisperse ellipses.

3.2 Packing fraction and contact statistics as a function of aspect ratio at shear jamming
The shear-jammed packing fraction is seen, see figure 6(a), to have a sharp increase from
$\unicode[STIX]{x1D719}_{c}\simeq 0.845$
at the disc limit
$\unicode[STIX]{x1D6FC}=1$
to maxima at around
$\unicode[STIX]{x1D719}_{c}^{max}\simeq 0.892$
at an aspect ratio of
$\unicode[STIX]{x1D6FC}^{max}\simeq 1.4$
for all of the lowest friction coefficient curves (
$\unicode[STIX]{x1D707}_{p}=0{-}0.01$
). These maxima are then followed by quasi-plateaus as the aspect ratio increases further. All the packing fractions, including the maxima, are still below the maximum packing fraction for monodisperse hard discs or ellipses of
$\unicode[STIX]{x1D719}^{max}\simeq 0.907$
(independent of
$\unicode[STIX]{x1D6FC}$
). For slightly higher friction coefficients (
$\unicode[STIX]{x1D707}_{p}=0.01{-}0.1$
) the curves tend to have the same general behaviour but with shifts to lower packing fractions for comparable
$\unicode[STIX]{x1D6FC}$
values. Above
$\unicode[STIX]{x1D707}_{p}=0.1$
the curves start, beside decreasing in value, to show non-monotonic behaviour as a function of
$\unicode[STIX]{x1D6FC}$
with clear maxima in packing fractions at approximately
$\unicode[STIX]{x1D6FC}^{max}\simeq 1.6$
with values dependent on the exact value of the microscopic friction coefficient. It is noteworthy that
$\unicode[STIX]{x1D707}_{p}=0.4$
shows the sharpest non-monotonic behaviour, a value standardly used in the modelling of granular material. Above approximately
$\unicode[STIX]{x1D707}_{p}=1$
all the curves show again monotonic behaviour, but now instead as a decreasing function of
$\unicode[STIX]{x1D6FC}$
. The high friction curves all collapse to one curve (i.e. saturate) above
$\unicode[STIX]{x1D707}_{p}=4$
, showing that we have reached the infinite friction limit (
$\unicode[STIX]{x1D707}_{p}=\infty$
). Similar non-monotonic behaviour was previously found by Guises et al. (Reference Guises, Xiang, Latham and Munjiza2009) for intermediate friction coefficients
$\unicode[STIX]{x1D707}_{p}=0.5$
, even if in their study they used monodisperse particles and rather studied the RCP (through sedimentation) than shear jamming. Our results for the frictionless are, however, in disagreement with their results, Guises et al. (Reference Guises, Xiang, Latham and Munjiza2009) did not find any sharp increase in packing fraction at low aspect ratios but rather a weak decreasing trend in packing fraction as a function of
$\unicode[STIX]{x1D6FC}$
from
$\unicode[STIX]{x1D719}^{max}\simeq 0.88$
at
$\unicode[STIX]{x1D6FC}^{max}=1$
. This could possibly be explained by the use of monodisperse ellipses in their study even if it is striking that this study’s and their study’s frictional results show the same trend. In that perspective our results should be more closely related to the ones of Delaney et al. (Reference Delaney, Weaire, Hutzler and Murphy2005), who also studied polydisperse samples of ellipses. They found a non-monotonic behaviour for the frictionless packings with a maximum
$\unicode[STIX]{x1D719}^{max}\simeq 0.895$
at
$\unicode[STIX]{x1D6FC}^{max}\simeq 1.7$
, which is again in disagreement with our results which show no significant maximum. Their packings were, however, generated by isotropic inflation of the particles and showed no directional ordering (at least none was reported) and which is in disagreement with our packings (see later). One possibility is that denser packings can be reached once the particles order themselves directionally, which can be induced by shear. More studies are needed to understand if these discrepancies are due to differences between isotropic (RCP) and anisotropic (shear) generated jammed states and/or polydispersity. We do, however, notice that our results agree well with the results of Donev et al. (Reference Donev, Connelly, Stillinger and Torquato2007) (for monodisperse) and Schreck et al. (Reference Schreck, Xu and O’Hern2010), Shen et al. (Reference Shen, Schreck, Chakraborty, Freed and O’Hern2012) (for bidisperse), even if they only reported up to
$\unicode[STIX]{x1D6FC}=1.6$
and
$\unicode[STIX]{x1D6FC}=1.9{-}2.1$
, respectively. A similar non-monotonic trend, with a maximum at
$\unicode[STIX]{x1D6FC}\sim 1.7$
(i.e. very close to the maximum for frictional ellipse case with
$\unicode[STIX]{x1D707}_{p}=0.4$
) was previously found for polydisperse and frictional (
$\unicode[STIX]{x1D707}_{p}=0.5$
) disc-caped rectangular particles (Azéma & Radjaï Reference Azéma and Radjaï2010).
The number of contacts (per particle), see figure 6(b), follows quite well the same trend as for the packing fractions with a sharp increase and then a plateau for the low friction coefficients and non-monotonic (first increase and then decrease) behaviour for the intermediate ones. The high friction cases do however show only a weak trend of decreasing, but rather seems to have saturated close to the theoretically predicted value of
$3=d_{f}$
for static packings, where
$d_{f}$
is the degrees of freedom (two translational and one rotational), for frictional particles (based on a counting argument). While the frictionless results are in agreement with previous results (Delaney et al.
Reference Delaney, Weaire, Hutzler and Murphy2005; Donev et al.
Reference Donev, Connelly, Stillinger and Torquato2007; Guises et al.
Reference Guises, Xiang, Latham and Munjiza2009; Schreck et al.
Reference Schreck, Xu and O’Hern2010; Shen et al.
Reference Shen, Schreck, Chakraborty, Freed and O’Hern2012) showing clear hypostatic packings (
$Z_{c}<6=2d_{f}$
), this is the first time one finds that the packings seem to give isostatic packings for frictional ellipses. The previous simulations of Guises et al. show the same trend as for our results for intermediate friction coefficients, even if our values seems to be slightly lower (
$Z_{c}\simeq 3.5$
for
$\unicode[STIX]{x1D707}_{p}=0.4$
compared to Guises et al.’s
$Z_{c}\simeq 4.3$
for
$\unicode[STIX]{x1D707}_{p}=0.5$
, both for large
$\unicode[STIX]{x1D6FC}$
values). Hence, from both
$\unicode[STIX]{x1D719}_{c}$
and
$Z_{c}$
we can identify possibly three different types of regions for particles depending on their friction coefficients: frictionless, intermediate (frictional sliding) and very frictional regime (frictional rolling).
Figure 6(c) shows the fraction of non-rattling particles, i.e. particles with two or more contacts. While the rattlers are rapidly vanishing as the aspect ratio is increased for the low friction cases, the opposite is found (even though with a weak tendency) for the highest frictional coefficients. Hence, increasing the aspect ratio decreases the number of rattlers for all
$\unicode[STIX]{x1D707}_{p}$
values except for
$\unicode[STIX]{x1D707}_{p}=4$
and 10. The fraction of non-rattlers is found to be smaller for higher friction coefficients for the same aspect ratio. This is maybe not unexpected as higher
$\unicode[STIX]{x1D707}_{p}$
yields lower
$Z_{c}$
values in general and hence, the gap closes between what is considered as a non-rattler and the average number of contacts (per particle).
To get more microscopic information we also monitor the preferential contact of two touching particles in the frame of the ellipse, measuring the angle between the (closest) major axis of the ellipse
$\boldsymbol{e}$
and contact vector
$\boldsymbol{c}$
(see the Numerical simulations section for more information about the definitions). While
$45^{\circ }$
indicates equal probability of having the contact close to either axis, one sees in figure 6(d) that the contact point is increasingly concentrated to the tip of the ellipses as the aspect ratio grows. This tendency is stronger for the frictionless or almost frictionless cases. This trend is maybe not unexpected as a rotation of a ellipse would eventually be limited by another particle and the point that first touches that particle would most likely be the tip (with the longest axis). Noteworthy is also that the lowest frictional cases show a preferential orientation at the short axis for the lowest aspect ratios with seemingly discontinuous jumps down to
$45^{\circ }$
at
$\unicode[STIX]{x1D6FC}=1$
. This has recently also been observed for spherocylinders (Nagy et al.
Reference Nagy, Claudin, Börzsönyi and Somfai2017; Marschall & Teitel Reference Marschall and Teitel2018).
3.3 Critical friction coefficients and normal stress difference as a function of aspect ratio at shear jamming
Equally important but less studied is the critical effective friction coefficient
$\unicode[STIX]{x1D707}_{c}$
for these packings. This critical friction is a (collective) macroscopic friction coefficient and is different from the microscopic one, even though there exists a relation between the two. It can be related to the angle of repose
$\unicode[STIX]{x1D703}_{r}$
of these ensembles as
$\unicode[STIX]{x1D703}_{r}=\arctan (\unicode[STIX]{x1D707}_{c})$
(if gravity is applied). For frictionless particles one finds a non-monotonic behaviour, see figure 6(e), in the critical friction coefficient as a function of
$\unicode[STIX]{x1D6FC}$
with a small increase from
$\unicode[STIX]{x1D707}_{c}\simeq 0.10$
at
$\unicode[STIX]{x1D6FC}=1$
up to
$\unicode[STIX]{x1D707}_{c}^{max}\simeq 0.12$
at
$\unicode[STIX]{x1D6FC}^{max}\simeq 1.4$
and then an almost linear decrease down to the smallest critical friction
$\unicode[STIX]{x1D707}_{c}\simeq 0.6$
at
$\unicode[STIX]{x1D6FC}=3$
. This behaviour is remarkably close to what has already been found for spherocylinders in three dimensions (Nagy et al.
Reference Nagy, Claudin, Börzsönyi and Somfai2017). The maximum in critical friction
$\unicode[STIX]{x1D707}_{c}^{max}\simeq 0.12$
for frictionless ellipses correlates quite well with the onset of the plateau in packing fraction. As the friction coefficient increases so does the effective friction. We find that the non-monotonic behaviour persists up to
$\unicode[STIX]{x1D707}_{p}=0.1$
, with a maximum at aspect ratios
$\unicode[STIX]{x1D6FC}\simeq 1.4$
or lower, while above the curves show monotonic behaviour. When the maxima exist one naturally finds that they grow with increasing
$\unicode[STIX]{x1D707}_{p}$
. Hence, up to around
$\unicode[STIX]{x1D707}_{p}=0.1$
ellipses with large aspect ratios yield more easily than discs while packings composed of ellipses with aspect ratios around 1.4 yield harder/tougher. This behaviour changes drastically as
$\unicode[STIX]{x1D707}_{p}$
is increased further. Here one finds a monotonic behaviour with an increased effective friction as the aspect ratio is increased with the highest values around
$\unicode[STIX]{x1D707}_{c}\sim 0.4{-}0.5$
for the largest aspect ratios. Surprisingly, one finds that the ‘infinite’ frictional cases have lower effective frictions and hence yield more easily compared to both
$\unicode[STIX]{x1D707}_{p}=0.4$
and 1. This could possibly be linked to the reduced orientational ordering in the ‘infinite’ friction limit compared to both
$\unicode[STIX]{x1D707}_{p}=0.4$
and 1, suggesting that the particles movement is less collective and/or less constrained in the former, leading to a lower effective friction.
The pressure-scaled normal stress differences
$N_{1}^{c}/P$
, see figure 6(f), show more noise than the other observables, but still a clear trend can be picked up where
$N_{1}^{c}/P$
increases with
$\unicode[STIX]{x1D6FC}$
. Surprisingly, one finds the lowest increase when
$\unicode[STIX]{x1D707}_{p}\simeq 0.4$
which is almost zero up to an aspect ratio of 2.
3.4 Alignment in the flow
We now turn our attention to the alignment of the particles. We monitor the directions of the particles through the preferred angle of the long axis in the frame of the shear cell. Figure 6(g) shows what happens to the directional angle (taken as the angle between the major axis in the upper half-space and the
$y$
-axis) as one elongates the particles, i.e. increase
$\unicode[STIX]{x1D6FC}$
. Zero angle signifies here that the ellipses preferentially face the walls in an orthogonal way but it could also be that they have a random directional orientation. In order to distinguish between these two cases we have chosen to also monitor the nematic order parameter
$S_{2}$
, where
$S_{2}=0$
signifies a totally random orientation and
$S_{2}=1$
is one in which all the particles point in the same direction. Naturally both
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\hat{\boldsymbol{y}}}\rangle _{c}$
(subscripts are for the direction of the major axis
$\boldsymbol{e}$
of the ellipse and unit vector in the
$y$
-direction of the cell
$\hat{\boldsymbol{y}}$
.) and
$S_{2}^{c}$
start at zero at
$\unicode[STIX]{x1D6FC}=1$
as the particles are isotropic (and inherently random in direction). Both
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\hat{\boldsymbol{y}}}\rangle _{c}$
and
$S_{2}^{c}$
then increase as
$\unicode[STIX]{x1D6FC}$
is increased, showing that the packings get an increased directional order as they get elongated. The positive values show that the ellipses prefer to lay in the extension direction (also seen in figure 2). One can clearly distinguish two behaviours, one for highly frictional particles and one for all the less frictional cases. For the latter the nematic order parameter slowly approaches unity as the aspect ratio increases with an alignment that gradually increases towards
$90^{\circ }$
(i.e. alignment in the
$x$
-direction) which surpasses the extension direction of
$45^{\circ }$
at around
$\unicode[STIX]{x1D6FC}\simeq 2$
. The highly frictional cases, in contrast, never surpass an angle of
${\sim}18^{\circ }$
or reach a nematic order parameter higher than 0.4 (both reached as
$\unicode[STIX]{x1D6FC}\rightarrow \infty$
).
3.5 Scaling behaviours for ellipses with ‘infinite’ friction coefficient
We now try to rationalise the findings by extracting scaling laws for the ‘infinite’ frictional cases. We fit all the critical quantities as function of the ‘effective’ elongation of the major axis relative the disc case (
$\sqrt{\unicode[STIX]{x1D6FC}-1}$
). (Other choices as
$(\sqrt{\unicode[STIX]{x1D6FC}}-1)$
and
$\sqrt{\unicode[STIX]{x1D6FC}^{2}-1}$
work equally satisfactorily, but yield different constants and exponents.) The motivation to such a form is based on a relation between the critical parameters and an excess excluded area of
$(\unicode[STIX]{x1D6FC}-1)$
of the ellipse compared to the disc upon a revolution of the ellipse. As a first guess and inspired by the power laws of (2.10) and (2.11) the critical quantities of the ellipses are either fitted according to
or as for
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\hat{\boldsymbol{y}}}\rangle _{c}$
and
$S_{2}^{c}$
to
where
$A/B$
represent the critical parameters for the ellipses, discs, or possibly saturation values as
$\unicode[STIX]{x1D6FC}\rightarrow \infty$
. The fits require the determination of the three unknowns: the constant
$c_{X}$
, the exponent
$\unicode[STIX]{x1D6FD}_{X}$
and either the saturation value or the disc value. The best-fit parameters are presented below in table 2 and shown in figure 6 as red dashed lines, which all show good fits. Even if the data are of limited experimental use (besides showing limiting behaviours) they serve as a theoretical benchmark. From the fits one can extract two different scaling exponents, one around
${\sim}1.5$
and the other around
${\sim}3$
. One finds especially that
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{c}\boldsymbol{\cdot }\boldsymbol{e}}^{o}\rangle _{c}$
and
$\unicode[STIX]{x1D707}_{c}$
scale in a similar fashion, suggesting that indeed there is a relation between these two quantities.
3.6 Observables at shear jamming as a function of the microscopic friction coefficient
Same as figure 6 but as a function of friction coefficient
$\unicode[STIX]{x1D707}_{p}$
for various aspect ratios
$\unicode[STIX]{x1D6FC}$
between 1 and 3. Dotted lines serve as guides to the eye.

We now look at the same observables but keeping the elongation
$\unicode[STIX]{x1D6FC}$
fixed and varying the friction coefficient. The results are summarised in figure 7. While
$\unicode[STIX]{x1D719}_{c}$
,
$Z_{c}$
and
$\unicode[STIX]{x1D701}_{c}$
all show monotonically decreasing behaviours towards plateau values as a function of
$\unicode[STIX]{x1D707}_{p}$
(at fixed
$\unicode[STIX]{x1D6FC}$
) the other observables (
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{c}\boldsymbol{\cdot }\boldsymbol{e}}\rangle _{c}$
,
$\unicode[STIX]{x1D707}_{c}$
,
$N_{1}^{c}/P$
,
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\hat{\boldsymbol{y}}}\rangle _{c}$
, and
$S_{2}^{c}$
) show clearly non-monotonic behaviours with increasing values up to around
$\unicode[STIX]{x1D707}_{p}\sim 0.1{-}1$
. The exceptions being the disc and the lowest studied elongation (
$\unicode[STIX]{x1D6FC}=1.025$
). Above around
$\unicode[STIX]{x1D707}_{p}\sim 0.1{-}1$
the observables
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\boldsymbol{c}}\rangle _{c}$
,
$\unicode[STIX]{x1D707}_{c}$
,
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\hat{\boldsymbol{y}}}\rangle _{c}$
and
$S_{2}^{c}$
all decrease to eventually flatten out to the ‘infinite’ frictional limit. Hence, these observables show a clear maximum around this range (0.1–1) of microscopic friction coefficients. The odd case being
$N_{1}/P$
which has both a maximum at around
$\unicode[STIX]{x1D707}_{p}=0.1$
followed by a minimum around just
$\unicode[STIX]{x1D707}_{p}\simeq 0.4$
and with a value close to zero. This minimum is then followed by an increase to the
$\unicode[STIX]{x1D707}_{p}=\infty$
plateau value. Taking a closer look at both
$\unicode[STIX]{x1D707}_{c}$
and
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{c}\boldsymbol{\cdot }\boldsymbol{e}}\rangle _{c}$
there seems to be a correlation between these two observables for large aspect ratios with two plateaus delimited by a maximum. This was further evident from the exponents (as a function of
$\unicode[STIX]{x1D6FC}$
) of these two quantities for the infinite frictional case. Indeed, one does expect a relationship between contact angle and the stress distributions (shear and normal stresses).
The critical packing fraction
$\unicode[STIX]{x1D719}_{c}$
is monotonically decreasing as a function of
$\unicode[STIX]{x1D707}_{p}$
for all aspect ratios. However, the small aspect ratio curves and the large aspect ratio curves cross each other in the frictional sliding regime (i.e. in the region
$\unicode[STIX]{x1D707}_{p}\sim 0.1{-}1$
), again illustrating that ellipses shear jam at higher packing fractions than discs if they are frictionless but that the reverse is true if they are very frictional. Similar crossings (as a function of
$\unicode[STIX]{x1D707}_{p}$
) are seen in the effective friction coefficient, again occurring in the frictional sliding region. Here the frictionless and large aspect ratio ellipses have a lower
$\unicode[STIX]{x1D707}_{c}$
compared to the disc case, while the reverse is true for very frictional particles.
One sees, furthermore, that very high microscopic friction coefficients are followed with a breakage of the directional ordering (seen both in
$\langle \unicode[STIX]{x1D703}_{\boldsymbol{e}\boldsymbol{\cdot }\hat{\boldsymbol{y}}}\rangle _{c}$
and
$S_{2}^{c}$
). This is further evident from figure 8 which shows three snapshots of for
$\unicode[STIX]{x1D6FC}=3$
at three different microscopic friction coefficients. This seems to have consequences for the normal stress differences which drastically decrease when the orientational order decreases and then to recover again at even higher friction coefficients.
Three typical snapshots for
$\unicode[STIX]{x1D6FC}=3$
and (a)
$\unicode[STIX]{x1D707}_{p}=0$
, (b)
$\unicode[STIX]{x1D707}_{p}=0.4$
and (c)
$\unicode[STIX]{x1D707}_{p}=10$
at
$I\sim 10^{-4}$
. Wall particles are coloured red and flowing particles blue. A darker blue indicates a larger number of contacts (colour bars to the right).

The theoretical value of 3 for the number of contacts per particle is reached for all cases, irrespective of the aspect ratio, as soon as the friction coefficient reaches a value above approximately 4. This is in contrast to the frictionless cases which tend to be in hypostatic conditions. The contact angle for the ellipses moves progressively towards the major axis as the aspect ratio increases. This is true both for frictional and non-frictional grains, but which is more enhanced for the former and least enhanced around friction coefficients of
$\unicode[STIX]{x1D707}_{p}\simeq 0.4$
.
3.7 Realistic friction coefficients
Measured microscopic friction coefficients between different materials usually range between 0.1 and 1 and can be extracted from macroscopic measurements (of discs or spheres). To better see the ‘realistic’ regime(s) all figures 6–7 have those zones shadowed in blue, showing the possible experimental spread.
Interpreting the microscopic friction can be done either from estimates of both the packing fraction or the effective friction (from e.g. avalanche angles) or both, ignoring the slight dependencies on polydispersity and the boundaries. Our data show that at least avalanching angles might be difficult to use for elongated ellipses, as one can not uniquely determine the microscopic friction from the angle in this realistic regime, see figure 7(e), and for the intermediate aspect ratios. This can, however, be cured if one also looks at the directional ordering after such an avalanching event.
4 Discussion and conclusion
We have shown that the
$\unicode[STIX]{x1D707}(I)$
-rheology can be extended to both frictionless and frictional ellipses yielding power laws in the various quantities similar to the isotropic disc and sphere cases. We find that the shear jamming and rheology of assemblies of ellipses heavily depend on the microscopic friction and that the difference between the critical packing fractions of frictionless and frictional ellipses increases with the elongation. While for frictionless ellipses one finds denser and denser shear-jammed packings as the elongation grows the same is not true for ellipses of intermediate or high microscopic friction. At realistic microscopic friction coefficients (
$\unicode[STIX]{x1D707}_{p}\simeq 0.4$
) the behaviour clearly becomes non-monotonic. Similar observations, following a pouring protocol and closer related to the RCP, have previously been found for monodisperse ellipses (Guises et al.
Reference Guises, Xiang, Latham and Munjiza2009). The packing fraction at shear jamming seems to be monotonic for low friction coefficients. Conversely, the same is not true for the effective friction which shows a non-monotonic behaviour with a maximum around
$\unicode[STIX]{x1D6FC}^{max}\simeq 1.4$
, which is also where the packing fractions saturate to plateau values (as a function of
$\unicode[STIX]{x1D6FC}$
). We find that the directional ordering increases with the aspect ratio for most cases and slowly converges (as a function of
$\unicode[STIX]{x1D6FC}$
) with the long axis along the direction of the
$x$
-axis of the cell. The exceptions are the very frictional cases which only reach a moderate ordering with a slight preference to align along the
$x$
-axis. For ellipses at fixed aspect ratios one finds the highest alignment for moderate friction coefficients of around
$\unicode[STIX]{x1D707}_{p}\simeq 0.1$
, an alignment that is followed by an increased normal stress difference.
The above results and behaviours can be rationalised by having three distinct regions of dissipation (by different motions): collisional, frictional sliding and frictional rolling. These regions have previously been identified both in granular and suspension flows of discs (DeGiuli et al.
Reference Degiuli, McElwaine and Wyart2016; Trulsson et al.
Reference Trulsson, DeGiuli and Wyart2017). As particles go from a collisional to a frictional sliding regime the dissipation increases and the effective friction increases. Simultaneously the number of contacts and hence also the packing fraction decreases as less contacts are needed as more and more of the contacts get stabilised by non-sliding tangential forces. This explains essentially the results for the low to moderate friction coefficients in figure 7 for all aspect ratios. At high friction coefficients, however, one sees a decrease in both effective critical friction and directional ordering for all ellipses (expect possibly for the tiniest elongations). This can be understood as a cross-over from a frictional sliding to a frictional rolling regime. The rolling regime has previously been found to occur at friction coefficients close to
$\unicode[STIX]{x1D707}_{p}\sim 1$
(DeGiuli et al.
Reference Degiuli, McElwaine and Wyart2016) and is consistent with our result
$\unicode[STIX]{x1D707}_{p}\sim 1$
for ellipses. While for discs (and ‘nearly’ discs) one does not see any signature of this in the macroscopic quantities, such as
$\unicode[STIX]{x1D707}_{c}$
, one clearly sees the signature in the ellipse cases due the fact that the ellipses have different axes and a rotation hence yields an excluded area different compared to that of the disc. A revolution of an ellipse would, hence, generate an effective excess excluded area proportional to
$(\unicode[STIX]{x1D6FC}-1)$
bigger than for the disc, but still have the same real area. One could, hence, expect a scaling of the packing fraction as the disc packing minus
$(\unicode[STIX]{x1D6FC}-1)$
times a constant (of order of unity). As neither of these are found it suggests that motions and orientations are correlated even for highly frictional particles. In the rolling regime we find that the number of contacts (per particle) approaches the theoretical value of 3, i.e. close to isostaticity, for static packings. Even though the sheared system is never strictly static (even in the quasistatic limit), the system seems to go from hypostatic to isostatic packings as friction is increased also for shear-jammed systems.
Time-averaged profiles (dotted lines) of (a) packing fraction
$\unicode[STIX]{x1D719}(y)$
, (b) velocity
$u(y)$
, (c) pressure
$P(y)$
and (d) stress
$\unicode[STIX]{x1D70E}_{xy}(y)$
(inset:
$\unicode[STIX]{x1D707}(y)$
) along the
$y$
-axis for
$\unicode[STIX]{x1D6FC}=3$
and
$\unicode[STIX]{x1D707}_{p}=0.4$
. The sampled region (excluding
$5d$
layers close to each wall) is shadowed in blue. Solid lines correspond to the spatial averages within the sampled region, except for the velocity profiles where they give the best slope.

In conclusion, we find that frictionless and very frictional ellipses behave differently at the yielding/jamming point. While one shear jam frictionless ellipses at denser packing fractions than discs the opposite is true for ellipses with high friction coefficients. We find that the explanation for this is the appearance of rolling motion of the ellipses rather than sliding while under shear. We furthermore find that the isostaticity seems to be fulfilled in the ‘infinite’ frictional regime which is in stark contrast to the case of hypostaticity in the frictionless cases. While more elongated ellipses lead to higher directional ordering this seems to be broken once the rolling motion sets in. For realistic friction coefficients one partly has frictional rolling and partly sliding motion. Increasing or decreasing the friction coefficient therefore changes the balance of these and leads to non-trivial non-monotonic behaviours in both effective frictions as well as directional ordering.
These results can both help us understand the onset of avalanching in granular material and the self-organised criticality (Frette et al. Reference Frette, Christensen, Malthe-Sorenssen, Feder, Jossang and Meakin1996) in those flows. It further gives insight to the directional ordering in flows (Börzsönyi et al. Reference Börzsönyi, Szabó, Törös, Wegner, Török, Somfai, Bien and Stannarius2012) and possible segregation effects in such flows (Windows-Yule et al. Reference Windows-Yule, Scheper, den Otter, Parker and Thornton2016).
It is natural to raise the question regarding how these critical packing properties at shear jamming are related to other protocols, e.g. isotropic compaction (i.e. RCP), as different protocols could either enhance or reduce these rotational or sliding motions and/or invoke various anisotropies in the system. This, as well as detailed studies of the force and contact networks, will be the subjects of future work.
Acknowledgements
We thank P. Claudin, J. J. Crassous and J. Stenhammar for useful discussions. This work was financial supported by the Swedish Research Council (621-2014-4387). The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the centre for scientific and technical computing at Lund University (LUNARC).
Spatially averaged profiles (in the sampled region) of (a) packing fraction
$\unicode[STIX]{x1D719}(t)$
, (b) shear rate
$\dot{\unicode[STIX]{x1D6FE}}(t)$
, (c) pressure
$P(t)$
and (d) stress
$\unicode[STIX]{x1D70E}_{xy}(t)$
(inset:
$\unicode[STIX]{x1D707}(t)$
) as a function of time
$t$
for
$\unicode[STIX]{x1D6FC}=3$
and
$\unicode[STIX]{x1D707}_{p}=0.4$
.

Appendix A. Supplementary information
A.1 Profiles
Spatial time-averaged profiles of packing fraction
$\unicode[STIX]{x1D719}$
, velocity along the
$x$
-axis (denoted
$u$
), pressure
$P$
and stress
$\unicode[STIX]{x1D70E}_{xy}$
can be seen in figure 9 for a few inertial numbers for
$\unicode[STIX]{x1D6FC}=3$
and
$\unicode[STIX]{x1D707}_{p}=0.4$
. The profiles are time averaged over ten strains, with strain
$\unicode[STIX]{x1D6FE}$
defined as
$\unicode[STIX]{x1D6FE}=\langle \dot{\unicode[STIX]{x1D6FE}}\rangle t$
. The
$y$
-axis is normalised by the average height
$h$
for each system, pressures and stresses by the corresponding imposed pressure at the walls
$P_{w}$
and the velocity
$u$
by the wall velocities
$\pm u_{w}$
. From figure 9 one sees that the packing fractions, pressures and stresses are all homogenous in the sampled region. The profiles have spatial fluctuations which are of the order of a few per cent of the means. Note that
$\unicode[STIX]{x1D707}$
-curves (inset) are much less noisy than both
$P$
and
$\unicode[STIX]{x1D70E}_{xy}$
(i.e.
$P$
and
$\unicode[STIX]{x1D70E}_{xy}$
are highly correlated). Velocity profiles are well described by linear profiles in the centre of the cell, i.e. the shear rate
$\dot{\unicode[STIX]{x1D6FE}}$
is to a good approximation constant over the whole sampled region irrespectively of
$I$
-value.
A.2 Time series
Figure 10 shows the time series of the instantaneously and spatially averaged
$\unicode[STIX]{x1D719}$
,
$\dot{\unicode[STIX]{x1D6FE}}$
,
$P$
, and
$\unicode[STIX]{x1D70E}_{xy}$
over 5 strains. Pressures and stresses are found to fluctuate more in the high
$I$
-region (where the particles have high kinetic energy) compared to low
$I$
-values, but again looking at
$\unicode[STIX]{x1D707}$
substantially reduces these fluctuations (i.e. fluctuations in pressure and stress are highly correlated also in time). The fluctuations in the estimated shear rates increase at lower
$I$
values compared to the higher ones. Similar behaviours are found for the other systems (other
$\unicode[STIX]{x1D6FC}$
and
$\unicode[STIX]{x1D707}_{p}$
-values). The data show no drift, validating that the simulations have indeed reached steady state. From the time series it is not evident if the system alternates between two (or more) shear states (i.e. shear localisation in time) or not.
A.3 Finite size effects
Finite size effects in
$\unicode[STIX]{x1D719}(I)$
and
$\unicode[STIX]{x1D707}(I)$
are shown in figure 11 for the case
$\unicode[STIX]{x1D6FC}=3$
and
$\unicode[STIX]{x1D707}_{p}=0.4$
. The two system sizes,
$N\sim 900$
and
$N\sim 1800$
, yield the same rheological curves,
$\unicode[STIX]{x1D719}(I)$
and
$\unicode[STIX]{x1D707}(I)$
, within numerical noise. The doubled system has the same dimension in the
$x$
-dimension as the smaller (
$N\sim 900$
), yielding an approximately doubled width (in
$y$
-direction) between the two rough walls compared to the latter.
Finite size effects in (a)
$\unicode[STIX]{x1D719}(I)$
and (a)
$\unicode[STIX]{x1D707}(I)$
for
$\unicode[STIX]{x1D6FC}=3$
and
$\unicode[STIX]{x1D707}_{p}=0.4$
.













