# Incremental shear strain chain: a mesoscale concept for slip lines in 2D granular materials

Jiaying Liu, Antoine Wautier, Wei Zhou, François Nicot, Félix Darve

► To cite this version:

Jiaying Liu, Antoine Wautier, Wei Zhou, François Nicot, Félix Darve. Incremental shear strain chain: a mesoscale concept for slip lines in 2D granular materials. Granular Matter, 2022, 24 (4), pp.1-39. 10.1007/s10035-022-01258-y. hal-03770142

HAL Id: hal-03770142

https://hal.inrae.fr/hal-03770142

Submitted on 6 Sep 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Incremental shear strain chain: a mesoscale
concept for slip lines in 2D granular materials

Jiaying Liu<fnref n="1" />,<fnref n="2" />, Antoine Wautier<fnref n="3" />*, Wei Zhou<fnref n="4" />, François Nicot<fnref n="5" />,<fnref n="6" /> and Félix Darve<fnref n="7" />

Abstract

Constitutive behaviors of granular materials are driven by both particle
interactions and geometric arrangements of contact network. To bridge
the gap between the grain scale and the sample scale, the mesoscale is of
great importance as it corresponds to the smallest scale at which geomet-
rical effect can be accounted for. Meso shear structures (sometimes called
microbands) have been observed frequently on incremental strain maps

<footnote n="1">¹Department of Civil Engineering, Zhejiang University City College, Huzhou Street, Hangzhou, 310015, China.</footnote>

<footnote n="2">²Zhejiang Engineering Research Center of Intelligent Urban Infrastructure, Huzhou Street, Hangzhou, 310015, China.</footnote>

<footnote n="3">³*INRAE, UMR RECOVER, Aix-Marseille University, Rte Cézanne, Aix-en-Provence, 13182, France.</footnote>

<footnote n="4">⁴State Key Laboratory of Water Resources and Hydropower
Engineering Science, Wuhan University, Bayi Road, Wuhan,
430072, China.</footnote>

<footnote n="5">⁵INRAE, UR ETGR, Université Grenoble Alpes, Rue de la Papeterie, St-Martin-d'Hères, 38402, France.</footnote>

<footnote n="6">⁶Laboratoire EDYTEM – USMB/CNRS, UMR 5204, Université Savoie Mont Blanc, Boulevard de la Mer Caspienne, Le Bourget-du-Lac, 73376, France.</footnote>

<footnote n="7">⁷CNRS, G-INP, Laboratoire 3SR, Université Grenoble Alpes,
UMR5521, Grenoble, 38000, France.</footnote>

<footnote>*Corresponding author(s). E-mail(s): antoine.wautier@inrae.fr;</footnote>

<footnote>Contributing authors: liujy@zucc.edu.cn; zw\_mxx@whu.edu.cn;</footnote>

<footnote>francois.nicot@inrae.fr; felix.darve@3sr-grenoble.fr;</footnote>
2 *Shear chain concept for granular materials*

in granular materials under shearing, while the Rudnicki and Rice local-
ization criterion for shear band is not fulfilled yet. These meso structures
are thin, quasi-linear and they involve a few grains as well as their sur-
rounding voids. This paper introduces the concept of “incremental shear
strain chain” (simply called “shear chain”) to provide a specific quantita-
tive definition of such mesostructures. “Shear chains” are defined based
on incremental deviatoric strain fields in 2D biaxial simulations. Particu-
lar attention is paid to demonstrate that the shear chain orientation is a
material scale property, insensitive to boundary conditions, loading paths
and sample densities. Since shear chains are shown to be closely related
to sliding mechanisms, they can stand for a mesoscale definition of the
concept of slip lines as defined in the standard elasto-plasticity theory.

**Keywords:** granular materials, DEM, mesostructures, slip line, incremental shear strain chain, granular plasticity

# 1 Introduction

Granular materials are found in a great variety of engineering problems and industrial applications. Despite their apparent simplicity, the understanding of their multiscale mechanical behaviors is a fascinating and very active domain of research. Constitutive behaviors of granular materials are complex and involve combined effects of local contact laws [1, 2], particle breakage [3, 4], particle size distribution [5, 6], fluid impact [7, 8] and contact network topology [9–12]. The contact-based and void-based fabric [13–15], the strong/weak contact subnetworks [12, 16, 17] and particle kinematics [18, 19] have been used to capture the influence of the geometrical arrangement of grains directly at the representative elementary volume scale. However, the links between the evolving microstructure and complex macroscopic mechanical behaviors are not yet well stabilized, even if local mechanisms responsible for failure and instability are increasingly more understood [20–25].

One of the most puzzling features for granular materials is the appearance
of shear bands. Even though shear bands relate to plastic dissipation and their
formation relies on Rudnicki and Rice localization criterion [26] in the con-
tinuum view, the micromechanical origin of shear band has not been clearly
understood yet. In recent years, the evolutions of particle rotations, local voids
and other micromechanical indices for granular materials have been tracked by
digital image correlation (DIC), X-ray tomography and synchrotron microto-
mography [27–30]. Besides, numerical tools, such as Discrete Element Method
(DEM), are proved to be relevant to perform numerical experiments on granu-
lar materials and access to an exhaustive view of their microstructure [31–35].
Such detailed descriptions of the deformation patterns in granular materials
have revealed the existence of small and short life shear patterns, long before a
persistent and well-marked shear band develops under deviatoric loads [30, 36–
38]. Kuhn [36] named these structures as “micro-bands” in DEM simulations,
and compared their orientation properties to other micromechanical evolu-
tions. Tordesillas et al. [39] correlated the “micro-bands” and the non-affine
deformations. In experimental observations, the “micro-bands” have also been
tracked and compared to the shear bands [40, 41]. However, the nature of
such structures is still debated and their links to shear bands remain to be
established.

At the macroscopic level, the shear band relates to obvious discontinuous
area in the strain field where failure occurs (Rudnicki and Rice localization
criterion [26]), which exhibits a size quite larger than that of slip lines. Slip
lines are accounted for in the standard continuum mechanics’ framework [42],
in which the finite thickness of shear band is frequently considered [43–45].
For granular materials, meso shear structures similar to slip lines exist under
shearing with a size of a few particles, which differ but correlate to the macro
shear bands [38, 45–48]. To distinguish shear bands from local shear patterns
in frictional granular materials, an intermediate scale, i.e., meso scale, between
grain and sample scales is needed. At this scale, slip structures can be analyzed
from both continuum and discrete points of views, provided that local stresses
and strains are consistently defined [49–51].

Researchers tried to model the small shear patterns observed in granular
samples. For example, [52–54] proposed micro-slip models of strain localization
in sand deformation and Shi and Horii [52] even formulated the evolution of
individual micro slips. Koenders [55] derived theoretical formulations of local
shear behaviors, and Gaspar and Koenders [56] proposed an interpretation of
microband involving heterogeneous elasticity theory. Zaiser and Aifantis [45]
proposed to account for shear band formation as avalanches of plastic slips. Le
Bouil et al. [37] characterized the length scale and persistence of microbands
thanks to digital image correlation in experimental setups. They found that
the final shear band does not arise from a coalescence of microbands, nor
is initiated by a single micro band that reaches the boundary and becomes
locked. Karimi and Barrat [40] used an elasto-plastic model of solid flow (based
on Mohr-Coulomb yield surface), to interpret microscopic correlations. These
models were established in the framework of continuum mechanics and did
not provide direct links to the underlying granular mesostructures. In a recent
work, Zhou et al. [38] used a new metric based on the concept of local intrinsic
dimensionality to provide quantitative measures of kinematic patterns in gran-
ular assemblies. Their results suggest that a complex symbiosis exists between
microbands and shear bands in granular material.

In this work we propose to adopt a mesomechanical point of view could
to relate the microband concept to mesostructures of a few grains instead of
using coarse graining approach [57] to work with continuum fields. At such a
scale, a few grains and associative voids consist of mesostructures (e.g., force
chains [58] and contact-based loops [36, 59, 60]) have been defined and inves-
tigated. Based on the contact and particle information, incremental deviatoric
strain fields can be defined and qualitatively visualize deformation features of
the whole sample[34, 36]. In previous work [34, 50], the authors have shown
4 Shear chain concept for granular materials

that meso shear structures are clearly visible in incremental deviatoric strain
fields under the early stage of biaxial shearing. Following these observations by
naked eyes, we introduce the concept of *incremental shear strain chain* (sim-
ply called *shear chain* in the following for the sake of brevity) to capture these
meso shear features from a quantitative point of view in this paper.

This paper is organized as follows. In Section 2, an example of 2D DEM sim-
ulation is presented and mesoscopic deformation quantities attached to grain
loops are shown. In Section 3, the concept of shear chain is introduced and the
detection algorithm is specified. Some statistical features of the shear chains
are also explored. In Section 4, a number of complementary DEM simulations
(proportional strain loading tests, drained tests of different confining pressures
and densities, drained tests of different sample aspect ratios) are exhibited in
order to demonstrate that the shear chain orientation is a material scale prop-
erty. Finally, in Section 5, a detailed discussion on the relationship between
the discrete shear chain concept and the slip line in continuum mechanics for
granular materials is given. Particular attention is given to sliding dissipation
and shear chain orientation with respect to slip line directions as derived from
non-associated plasticity theory.

2 From discrete to continuum views of granular materials

In this section, a DEM biaxial test is presented and the construction of the
incremental strain maps during the loading process is recalled. This simple
biaxial DEM test shows how meso shear structures appear and develop under
shearing, until the formation of a persistent shear band.

## 2.1 Drained biaxial test numerical modeling

Quasi-2D<fnref n="1" /> biaxial tests are carried out numerically with the use of the open-source DEM software YADE [61]. Granular samples are generated within a rectangular box of aspect ratio 1.5 with rigid boundaries (as shown in Fig. 1(a), *L*/*W* = 1.5, *L* represents the length and *W* denotes the width) containing a single layer of 20,000 spheres with a uniform distribution of diameters. The average diameter *D*₀₅₀ = 0.008 m and *D*<sub>max</sub>/*D*<sub>min</sub> = 1.98. Then they are slowly compressed to reach an isotropic stress state under the confining pressure of 100 kPa<fnref n="2" />. During the compression process, the inter-granular friction angle is temporarily reduced to 2° to reach a relative dense state. After the isotropic loading, the drained loading condition is considered by keeping the lateral confining pressure constant to 100 kPa and imposing a constant compression strain rate *ε*̇₁ = 0.01 s⁻<fnref n="1" /> in the axial direction (the direction of Y axis). The parameters used in DEM simulations are listed in Table 1, where *ρ* is the grain

<footnote n="1">¹Particles with a 3D shape are considered (spheres in the present case) but their motion is restricted to 2D (X-Y plane, two translations and one rotation). The difference with fully 2D simulations lies in the definition of the inertia of the particles.</footnote>

<footnote n="2">²Please note that despite the fact that samples are in 2D, stresses are derived by use of an out of plane arbitrary dimension taken equal to the maximum particle diameters.</footnote>
density, *n* is the porosity of the granular assembly, *Z*<sub>m</sub> is the average initial coordination number, *k*<sub>n</sub> and *k*<sub>t</sub> are the normal and tangential stiffness of the contact model respectively, *R*<sub>s</sub> = 2*R*<sub>1</sub>*R*<sub>2</sub>/*(R*<sub>1</sub> + *R*<sub>2</sub>)* is the harmonic average of the radii of the particles in contact and *φ*<sub>*g*</sub> is the contact friction angle between spheres. The numerical damping is used in the simulations, which will adjust forces on particles by considering the current acceleration and particle velocity, details can be found in Šmilauer et al. [61].

<table><caption>**Table 1** Initial isotropic states and loading parameters of the biaxial test example of dense sample.</caption><thead><tr><th>ρ (kg·m<sup>-3</sup>)</th><th>n</th><th>Z<sub>m</sub></th><th>k<sub>n</sub>/R<sub>s</sub> (MPa)</th><th>k<sub>t</sub>/k<sub>n</sub></th><th>φ<sub>g</sub></th><th>Damping</th></tr></thead><tbody><tr><td>3,000</td><td>0.161</td><td>4.01</td><td>300</td><td>0.5</td><td>35°</td><td>0.25</td></tr></tbody></table>

Figure 1(b) presents the macroscopic stress/strain responses of the biaxial test (rigid boundaries), where *q* = 0.5(*σ*<sub>1</sub> − *σ*<sub>2</sub>) and *ε*<sub>v</sub> = *ε*<sub>1</sub> + *ε*<sub>2</sub><sup>3</sup>. The biaxial loading can be considered as quasi-static as the dimensionless ratio between the mean resultant force acting on grains and the mean contact force remains always below 0.01 [62]<sup>4</sup>. Key states are labeled under shearing: the initial state (A), the end of the linear-elastic period (B), the most contractive state (C), the peak of deviatoric stress (D), the onset of the single shear band forming (E), the softening regime (F)<sup>5</sup>. At State D (stress peak or limit stress state), the macro shear band firstly occurs and then develops until the ultimate state. In this paper, we mostly focus on the period before State D.

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Biaxial loading condition (a) and deviatoric stress and volumetric strain evolutions of an example of the biaxial test (b). In (b), different states are marked: A—the initial state, B—the end of the linear-elastic period, C—the most contractive state, D—the peak of deviatoric stress, E—the onset of the single shear band forming, F—the softening regime. Thick parts of the curves are key periods concerned in this paper.</figcaption></figure>

<footnote n="3"><sup>3</sup>Note that soil mechanics conventions are used with compressions counted positive.</footnote>

<footnote n="4"><sup>4</sup>This condition was also checked in the present study by carrying a sensitivity analysis of the stress strain response to the loading parameters</footnote>

<footnote n="5"><sup>5</sup>Details can be found in Liu et al. [34].</footnote>
6 Shear chain concept for granular materials

## 2.2 Incremental deviatoric strain fields attached to grain loops

Incremental deviatoric strain maps are constructed from the unique loop tessellation of the granular assembly. The grain loops in 2D granular materials are defined by dividing the contact network into enclosed polygons containing deformable or undeformable voids [33, 36, 39, 59]. An example of the 2D loop tessellation is shown in Figure 2(a). The use of grain loops to construct incremental strain fields (instead of coarse graining approach [57, 63] for instance) is motivated by existing results highlighting that grain loops are relevant mesoscale structures to capture the mechanical strength of granular materials [33, 36, 49, 50, 64]. It enables keeping alive the link with the underlying contact network at the microscopic scale, and correlating meso and micro-scales (see Section 5.1).

<figure><img src="image_3.png" /><figcaption>Fig. 2 Loop tessellation within the granular assembly (a) and schematic diagram for the computation of the strain tensor in loops (b).</figcaption></figure>

Based on the loop tessellation, different local incremental strain tensors can be attached to each loop cell in order to draw an incremental strain map [36, 65–67]. In this paper, loop incremental strain is defined as the average of the symmetric part of the displacement gradient (neglecting grain rotations) over the loop domain:

$$ \epsilon_{ij}^{loop} = \frac{1}{|L|} \int_{L} \frac{u_{i,j} + u_{j,i}}{2} \mathrm{dS} \quad (1) $$

where $|L|$ is the area of the loop domain $L$ and $u_{i,j} = \frac{\partial u_i}{\partial x_j}$ is the gradient of displacement field within $L$.

Thanks to Gauss theorem, the integration can be changed to the loop boundary $\partial L$:

$$ \epsilon_{ij}^{loop} = \frac{1}{|L|} \int_{\partial L} \frac{u_i n_j + u_j n_i}{2} \mathrm{d}l \quad (2) $$
in which **n** is the outer normal to ∂L. At the microscale, the notion of continuous displacement field is meaningless and only grain displacements are known. By assuming a linear interpolation of the displacement along the loop edges, the incremental strain tensor of loop L is eventually defined as:

$$ \epsilon_{ij}^{loop} = \frac{1}{|L|} \sum_{k=1}^{c} \frac{1}{2} l^k \left( n_j^k \frac{u_i^{1k} + u_i^{0k}}{2} + n_i^k \frac{u_j^{1k} + u_j^{0k}}{2} \right) \quad (3) $$

where **u**<sup>1k</sup> and **u**<sup>0k</sup> are the displacement vectors of the edge k vertices. The notations used in Equation (3) are summarized in Figure 2(b). As discussed in Deng et al. [68] for instance, this strain definition is consistent with the macroscopic strain definition through spatial averaging.

During the biaxial loading, incremental strain tensors dε<sup>loop</sup> for loop cells within the granular assembly are computed from the incremental grain displacements corresponding to macroscopic strain increments of dε<sub>1</sub> = 2.0*10<sup>-4</sup><fnref n="6" />. The strain of the granular sample has been tested, which is quite close to the average value of meso strains within the domain. From the computation of the incremental strain tensors according to Equation (3), incremental deviatoric strains are defined as the absolute difference between the two principal incremental strain values (dε<sup>loop</sup><sub>d</sub> = |dε<sup>loop</sup><sub>1</sub> - dε<sup>loop</sup><sub>2</sub>| where dε<sup>loop</sup><sub>1</sub> and dε<sup>loop</sup><sub>2</sub> are the eigenvalues of dε<sup>loop</sup>). For State B and C, small and non-persisting meso shear structures can be observed in the incremental deviatoric strain fields of the granular sample, as shown in Fig. 3 (the incremental deviatoric strain fields of other states can refer to Liu et al. [34]). The observed meso shear structure is similar to existing research of microband [36], in this paper we will have a specific definition of this structure based on the loop tessellation (Section 3).

In Fig. 3, meso shear patterns are visible and distribute uniformly in the sample domain with two symmetric preferred orientations. By zooming in the selected area of State B in Fig. 3(a), Fig. 4 highlights that the meso shear patterns consist in thin line like structures composed of loops with relatively large incremental deviatoric strain (mainly green loops). In Fig. 4(a), the direction of the incremental displacement of the particles is shown with unit arrows, demonstrating very little change when crossing the meso shear structures colored green. And in Figure 4(b), dark particles are force chain members (details of the definition in Appendix A) which are mostly excluded by meso shear structures.

# 3 The shear chain concept

From the observations in Fig. 4, the width of the meso slip structures corresponds to the size of one meso-loop. Therefore, it is tempting to assume that the meso shear patterns (microband) can be regarded as connected chains of

<footnote n="6">⁶Note that the results obtained in the present study are not impacted by the particular choice for the macroscopic strain increment. Indeed, all strain increments presented in the study can be rescaled in terms of strain rate. The use of grain displacements instead of velocities to compute strain rates enable to get rid of small oscillations resulting from the elastic stiffness at contacts.</footnote>
8 *Shear chain concept for granular materials*

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Examples of shear structure distribution within the 2D granular assembly [34]: (a) State B, $\epsilon_1$ = 0.004; (b) State C, $\epsilon_1$ = 0.009. The yellow rectangle domain in (a) is enlarged in Figure 4</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Zoomed images of a mesoscopic shear structure from State B (selected area is marked in Fig. 3(a)). The direction of the displacement field is represented in subfigure (a) while force chain particles are displayed in subfigure (b). The mesoscopic shear structures are mostly colored in green which correspond to large deviatoric incremental strain levels.</figcaption></figure>

loops of large incremental deviatoric strains. In the following sections, these structures are named *incremental shear strain chain* or simply *shear chains* for brevity, and are expounded with a rational definition and algorithm. Inspired
by force chain definition [58] recalled in Appendix A, this section will intro-
duce the concept of *incremental shear strain chain* from a quantitative way, to
describe the meso shear structures observed in incremental deviatoric strain
field before shear banding.

**3.1 Definition of incremental shear strain chain**

To capture the incremental shear strain chains, the maximum shearing direc-
tion of the incremental deviatoric strain of each loop is considered. Fig. 5 shows
the schematic drawing of a shear chain, and the procedure is listed in Fig. 6.
In detail, the steps to define shear chains are as follows:

1. For each mesoloop, compute the incremental strain tensor $d\mathbf{e}^{\text{loop}}$ using the strain definition of loops in Section 2, and identify meso loops with incremental deviatoric strains $d\varepsilon_d^{\text{loop}} = |d\varepsilon_1^{\text{loop}} - d\varepsilon_2^{\text{loop}}|^7$ larger than the average incremental deviatoric strain $\langle d\varepsilon_d^{\text{loop}} \rangle$.

2. Within this group of loops, randomly select one head loop (only one connection with other loops of large incremental deviatoric strain) and find the possible chain members. This is done by computing the maximum shearing directions<fnref n="8" /> of the chosen loop (called Loop A) and its connection loop (called Loop B) and then checking whether the geometric directions joining the 2 barycenters are close to the shearing directions of the two loops. In Fig. 5, $\beta$ denotes this deviation angle. The connection between the two loops are built if $\beta$ fulfills a given limitation, in this paper it is set as $\beta \le \beta_{th} = 22.5^\circ$. After one connection built, a new cycle of identifying chain members begins with the Loop B of the last step, until the threshold $\beta \le \beta_{th} = 22.5^\circ$ is not fulfilled or no connection loops are found.

3. If a chain of loops is composed of at least three loops, then it is regarded as a shear chain.

It is important to note that the shear chain does not correspond to the
shear band (as explained in Liu et al. [34]) at the macroscopic level because
it cannot be thicker than a mesoloop, and it does not correspond to sliding
contacts at the microscopic scale as it is defined on averaged incremental grain
displacements only (no interparticle slip is considered).

Similarly to force chain definition, the shear chain definition relies on three
thresholds, the combination of which prevents the identification of shear chains
in arbitrary deformation field. Because the shearing directions consist in two
perpendicular directions, the linearity threshold is set to 45°/2 = 22.5°. The
choice to use the mean incremental deviatoric strain is justified by sensitivity
analyses presented in B.1. Moreover, the robustness of this definition is care-
fully shown in B.2 (in particular the random selection of the head loops is
shown not to significantly affect the identified shear chain population). Besides,

<footnote>$^7$d$\varepsilon_1^{\text{loop}}$ and $d\varepsilon_2^{\text{loop}}$ are the eigenvalues of $d\mathbf{e}^{\text{loop}}$</footnote>

<footnote n="8">⁸The shearing directions correspond to the frame in which the deviatoric part of $d\boldsymbol{\epsilon}^{\text{loop}}$ has zero diagonal terms. They are oriented to $\pm 45^\circ$ from the principal strain directions. In Fig. 5 the two shearing directions are shown with red dashed lines</footnote>
10 Shear chain concept for granular materials

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Schematic drawing of an incremental shear strain chain. Shearing directions are displayed in relation with the principal incremental strain directions and the geometric direction joining the loop barycenters.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6** The flow chart of searching incremental shear strain chains.</figcaption></figure>
the shear chain definition proposed in this paper relies on the deviatoric part of the incremental strain to focus on shearing mechanisms <fnref n="9" />.

An example of the application of the shear chain detection algorithm is given in Fig. 7 between State A and State B for the drained biaxial test presented in Section 2. It can be seen that meso slip features are captured by the shear chain definition, and a number of shear chains with different colors are evenly distributed within the sample.

<figure><img src="image_8.png" /><figcaption>Fig. 7 Shear chains within the granular assembly (between State A and State B). Quasi-linear structures of different colors denote various shear chains.</figcaption></figure>

By definition, shear chains are line-like structures. A representation of shear chains is achieved by joining the geometric barycenters of adjacent cells (dark thick line in Fig. 5). For a shear chain of n adjacent mesoloops, the least square method defines the mean direction θ as:

$$ \theta = \arctan \left( \frac{\sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^{n} (x_i - \bar{x})^2} \right) \quad (4) $$

where $x_i$ and $y_i$ corresponds to the x and y coordinates of the geometric barycenter of loop cell i, and $\bar{x}$ and $\bar{y}$ are the averages of the $x_i$ and $y_i$, respectively.

<footnote n="9">⁹Should the total incremental strain be considered instead of its deviatoric part (to account for volume variations), an alternative definition for shear chains might be sought by considering the perpendicular direction to the zero-extension direction (if such a direction exists)</footnote>
## 3.2 Statistical features of shear chains

Following the algorithm suggested in Section 3.1, all the states of the biaxial test presented in Section 2 are considered to identify shear chains within the granular assembly, especially for the states before the shear banding (state D). The spatial distribution of shear chains is shown in Fig. 7 quite near the initial isotropic stress state. Hundreds of shear chains are randomly distributed within the granular assembly, demonstrating that the mesoscopic shearing occurrence is diffuse and disordered in space at the first stage of deviatoric loading.

The probability distribution of the shear chain orientation $\theta$ corresponding to the state of Fig. 7 is shown in Fig. 8(a). Shear chains are more likely to orientate in one of the two symmetric directions $\pm 45^\circ$, which is similar to the microband observation at the first stage of deviatoric loading in literature [36, 41].

<figure><img src="image_9.png" /><figcaption>**Fig. 8** Probability distributions of shear chain orientation $\theta$ for $\epsilon_1$ = 0.001 (close to State A) (a) and for $\epsilon_1$ = 0.005 (close to state B) (b) of the granular sample in Section 2.1.</figcaption></figure>

As the biaxial loading is pursued, the probability distribution of the shear chain orientation $\theta$ keeps the same double-peak feature but the preferred orientation seems to change and larger fluctuations are observed, which can be found in Fig. 8(b) showing the probability distribution of $\theta$ at the state $\epsilon_1$ = 0.005. These observations are consistent with the results of Houdoux et al. [41]. In their experimental work, an increase in the size of the error bars is observed for the orientations of the fast fluctuation part of the spatial autocorrelation before shear banding occur (see Fig. 7 in [41]).

Since the preferred orientations of shear chains are symmetric, then the absolute value $|\theta|$ could be a proper indicator to track the orientation evolution. The evolution of its mean value $<|\theta|>$ is plotted in Fig. 9 for the period before the stress peak (state D, the onset of macroscopic shear band). At the beginning of the loading, $<|\theta|>$ is around $45^\circ$. As the deviatoric loads continue to be applied, $<|\theta|>$ increases gently (shear chains get closer to the principal compression direction). As illustrated qualitatively in Fig. 8, the standard deviation around $<|\theta|>$ is increasing along the drained biaxial
test. Nevertheless, a statistical analysis detailed in C demonstrates that the increasing trend visible in Fig. 9 is statistically representative.

<figure><img src="image_10.png" /><figcaption>Fig. 9 Evolution of the average shear chain direction $<|\theta|>$ during biaxial drained loading.</figcaption></figure>

As a result, the shear chain orientation is influenced by the stress state and/or possibly by the type of loading conditions and the change in geometry of the deformed sample. Before discussing the change in orientation of shear chains with the framework of non-associated elasto-plasticity theory in Section 5, the invariance of the shear chain orientation with respect to structural effects is investigated in Section 4.

# 4 Shear chain orientation under different loading conditions

In Section 3, we have defined the shear chain and analyzed the orientation evolution considering a standard drained biaxial test. In this section, different loading conditions are considered to highlight a unique relationship between shear chain orientation and the stress state of the granular assembly$^{10}$. Results in this section provide evidences that the shear chain orientation is a material scale property.

## 4.1 Drained biaxial tests with varying sample aspect ratios

Similar to the example we analyze in Section 2, a set of quasi-2D biaxial drained tests of samples with different aspect ratios ($L/W$ in Fig. 1(a)) are conducted.

<footnote>$^{10}$There might be an additional dependence to some constitutive microstructure features too (e.g. the fabric tensor). Such a study is however out of the scope of the present study.</footnote>
14 Shear chain concept for granular materials

Firstly, all specimens are randomly generated within rectangular domains of prescribed aspect ratios and compacted to isotropic confining state of 100 kPa.

The same uniform distribution of DEM sample in Section 2, as well as other DEM parameters given in Table 1 are used for biaxial samples of different L/W. Additional specific parameters are summarized in Table 2. For the sixteen samples, only the aspect ratio L/W and the particle number are different at the initial compacted states. The slight changes in void ratio $V_r$ can be neglected. When conducting the drained biaxial loading, a compression is imposed in the vertical direction with a constant strain rate of $\varepsilon_1 = 0.01 \text{ s}^{-1}$ and the lateral pressure is maintained constant to $\sigma_2 = 100 \text{ kPa}$. The overall stress/strain responses for different boundary conditions (different L/W) are shown in Fig. 10. The sample in Section 2 corresponds to Sample DF here.

<table><caption>**Table 2** Sample parameters for the series of drained DEM tests. The sample presented in Section 2 is Sample DF.</caption><thead><tr><th>Sample</th><th>L/W</th><th>N<sub>p</sub></th><th>Initial void ratio V<sub>r</sub></th></tr></thead><tbody><tr><td>DA</td><td>1.0</td><td>13333</td><td>0.0162</td></tr><tr><td>DB</td><td>1.1</td><td>14667</td><td>0.0161</td></tr><tr><td>DC</td><td>1.2</td><td>16000</td><td>0.0161</td></tr><tr><td>DD</td><td>1.3</td><td>17333</td><td>0.0161</td></tr><tr><td>DE</td><td>1.4</td><td>18667</td><td>0.0161</td></tr><tr><td>DF</td><td>1.5</td><td>20000</td><td>0.0161</td></tr><tr><td>DG</td><td>1.6</td><td>21333</td><td>0.0161</td></tr><tr><td>DH</td><td>1.7</td><td>22667</td><td>0.0160</td></tr><tr><td>DI</td><td>1.8</td><td>24000</td><td>0.0161</td></tr><tr><td>DJ</td><td>1.9</td><td>25333</td><td>0.0160</td></tr><tr><td>DK</td><td>2.0</td><td>26667</td><td>0.0160</td></tr><tr><td>DL</td><td>2.1</td><td>28000</td><td>0.0160</td></tr><tr><td>DM</td><td>2.2</td><td>29333</td><td>0.0160</td></tr><tr><td>DN</td><td>2.3</td><td>30667</td><td>0.0160</td></tr><tr><td>DO</td><td>2.4</td><td>32000</td><td>0.0159</td></tr><tr><td>DP</td><td>2.5</td><td>33333</td><td>0.0159</td></tr></tbody></table>

Although the sizes of samples are different for the drained biaxial tests in Fig. 10, the stress-strain curves before $\varepsilon_1 = 0.01$ are identical. The stage $\varepsilon_1 = 0.01$ corresponds to a bifurcation point after which the response of the sample is strongly influenced by the boundary conditions. During the drained biaxial loading, the incremental deviatoric strain field of specimens experiences a transition from random-distributed meso shear chains to the macro shear band. Fig. 11(a) shows the incremental deviatoric strain fields at $\varepsilon_1 = 0.005$ for the samples DC (L/W=1.2) and DK (L/W=2.0) respectively, which exhibit almost the same spatial distribution of incremental deviatoric strain as for sample DF in Fig. 3. Similar incremental strain fields were obtained for the other samples (not shown here for the sake of brevity).

For each specimen with a given aspect ratio L/W, shear chains exist from the beginning of the loading until the shear band forms. The probability density distributions of shear chain orientations for sample DC (L/W=1.2) and
<figure><img src="image_11.png" /><figcaption>**Fig. 10** Evolutions of deviatoric stress $q$ (a) and volumetric strain $\varepsilon_v$ (b) during the drained biaxial tests with axial strain rate $\dot{\varepsilon}_1 = 0.01 \text{ s}^{-1}$ and lateral confinement $\sigma_2 = 100 \text{ kPa}$. Different curves correspond to different aspect ratios $L/W$ ranging from 1.0 to 2.5.</figcaption></figure>

DK (L/W=2.0) at the same state of Fig. 11(a) are given in Fig. 11(b). As sample DF, the preferred orientations of shear chains in Fig. 11(b) are also close to $\pm 45^\circ$.

The evolutions of the mean direction $|\theta|$ are shown for all the biaxial tests in Fig. 12. Despite some fluctuations, the orientation of shear chain is not affected by the boundary conditions: the preferred direction of shear chains increases from $45^\circ$ to nearly $50^\circ$ just before the stress peak and the onset of persistent shear bands.

## 4.2 Proportional strain loading path

In order to broaden the scope of the results of this paper, another type of loading path is considered in the form of proportional strain loading tests. We use a parameter $\alpha$ to control the change in volume of the sample during the loading as in Jrad et al. [69]. The incremental volumetric strain is set as:

$$d\varepsilon_v = d\varepsilon_1 + d\varepsilon_2 = \alpha d\varepsilon_1 \quad (5)$$

Under an axial compression ($d\varepsilon_1 > 0$), $\alpha < 0$ corresponds to an increase in the sample volume, while $\alpha > 0$ corresponds to an imposed contraction. $\alpha = 0$ corresponds to the special case of the undrained loading path (no volume change). Table 3 summarizes the values of $\alpha$ considered in the different tests in this paper. Except the loading path, all the initial conditions for samples in Table 3 are the same as for Sample DF in Table 2.

<table><caption>**Table 3** Sample parameters used for the proportional strain loading paths.</caption><thead><tr><th>Sample</th><th>PA</th><th>PB</th><th>PC</th><th>PD</th><th>PE</th><th>PF</th><th>PG</th></tr></thead><tbody><tr><td>α</td><td>-1.0</td><td>-0.2</td><td>-0.5</td><td>0.</td><td>0.2</td><td>0.5</td><td>1.0</td></tr></tbody></table>
16 Shear chain concept for granular materials

<figure><img src="image_12.png" /><figcaption>**Fig. 11** Examples of the incremental deviatoric strain field (a) and the probability distribution density (b) for Sample DC (L/W=1.2) and Sample DK (L/W=2.0), at the state of axial strain $\epsilon_1$ = 0.005.</figcaption></figure>

For the different proportional strain loading paths, the evolution of the generalized deviatoric stress $\sigma_1 - (1-\alpha)\sigma_2$<fnref n="11" /> is shown in Fig. 13. In all proportional loading paths, shear chains are clearly visible in incremental deviatoric strain fields and for some values of $\alpha$, and shear banding eventually occurs. Points in Fig. 13 denote the onset of shear banding. In Fig. 14, incremental deviatoric strain fields are shown for Sample PC ($\alpha$ = -0.2) and PG ($\alpha$ = 1) at $\epsilon_1$ = 0.002.

<footnote n="11">¹¹See details in Nicot et al. [70]. Under the constraint $d\epsilon_1 + d\epsilon_2 = \alpha d\epsilon_1$, the first order variation of the energy balance writes $dW = d\epsilon_1 [\sigma_1 - (1-\alpha)\sigma_2]$. Thus the generalized deviatoric stress is the conjugate variable to the axial strain. In case $\alpha$ = 0 it corresponds to the usual deviatoric stress.</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 12** Evolution of average preferred shear strain direction for specimens of different aspect ratios.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13** Evolutions of generalized deviatoric stresses $\sigma_1 - (1-\alpha)\sigma_2$ during the proportional strain loading with axial strain rate $\dot{\epsilon}_1 = 0.01 \text{ s}^{-1}$ and volumetric strain rate $\dot{\epsilon}_v = \alpha\dot{\epsilon}_1$. The points represent the beginning of strain localization.</figcaption></figure>

Since stress states diverge at a given axial strain for different proportional strain loading paths, the stress ratio $q/p = (\sigma_1 - \sigma_2)/(\sigma_1 + \sigma_2)$ is considered$^{12}$ to compare the mean orientation of shear chains in Fig. 15. In addition, the data corresponding to a drained biaxial test (Sample DF) is also presented in Fig. 15 for comparison. Although $q/p$ is the macroscopic description and $|\theta|$ is a mesoscale character of shear chains, they can be compared since the macro stress and strain of granular materials are statistical average of meso or micro

<footnote>$^{12}$The choice of such state variable is driven by the fact that granular materials are known to follow Mohr-Coulomb plasticity as a first approximation.</footnote>
18 Shear chain concept for granular materials

<figure><img src="image_15.png" /><figcaption>**Fig. 14** Examples of the incremental deviatoric strain field (a) and the probability distribution density (b) for Sample PC ($α = -0.2$) and Sample PK ($α = 1$), at the state of axial strain $ε_1 = 0.002$.</figcaption></figure>

stress and strains. The data corresponding to all the tests collapse on a single master curve which indicates that the shear chain orientation does not depend on the type of loading path but only on the stress state $q/p$.

## 4.3 Effect of density and confining stress

In this subsection, the impact of the initial density and the confining stress level are assessed by considering 9 drained biaxial tests for different initial porosities and confining pressures. Dense, medium and loose samples under the confining pressures of 100 kPa, 200 kPa and 300 kPa are named as “d100”, “d200”, “d300”, “m100”, “m200”, “m300”, “l100”, “l200” and “l300” respectively. Except the initial porosities and the confining pressures summarized in Table 4 parameters for the simulations are the same as Sample DF presented in Section 2.1.
<figure><img src="image_16.png" /><figcaption>**Fig. 15** Evolution of average preferred shear strain direction for samples under different proportional strain loading paths. Pink points corresponds to the drained biaxial test on sample DF (Section 2.1)</figcaption></figure>

<table><caption>**Table 4** Sample parameters of drained DEM tests of different densities and confining pressures.</caption><thead><tr><td>Sample</td><td>Initial porosity (2D)</td><td>Confining pressure</td></tr></thead><tbody><tr><td>d100</td><td>0.160</td><td>100 kPa</td></tr><tr><td>d200</td><td>0.154</td><td>200 kPa</td></tr><tr><td>d300</td><td>0.147</td><td>300 kPa</td></tr><tr><td>m100</td><td>0.185</td><td>100 kPa</td></tr><tr><td>m200</td><td>0.178</td><td>200 kPa</td></tr><tr><td>m300</td><td>0.172</td><td>300 kPa</td></tr><tr><td>l100</td><td>0.207</td><td>100 kPa</td></tr><tr><td>l200</td><td>0.198</td><td>200 kPa</td></tr><tr><td>l300</td><td>0.190</td><td>300 kPa</td></tr></tbody></table>

Samples of different initial densities and confining pressures exhibit different stress/strain relationships. Dense samples are characterized by strain hardening/softening and contractancy/dilatancy, while loose samples are related to the hardening and contractancy only. As conducted in the previous sections, shear chains are detected and analyzed for these nine drained tests. The evolutions of the mean $|\theta|$ versus the stress ratio $q/p$ are presented in Fig. 16. Despite several outliers, samples of different densities and confining pressures exhibit the almost the same relationship between the preferred direction of shear chains and the stress ratio. Combined to the results obtained in Sections 4.1 and 4.2, shear chain orientation is shown to be a material scale property of granular materials influenced by the stress state only (at least as a first approximation).

In addition, the stress-strain relationship and shear chain orientation evolution for samples of different coordination numbers but under the same
<figure><img src="image_17.png" /><figcaption>Fig. 16 Evolution of average preferred shear strain direction for drained biaxial tests with different initial porosities and confining pressures.</figcaption></figure>

confining pressure and porosity [71] should also be interesting to explore, to further capture the material property of the shear chain features.

# 5 Shear chains relationships to granular plasticity

The shear chain concept has been defined in Section 3, and the orientation of the shear chain has been shown to be a material scale property mainly influenced by the stress ratio $q/p$ in Section 4. In this section, we intend to correlate discrete plastic dissipation and continuum views of plasticity in granular materials in terms of the defined shear chain. The spatial correlation between contact sliding and shear chain is investigated and the shear chain orientation is analyzed in relation with the slip line orientations predicted by the standard elasto-plastic continuum formalism for non-associated materials.

## 5.1 Spatial correlation between shear chains and contact sliding

At microscale, energy dissipation can occur through contact sliding. According to the contact model used in the DEM simulations, a sliding index $I_p$ can be constructed for each contact as:

$$I_p = \frac{|f_t|}{|f_n| \tan \phi_g} \qquad (6)$$
where $|f_n|$ and $|f_t|$ are the normal and tangential force magnitudes respectively, and $\phi_g$ is the contact friction angle. As a consequence of local Coulomb criterion, contact sliding occurs when $I_p = 1$. Fig. 17 presents the sliding contact distribution for the 5 key states identified in Section 2.1.

<figure><img src="image_18.png" /><figcaption>Fig. 17 Spatial distributions of the contact sliding for different loading states (red points denote the positions of sliding, blue arrows show the tangential direction of the contact).</figcaption></figure>

Visually, the spatial distribution of sliding contacts can reflect the shear band characteristic (State D and E), however the relationship between contact sliding and meso shear features (in Fig. 3) is not obvious and has to be quantified from a statistical point of view. Note that the shear chain definition relies only on grain kinematics and thus the relation with contact sliding needs to be established.

Based on the shear chain definition, contact sliding analysis can be categorized to 3 groups:

- All contacts within the granular assembly;

- Shear chain contacts, i.e. all contacts of mesoloops composing the shear chains;

- Connection contacts in shear chains, denoting contacts connecting two adjacent loops in shear chains, which are marked in red in Fig. 5.

Evolutions of the sliding fraction $S_r$ and average sliding index $<I_p>$ for these three groups are shown in Fig. 18. The sliding fraction is defined as $S_r = N_s/N_c$, where $N_s$ denotes the number of sliding contacts among the set of contacts $N_c$ considered (all contacts, shear chain contacts or connection contacts).

Both $S_r$ and $<I_p>$ have the largest values for connection contacts. Shear chain contact population also exhibits larger $S_r$ and $<I_p>$ than the set of all contacts. Sliding occurs predominantly inside shear chains, which demonstrates that shear chains are mesostructures characterizing sliding mechanisms in granular materials. The fact that connection contacts are even more likely to slide than shear chain contacts agrees well with the existence of mesoscopic sliding lines oriented in the direction of the shear chain. Nevertheless, it should be acknowledged that only 15 % of the contacts are sliding, which means that
<figure><img src="image_19.png" /><figcaption>Fig. 18 Average sliding fraction $S_r$ (a) and average sliding index $<I_p>$ (b) for connection contacts, shear strain contacts and all contacts.</figcaption></figure>

the kinematic patterns defining shear chains do not result only from sliding in contacts but also probably from contact opening or simply loop deformation. Further micromechanical analysis will be carried out in forthcoming studies to clarify this point.

## 5.2 Upper and lower bounds for slip line orientation derived from elasto-plasticity theory

As a first approximation, granular plasticity is known to follow Mohr-Coulomb yield criterion. According to this yield criterion, plasticity activates at a given point if the maximum ratio between the shear and normal stresses is larger than a limit value. For a 2D stress tensor with principal stresses $\sigma_1 > \sigma_2 > 0$, the geometric direction giving rise to the maximum ratio between the shear stress and the normal stress is expressed as

$$ \theta_{stat} = \pm (45^{\circ} + \phi_{mob}/2), \qquad (7) $$

where $\phi_{mob}$ is the mobilized friction angle. It is related to the stress state through the relationship:

$$ \sin \phi_{mob} = \frac{\sigma_1 - \sigma_2}{\sigma_1 + \sigma_2} \qquad (8) $$

In the above equation, $\theta_{stat}$ is given with respect to the minor stress direction (horizontal direction in the present case). It defines the first direction in which sliding occurs when the mobilized friction $\phi_{mob}$ equals the internal friction angle $\phi$ (the Mohr-Coulomb yield criterion) <fnref n="13" />. In the present case, the major and minor stress directions coincide with the axes $y$ and $x$ respectively.

<footnote n="13">¹³Since the elastic domain is small for granular materials, the mobilized friction angle rapidly equals the internal friction angle. Then, its value vary along the loading path according to strain hardening/softening behavior.</footnote>
Thus, the angle $\theta_{stat}$ is given here with respect to the horizontal axis (like $\theta$ in Fig. 5).

If the yield function controls plasticity activation (the static condition), plastic strain develops according to a the plastic potential (the kinematic condition, aka the flow rule) that can be a different function in case of non-associated plasticity. For granular materials, the flow rule is usually characterized by the dilatancy angle $\psi$ (where $\psi > 0$ corresponds to dilatancy). This angle relates to the principal values of the incremental plastic strain as:

$$ \sin \psi = - \frac{d\varepsilon_1 + d\varepsilon_2}{d\varepsilon_1 - d\varepsilon_2} \qquad (9) $$

where $d\varepsilon_1$ and $d\varepsilon_2$ are the principal incremental strain values<fnref n="14" />. Provided that $d\varepsilon_1 > 0 > d\varepsilon_2$, there exists two directions of zero length variation. The perpendicular directions can be seen as slipping directions in the sense that no length variation is observed across the slipping direction. Expressed with respect to the minor stress direction, the slipping angle $\theta_{kinemat}$ satisfy the equation

$$ \mathbf{n} \cdot \begin{pmatrix} d\varepsilon_1 & 0 \\ 0 & d\varepsilon_2 \end{pmatrix} \cdot \mathbf{n} = 0 \qquad (10) $$

with $\mathbf{n} = \begin{pmatrix} \cos(\theta_{kinemat} - \pi/2) \\ \sin(\theta_{kinemat} - \pi/2) \end{pmatrix}$. This gives

$$ \cos(2\theta_{kinemat} - \pi) = \sin(-\frac{\pi}{2} \pm 2\theta_{kinemat}) = -\frac{d\varepsilon_1 + d\varepsilon_2}{d\varepsilon_1 - d\varepsilon_2} \qquad (11) $$

With use of Equation (9), the slipping directions are expressed with respect to the minor direction (horizontal direction here) as

$$ \theta_{kinemat} = \pm (45^{\circ} + \psi/2) \qquad (12) $$

For a pure incremental shear deformation, $d\varepsilon_1 + d\varepsilon_2 = 0$, and $\theta_{kinemat} = \pm 45^{\circ}$. In this case, zero extension directions coincide with slipping directions and the slipping directions are perpendicular. In the general case ($\psi \neq 0$), the slipping directions are not perpendicular.

When material failure occurs (i.e. when plastic strain develops), $\theta_{stat}$ and $\theta_{kinemat}$ “compete” to drive the failure line orientation. Indeed, the static condition (given by the yield surface) requires that sliding occurs along direction $\theta_{stat}$, while the kinematic condition (given by the plastic potential) predicts that slipping develops along direction $\theta_{kinemat}$. In case of associated plasticity $\theta_{stat} = \theta_{kinemat}$ and both conditions are equivalent, but this is not the case for non-associated plasticity.

At the meso scale, depending on whether the yield criterion or the flow rule controls plastic development, slip lines are also expected to form along

<footnote n="14">¹⁴Rigorously speaking, only the plastic part of the incremental strain should be considered, but in practice, the elastic incremental strain can be neglected as soon as the plasticity activates.</footnote>
24 Shear chain concept for granular materials

directions close to $\theta_{stat}$ or $\theta_{kinemat}$ [48]. In this paper, we have shown that the shear chain orientation is a material property. As a result, its orientation might be related to the yield criterion or flow rule at the macroscale (and not Rudnicki and Rice criterion). In Fig. 19, the mean orientation of shear chains found in this paper is compared to the slip line directions $\theta_{stat}$ and $\theta_{kinemat}$ estimated from the macroscopic stress and incremental strain tensors.

<figure><img src="image_20.png" /><figcaption>**Fig. 19** Comparison between the preferred orientation of shear chains with the slip line directions $\theta_{stat} = 45^{\circ} + \phi/2$ or $\theta_{kinemat} = 45^{\circ} + \psi/2$ predicted by plasticity theory.</figcaption></figure>

In Fig. 19, neither $\theta_{stat}$ nor $\theta_{kinemat}$ accounts alone for the shear chain direction. $\theta_{stat}$ overestimates the shear chain orientation while $\theta_{kinemat}$ underestimates the shear chain orientation.

Nevertheless, it can be underlined that the shear chain orientation is very close to $\theta_{stat}$ in the beginning of the test, while it becomes closer to $\theta_{kinemat}$ just before the onset of persisting shear band (State D). This observation is consistent with the fact that first slip lines appear while very few grain rearrangement has occurred, demonstrating that mesoscale failure is expected to be stress driven ($\theta_{stat}$). On the other hand, as the test is carried on, grains rearrange more and more, and the failure pattern is increasingly driven by the plastic strain development. Between the initial and failure states, the shear chain orientation remains in the range of [$\theta_{stat}$, $\theta_{kinemat}$]. It should also be underlined that the slipping directions $\theta_{stat}$ and $\theta_{kinemat}$ consider the macroscopic stress or strain whereas shear chain is a mesoscopic concept. As a result, it will be interesting to define and analyze meso stress and strain in further studies before being able to derive a constitutive relationship for the shear chain orientation.
## 5.3 Meso shear behaviors and Eshelby theory

In the previous subsection we show that, standard non-associated elasto-plastic formalism provides a lower bound and an upper bound for the shear chain orientation depending on whether the stress state or the plastic strain imposes the slip line direction. These two directions are related to the friction angle and the dilatancy angle separately (equations (7) (12)). Following the ideas proposed by Mc Namara et al. [72] and Karimi and Barrat [40], the orientation of slip lines can be derived from an Eshelby problem to relate the microband direction to both the friction and dilatancy angles. In their papers, these authors assume that a grain reorganization occurs in a given location in the granular assembly. This small reorganization is a consequence of a local collapse, and generates an incremental strain located in a small surrounding domain, i.e. in an Eshelby inclusion. This local reorganization generates a perturbation of the stress field and can result in several secondary grain reorganizations in the vicinity of the first reorganization.

By solving analytically this problem for Mohr-Coulomb type materials, the geometrical directions corresponding to the largest increase of the shear stress to normal stress ratios are recovered. Such directions are given by a quadrupolar function that predicts slip line directions close to $\pm 45^\circ$. The slip line direction is larger or smaller than $45^\circ$ depending on whether the triggering event is contractive or dilative respectively. In this approach, it is interesting to see that the slip line direction combines stress and strain ingredients.

In Mc Namara et al. [72], the authors simulated the Eshelby inclusion problem from 2D DEM computations for different intergranular friction angles. Whatever the friction used, they observed a slight increase in the slip line direction when the axial strain is increased in a biaxial loading test. Compared with the present results in Fig. 19, the qualitative convex shape of the curve is also observed in Mc Namara et al. [72].

In addition, the analytical work achieved in Mc Namara et al. [72] underlines that the slip line orientation can also be influenced by the underlying microstructure anisotropy, especially for large friction angles. As stated by the authors, the anisotropy effect is responsible for additional increase in the slip line orientations as well as in the fluctuations around the mean value. Such a result is consistent with the observed increase in the fluctuations of the shear chain orientations shown in Fig. 8 for instance.

In this section, we discuss how the meso shear chains evolve according to plasticity theories. Meso shear chains defined in this paper is a meso representation of slips for granular materials, and they are different from the well-known shear band in terms of sizes, duration, etc. The detailed comparison of these two sheared structures will be explained in future work of the authors. La Ragione et al. [73] developed a formulation to describe the shear band orientation, will also be compared in our future papers.
# 6 Conclusion

Based on series of 2D biaxial DEM tests, this paper investigates the development of meso shear structures (the new-defined incremental shear strain chains) in granular materials under different conditions. A quantitative definition of incremental shear strain chain (briefly “shear chain”) is proposed and the detection algorithm is detailed. Based on the kinematic definition of shear chains we have shown that shear chains are spatially correlated with microscopic dissipation and that their orientation is consistent with slip line theory in continuum mechanics. Consequently, we believe that shear chains correspond indeed to slip lines but defined in discrete mechanics instead of continuum mechanics framework. Main points can be summarized as follows:

1. According to the observations of incremental deviatoric strain field of a dense granular assembly under biaxial loading, the concept of incremental shear strain chain is proposed to characterize meso shear patterns (i.e., microband) from a quantitative point of view. The geometric features of meso shear structures are well captured with the definition of shear chain by demonstrating orientation distribution with marked double peaks close to $\pm45^\circ$. Such a rational concept is defined as a dual concept to the well-known “force chain”, to analyze the mesoscopic shear deformation features of granular materials.

2. The orientation of the shear chains in granular materials is shown to be a material-scale quantity which depends essentially on the macro stress ratio $q/p$. This conclusion was reached by considering different boundary conditions (drained or proportional strain tests), different confining pressures and different initial sample densities.

3. The shear chain concept is related to plastic dissipation at contact scale, by demonstrating that the microscopic contact sliding is more concentrated in shear chain connections. In the framework of non-associated elasto-plastic theory, the slip line is considered as internal shear events. The shear chain orientation is consistent with the slip line direction which is influenced by both the yield surface and the plastic flow rule. The more microstructure rearranges, the more the flow rule drives the shear chain direction. In addition, it is conjectured that the shear chain direction is probably influenced by the underlying microstructural anisotropy.

4. In our future work, more loading conditions for granular materials should be considered to enhance our conclusion and a better algorithm for searching shear chains should be developed to incorporate all shear chains. As the shear chain could be related to the micro slip, the comparison between shear chain assumption and mathematical formulation (e.g., Gaspar and Koenders [56], Karimi and Barrat [40], La Ragione et al. [53]) is necessary to further develop the framework of shear chain definition. Besides, shear bands and shear chains are shear behaviors at different scales (time and space), how to correlate them will be presented in our following papers.
Acknowledgements

This work was supported by National Natural Science Foundation of China
under Grant No.51909194 and Zhejiang Provincial Natural Science Foundation
of China under Grant No. LY22E090002. The authors express their sincere
thanks to the International Research Network GeoMech (IRN CNRS) for pro-
moting positive and convivial interactions among researchers. We thank the
editors and reviewers for their helpful suggestions on the quality improvement
of our paper.

Compliance with ethical standards

The authors declare that they have no conflict of interest.

**Appendix A Force chain definition**

In the understanding of the micromechanics of granular materials, grain loops and force chains are important mesostructures. Force chains account for the macroscopic mechanical strength of granular materials, while loops play an important role in volumetric and anisotropic evolutions. Moreover, loops surrounding force chains are known to strongly influence the mechanical stability of granular materials [24, 50, 74–76].

The definition of a force chain [58] <fnref n="15" /> is briefly recalled here and illustrated
in Fig. A1:

- A particle within a force chain fulfils that its major principal stress is larger than the mean major principal stress of all grains ($\sigma_1 \ge \langle\sigma_1\rangle$).

- The major principal stress direction of chained particles is aligned with the geometrical direction of contact (less than 45° deviation).

- A force chain contains at least 3 contacting particles.

Force chain definition relies on three thresholds that involve some arbitrari-
ness. However a minimum of three particles is necessary to define a chained
structure, and two directions are correlated/anti-correlated if they differ by an
angle smaller or larger than 45°. Besides, the threshold ($\sigma_1/ \langle \sigma_1 \rangle$) $\ge$ 1 is
justified by the shape of the distribution function of the principal stresses in
grains (as long as grain polydispersity is not too large) [16].

<footnote n="15">¹⁵Note that force chains, as defined by Peters et al. [58], should be better renamed as stress chained. However we keep the name force chains as it is widely used in the literature.</footnote>
28 Shear chain concept for granular materials

<figure><img src="image_21.png" /><figcaption>Fig. A1 Schematic drawing of a force chain.</figcaption></figure>

# Appendix B Additional details for incremental shear strain chain

## B.1 Sensitivity analysis of parameters for defining shear chains

When defining an incremental shear strain chain in Section 3, 3 parameters should be appointed:

- the minimum number of loops within a shear chain;

- the partition threshold of selected area of large incremental deviatoric strain;

- the deviation angle between the shear direction of strain tensor and the geometrical direction of two connected loop cells, i.e., $\beta_{th}$ in Fig. 5.

For the minimum number of loops within a shear chain, 3 is chosen for the reason that connecting structures containing elements below 3 could not be regarded as chains.

For the selected area for searching shear chains, the average incremental deviatoric strain $<d\epsilon_d^{loop}>$ is used to distinguish the large incremental deviatoric strain. Here we compare the patterns of different partitions of large incremental deviatoric strain for State B of the sample in Section 2, by using 0.5 $<d\epsilon_d^{loop}>$, 0.8 $<d\epsilon_d^{loop}>$, $<d\epsilon_d^{loop}>$, 1.5 $<d\epsilon_d^{loop}>$, 2 $<d\epsilon_d^{loop}>$ respectively in Fig. B2. It can be seen that the smaller threshold does not result in a clear shearing pattern and a too large threshold will ignore some meso shear structures. Considering the shear chain orientation evolution, we also plot Fig. B3, and it can be demonstrated that the average orientation of shear chains are relatively robust with a proper meso strain threshold. In this paper, we use $<d\epsilon_d^{loop}>$.

For the value of $\beta_{th}$, i.e., the deviation angle between the shear direction of strain tensor and the geometrical direction of two connected loop cells in Fig. 5, some skewing should be accepted to gain more shear chains, but $\beta_{th}$ should not be too large when the chains do not correspond to shear features. Here
<figure><img src="image_22.png" /><figcaption>**Fig. B2** Spatial patterns of selected large strain area defined by different thresholds: 0.5 < $d\epsilon_d^{loop}$>, 0.8 <$d\epsilon_d^{loop}$>, <$d\epsilon_d^{loop}$>, 1.5 <$d\epsilon_d^{loop}$>, 2 <$d\epsilon_d^{loop}$>.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. B3** Evolution of average orientation of shear chains with respect to different strain thresholds.</figcaption></figure>

we compared seven values: $\beta_{th} = 10^{\circ}$, $\beta_{th} = 15^{\circ}$, $\beta_{th} = 19^{\circ}$, $\beta_{th} = 22.5^{\circ}$ (chosen value in the paper), $\beta = 25^{\circ}$, $\beta_{th} = 27^{\circ}$, $\beta_{th} = 30^{\circ}$. The shear chains captured by $\beta_{th} = 10^{\circ}$, $\beta_{th} = 22.5^{\circ}$ and $\beta_{th} = 30^{\circ}$ are shown in Fig. B4, compared to selected area of large incremental deviatoric strain. In Fig. B5, the $|\theta|$ evolutions before shear banding are compared for different $\beta_{th}$ thresholds. Different options of $\beta_{th}$ will affect the shear chain orientations. As the meso structure initially orientate 45° in some theories [40, 41] and small differences are found between $\beta_{th} = 19^{\circ}$ and $\beta_{th} = 25^{\circ}$ in Fig. 5(b), the chosen $\beta_{th} = 22.5^{\circ}$ in this paper should be assumed to be rational.

## B.2 Robustness of the shear chain detection algorithm

In a similar way to the force chain detection algorithm [58], the shear chain detection algorithm does not account for possible shear chain branches which may occur if the shear zone appears to be thicker than a grain loop. Since the order in which shear chains are detected is random, two independent runs of the algorithm gives slightly different sets of shear chains. In order to quantify the amount of shear chains that are always detected whatever the run, we choose 5 states of the biaxial test to run 20 times the algorithm, and then
30

Shear chain concept for granular materials

<figure><img src="image_24.png" /><figcaption>Fig. B4 Shear chain distribution captured by different $\beta_{th}$: $\beta_{th} = 10^{\circ}$, $\beta_{th} = 22.5^{\circ}$, $\beta_{th} = 30^{\circ}$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. B5 Evolution of average orientation of shear chains with respect to different $\beta_{th}$: (a) all $\beta_{th}$ values; (b) proper $\beta_{th}$ values.</figcaption></figure>

compare the loop Ids composing the sets of shear chains. Table B1 summarizes the results.

<table><caption>Table B1 Proportions of overlap shear chain loops in 5, 10 or 20 independent runs of the shear chain detection algorithm.</caption><thead><tr><th>States chosen</th><th>5 trials</th><th>10 trials</th><th>20 trials</th></tr></thead><tbody><tr><td>ε<sub>1</sub> = 0.001</td><td>74.5%</td><td>71.6%</td><td>70.9%</td></tr><tr><td>ε<sub>1</sub> = 0.003</td><td>74.1%</td><td>71.3%</td><td>70.6%</td></tr><tr><td>ε<sub>1</sub> = 0.008</td><td>77.6%</td><td>74.3%</td><td>73.9%</td></tr><tr><td>ε<sub>1</sub> = 0.014</td><td>77.2%</td><td>72.6%</td><td>71.3%</td></tr><tr><td>ε<sub>1</sub> = 0.025</td><td>77.5%</td><td>70.5%</td><td>69.9%</td></tr></tbody></table>

As shown in Table B1, the overlap part of loops for 5 trials, 10 trials and 20 trials decreases slightly but we can conclude that around 70% of shear chain loops are always found whatever the particular run of the algorithm. From the authors' experience, such result is quite similar to the one obtained with force chains detection. As a result, we can conclude that shear chain statistics
computed on a single run are relevant and that the current algorithm for the
shear strain chain in Section 3.1 is rational.

# Appendix C Statistical representativeness of the increasing trend in shear chain orientation during biaxial tests

As illustrated in Fig. 8, as the deviatoric stress increases, the orientation of the shear chains exhibit larger fluctuations around its mean value. As a result, one could argue that the evolution of the mean $|\theta|$ angle plotted in Fig. 9 may not be representative as the standard deviation becomes too large.

To identify whether the trend in Fig. 9 is meaningful, we use ANOVA (Analysis of Variance) to check whether the probability distribution of shear chains of a given state is similar to its previous ones. The orientations of the shear chains are regarded as realizations of different random variables $|\theta|(\varepsilon_1)$, and ANOVA test is used to assess which extent these random variables are similar. Several states before the stress peak are chosen for the ANOVA investigation, and each distribution $|\theta|(\varepsilon_1^t)$ is compared with its previous ones $|\theta|(\varepsilon_{yy}^t, t \in [1, i-1])$. After several ANOVA calculations, the corresponding $F$ and $p$ for ANOVA are shown in Fig. C6.

<figure><img src="image_26.png" /><figcaption>**Fig. C6** $F$ (a) and $p$ (b) evolutions for ANOVA of shear chain orientations along biaxial loading.</figcaption></figure>

According to the law of ANOVA, the coefficient $p$ characterize the degree of similarity of the selected distributions. If $p$ is small enough (typically $p < 0.05$), we can estimate that the set of considered distributions are not the same. Fig. 6(b) shows that the probability $p$ drops when the axial strain is larger than 0.006. In other words, the distribution of the shear chain orientation experiences important changes during the deviatoric loading especially when the axial strain is larger than 0.006. From this point, the increasing trend for
$|\theta|$ observed in Fig. 9 is meaningful and confirms the relevance to consider the mean value $<|\theta|>$ in the present paper.

## References

[1] Ucgul, M., Fielke, J. M. & Saunders, C. Three-dimensional discrete element modelling of tillage: Determination of a suitable contact model and parameters for a cohesionless soil. *Biosystems Engineering* **121**, 105–117 (2014).

[2] Zhao, C.-F., Kruyt, N. P. & Millet, O. Capillary bridges between spherical particles under suction control: Rupture distances and capillary forces. *Powder technology* **360**, 622–634 (2020).

[3] Daouadji, A. & Hicher, P.-Y. An enhanced constitutive model for crushable granular materials. *International Journal for Numerical and Analytical Methods in Geomechanics* **34** (6), 555–580 (2010). URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nag.815. https://doi.org/10.1002/nag.815, https://arxiv.org/abs/https://onlinelibrary.wiley.com/doi/pdf/10.1002/nag.815.

[4] Zhou, W. *et al*. Discrete element modeling of particle breakage considering different fragment replacement modes. *Powder Technology* **360**, 312–323 (2020).

[5] Shi, X. & Zhao, J. Practical estimation of compression behavior of clayey/silty sands using equivalent void-ratio concept. *Journal of Geotechnical and Geoenvironmental Engineering* **146** (6), 04020046 (2020).

[6] Shi, X. S., Liu, K. & Yin, J. Analysis of mobilized stress ratio of gap-graded granular materials in direct shear state considering coarse fraction effect. *Acta Geotechnica* (2021).

[7] Xiong, H., Yin, Z. Y., Zhao, J. & Yang, Y. Investigating the effect of flow direction on suffusion and its impacts on gap-graded granular soils. *Acta Geotechnica* **16** (3) (2021).

[8] Zhou, W., Ma, Q., Ma, G., Cao, X. & Cheng, Y. Microscopic investigation of internal erosion in binary mixtures via the coupled lbm-dem method. *Powder Technology* **376** (2020).

[9] Papadopoulos, L., Porter, M. A., Daniels, K. E. & Bassett, D. S. Network analysis of particles and grains. *Journal of Complex Networks* **6** (4), 485–565 (2018). URL https://doi.org/10.1093/comnet/cny005. https://doi.org/10.1093/comnet/cny005, https://arxiv.org/abs/http://oup.prod.sis.lan/comnet/article-pdf/6/4/485/25451976/cny005.pdf.
[10] Walker, D. M. & Tordesillas, A. Topological evolution in dense granular materials: a complex networks perspective. *International Journal of Solids and Structures* **47** (5), 624–639 (2010).

[11] Wan, R. & Pouragha, M. Fabric and connectivity as field descriptors for deformations in granular media. *Continuum Mechanics and Thermodynamics* **27** (1-2), 243–259 (2015).

[12] Liu, J., Zhou, W., Ma, G., Yang, S. & Chang, X. Strong contacts, connectivity and fabric anisotropy in granular materials: a 3d perspective. *Powder Technology* (2020).

[13] Zhao, J. & Guo, N. The interplay between anisotropy and strain localization in granular soils: a multiscale insight. *Géotechnique* **65** (8), 642–656 (2015). URL https://doi.org/10.1680/geot.14.P.184. https://doi.org/10.1680/geot.14.P.184, https://arxiv.org/abs/https://doi.org/10.1680/geot.14.P.184.

[14] Li, X. & Yu, H.-S. On the stress–force–fabric relationship for granular materials. *International Journal of Solids and Structures* **50** (9), 1285–1302 (2013).

[15] Kruyt, N. P. Micromechanical study of fabric evolution in quasi-static deformation of granular materials. *Mechanics of materials* **44**, 120–129 (2012).

[16] Radjai, F., Wolf, D. E., Jean, M. & Moreau, J.-J. Bimodal character of stress transmission in granular packings. *Phys. Rev. Lett.* **80**, 61–64 (1998). URL https://link.aps.org/doi/10.1103/PhysRevLett.80.61. https://doi.org/10.1103/PhysRevLett.80.61.

[17] Guo, N. & Zhao, J. The signature of shear-induced anisotropy in granular media. *Computers and Geotechnics* **47**, 1–15 (2013).

[18] Alshibli, K. A. & El-Saidany, H. A. Quantifying void ratio in granular materials using voronoi tessellation. *Journal of Computing in Civil Engineering* **15** (3), 232–238 (2001).

[19] Cil, M. B. & Alshibli, K. A. 3d analysis of kinematic behavior of granular materials in triaxial testing using dem with flexible membrane boundary. *Acta Geotechnica* **9** (2), 287–298 (2014). URL https://doi.org/10.1007/s11440-013-0273-0. https://doi.org/10.1007/s11440-013-0273-0.

[20] Radjai, F., Roux, J.-N. & Daouadji, A. Modeling granular materials: Century-long research across scales. *Journal of Engineering Mechanics* **143** (4), 04017002 (2017). URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%29EM.1943-7889.0001196. https://doi.org/10.1061/(ASCE)
34 *Shear chain concept for granular materials*
EM.1943-7889.0001196.

[21] Rechenmacher, A. L., Abedi, S., Chupin, O. & Orlando, A. D. Characterization of mesoscale instabilities in localized granular shear using digital image correlation. *Acta Geotechnica* **6** (4), 205–217 (2011). URL https://doi.org/10.1007/s11440-011-0147-2. https://doi.org/10.1007/s11440-011-0147-2.

[22] Nicot, F., Hadda, N. & Darve, F. Second-order work analysis for granular materials using a multiscale approach. *International Journal for Numerical and Analytical Methods in Geomechanics* **37** (17), 2987–3007 (2013). URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nag.2175. https://doi.org/10.1002/nag.2175, https://arxiv.org/abs/https://onlinelibrary.wiley.com/doi/pdf/10.1002/nag.2175.

[23] Zhao, C.-F., Yin, Z.-Y. & Hicher, P.-Y. A multiscale approach for investigating the effect of microstructural instability on global failure in granular materials. *International Journal for Numerical and Analytical Methods in Geomechanics* **42** (17), 2065–2094 (2018).

[24] Wautier, A., Bonelli, S. & Nicot, F. Micro-inertia origin of instabilities in granular materials. *International Journal for Numerical and Analytical Methods in Geomechanics* **42** (9), 1037–1056 (2018).

[25] Xiong, H., Yin, Z. Y., Nicot, F., Wautier, A. & Philippe, P. A novel multiscale large deformation approach for modelling of granular collapse. *Acta Geotechnica* 1–18 (2021).

[26] Rudnicki, J. W. & Rice, J. Conditions for the localization of deformation in pressure-sensitive dilatant materials. *Journal of the Mechanics and Physics of Solids* **23** (6), 371–394 (1975).

[27] Oda, M. & Kazama, H. Microstructure of shear bands and its relation to the mechanisms of dilatancy and failure of dense granular soils. *Géotechnique* **48** (4), 465–481 (1998). URL https://doi.org/10.1680/geot.1998.48.4.465. https://doi.org/10.1680/geot.1998.48.4.465, https://arxiv.org/abs/https://doi.org/10.1680/geot.1998.48.4.465.

[28] Rechenmacher, A. L. Grain-scale processes governing shear band initiation and evolution in sands. *Journal of the Mechanics and Physics of Solids* **54** (1), 22–45 (2006). URL http://www.sciencedirect.com/science/article/pii/S0022509605001481. https://doi.org/10.1016/j.jmps.2005.08.009.

[29] Hasan, A. & Alshibli, K. Experimental assessment of 3d particle-to-particle interaction within sheared sand using synchrotron microtomography. *Géotechnique* **60** (5), 369–379 (2010). URL https://doi.org/
10.1680/geot.2010.60.5.369. https://doi.org/10.1680/geot.2010.60.5.369, https://arxiv.org/abs/https://doi.org/10.1680/geot.2010.60.5.369.

[30] Desrues, J., Andò, E., Mevoli, F. A., Debove, L. & Viggiani, G. How does strain localise in standard triaxial tests on sand: Revisiting the mechanism 20 years on. *Mechanics Research Communications* **92**, 142 – 146 (2018). URL http://www.sciencedirect.com/science/article/pii/S0093641318300247. https://doi.org/https://doi.org/10.1016/j.mechrescom.2018.08.007.

[31] Iwashita, K. & Oda, M. Micro-deformation mechanism of shear banding process based on modified distinct element method. *Powder technology* **109** (1-3), 192–205 (2000).

[32] Gu, X., Huang, M. & Qian, J. Discrete element modeling of shear band in granular materials. *Theoretical and Applied Fracture Mechanics* **72**, 37–49 (2014).

[33] Zhu, H., Nguyen, H. N., Nicot, F. & Darve, F. On a common critical state in localized and diffuse failure modes. *Journal of the Mechanics and Physics of Solids* **95**, 112–131 (2016).

[34] Liu, J., Nicot, F. & Zhou, W. Sustainability of internal structures during shear band forming in 2d granular materials. *Powder Technology* **338**, 458 – 470 (2018). URL http://www.sciencedirect.com/science/article/pii/S003259101830490X. https://doi.org/https://doi.org/10.1016/j.powtec.2018.07.001.

[35] Ma, G., Regueiro, R. A., Zhou, W., Wang, Q. & Liu, J. Role of particle crushing on particle kinematics and shear banding in granular materials. *Acta Geotechnica* **13** (3), 601–618 (2018).

[36] Kuhn, M. R. Structured deformation in granular materials. *Mechanics of Materials* **31** (6), 407 – 429 (1999). URL http://www.sciencedirect.com/science/article/pii/S0167663699000101. https://doi.org/https://doi.org/10.1016/S0167-6636(99)00010-1.

[37] Le Bouil, A., Amon, A., McNamara, S. & Crassous, J. Emergence of cooperativity in plasticity of soft glassy materials. *Phys. Rev. Lett.* **112**, 246001 (2014). URL https://link.aps.org/doi/10.1103/PhysRevLett.112.246001. https://doi.org/10.1103/PhysRevLett.112.246001.

[38] Zhou, S., Tordesillas, A., Pouragha, M., Bailey, J. & Bondell, H. On local intrinsic dimensionality of deformation in complex materials. *arXiv preprint arXiv:2104.01775* (2021).
36 Shear chain concept for granular materials

[39] Tordesillas, A., Muthuswamy, M. & Walsh, S. D. Mesoscale mea-
sures of nonaffine deformation in dense granular assemblies. *Journal of
engineering mechanics* **134** (12), 1095–1113 (2008).

[40] Karimi, K. & Barrat, J.-L. Correlation and shear bands in a plastically deformed granular medium. *Scientific Reports* **8** (1), 4021 (2018).

[41] Houdoux, D., Nguyen, T. B., Amon, A. & Crassous, J. Plastic flow and localization in an amorphous material: experimental interpretation of the fluidity. *Physical Review E* **98** (2), 022905 (2018).

[42] Ibrahimbegovic, A. & Brancherie, D. Combined hardening and softening constitutive model of plasticity: precursor to shear slip line failure. *Computational Mechanics* **31** (1-2), 88–100 (2003).

[43] Hashiguchi, K. & Tsutsumi, S. Shear band formation analysis in soils by the subloading surface model with tangential stress rate effect. *International Journal of Plasticity* **19** (10), 1651–1677 (2003).

[44] Voyiadjis, G. Z., Alsaleh, M. I. & Alshibli, K. A. Evolving internal length scales in plastic strain localization for granular materials. *International journal of plasticity* **21** (10), 2000–2024 (2005).

[45] Zaiser, M. & Aifantis, E. C. Randomness and slip avalanches in gradient plasticity. *International journal of plasticity* **22** (8), 1432–1455 (2006).

[46] Aifantis, E. C. The physics of plastic deformation. *International journal of plasticity* **3** (3), 211–247 (1987).

[47] Ostoja-Starzewski, M. Scale effects in plasticity of random media: sta-
tus and challenges. *International journal of plasticity* **21** (6), 1119–1160
(2005).

[48] Darve, F., Nicot, F., Wautier, A. & Liu, J. Slip lines versus shear bands: Two competing localization modes. *Mechanics Research Communications* 103603 (2020).

[49] Nguyen, N.-S., Magoariec, H. & Cambou, B. Local stress analysis in granular materials at a mesoscale. *International Journal for Numerical and Analytical Methods in Geomechanics* **36** (14), 1609–1635 (2012).

[50] Liu, J., Wautier, A., Bonelli, S., Nicot, F. & Darve, F. Macroscopic softening in granular materials from a mesoscale perspective. *International Journal of Solids and Structures* (2020).

[51] Clerc, A., Wautier, A., Bonelli, S. & Nicot, F. Meso-scale signatures of inertial transitions in granular materials. *Granular Matter* **23** (2), 1–14
(2021).

[52] Shi, Z.-H. & Horii, H. Microslip model of strain localization in sand deformation. *Mechanics of Materials* **8** (2), 89 – 102 (1989). URL http://www.sciencedirect.com/science/article/pii/0167663689900094. https://doi.org/https://doi.org/10.1016/0167-6636(89)90009-4.

[53] La Ragione, L., Prantil, V. C. & Sharma, I. A simplified model for inelastic behavior of an idealized granular material. *International Journal of Plasticity* **24** (1), 168–189 (2008).

[54] Zhu, Q., Shao, J.-F. & Mainguy, M. A micromechanics-based elasto-plastic damage model for granular materials at low confining pressure. *International Journal of Plasticity* **26** (4), 586–602 (2010).

[55] Koenders, M. Evolution of spatially structured elastic materials using a harmonic density function. *Physical Review E* **56** (5), 5585 (1997).

[56] Gaspar, N. & Koenders, M. Micromechanic formulation of macroscopic structures in a granular medium. *Journal of engineering mechanics* **127** (10), 987–993 (2001).

[57] Zhang, J., Behringer, R. P. & Goldhirsch, I. Coarse-graining of a physical granular system. *Progress of Theoretical Physics Supplement* **184**, 16–30 (2010).

[58] Peters, J. F., Muthuswamy, M., Wibowo, J. & Tordesillas, A. Characterization of force chains in granular material. *Phys. Rev. E* **72**, 041307 (2005). URL https://link.aps.org/doi/10.1103/PhysRevE.72.041307. https://doi.org/10.1103/PhysRevE.72.041307.

[59] Kruyt, N. P. & Rothenburg, L. Micromechanical definition of the strain tensor for granular materials. *Journal of Applied Mechanics* **63** (3), 706–711 (1996).

[60] Nguyen, N.-S., Magoariec, H., Vincens, E. & Cambou, B. On the definition of a relevant meso-scale for upscaling the mechanical behavior of 3d granular materials. *Granular Matter* **22** (1), 1–5 (2020).

[61] Šmilauer, V. et al. Yade Documentation 2nd ed (The Yade Project, 2015). Http://yade-dem.org/doc/.

[62] Zhang, L. & Evans, T. M. Boundary effects in discrete element method modeling of undrained cyclic triaxial and simple shear element tests. *Granular Matter* **20** (4), 1–23 (2018).
[63] Artoni, R. & Richard, P. Coarse graining for granular materials: micro-polar balances. *Acta Mechanica* **230** (9), 3055–3069 (2019).

[64] Deng, N. et al. On the attraction power of critical state in granular materials. *Journal of the Mechanics and Physics of Solids* **149**, 104300 (2021). URL https://www.sciencedirect.com/science/article/pii/S0022509621000089. https://doi.org/https://doi.org/10.1016/j.jmps.2021.104300.

[65] Bagi, K. Stress and strain in granular assemblies. *Mechanics of materials* **22** (3), 165–177 (1996).

[66] Dedecker, F., Chaze, M., Dubujet, P. & Cambou, B. Specific features of strain in granular materials. *Mechanics of Cohesive-frictional Materials* **5** (3), 173–193 (2000).

[67] Cambou, B., Jean, M. & Radjai, F. *Micromechanics of granular materials* (John Wiley & Sons, 2013).

[68] Deng, N. et al. Lifespan dynamics of cluster conformations in stationary regimes in granular materials. *Physical Review E* **105** (1), 014902 (2022).

[69] Jrad, M., Sukumaran, B. & Daouadji, A. Experimental analyses of the behaviour of saturated granular materials during axisymmetric proportional strain paths. *European Journal of Environmental and Civil Engineering* **16** (1), 111–120 (2012). URL https://doi.org/10.1080/19648189.2012.666900. https://doi.org/10.1080/19648189.2012.666900, https://arxiv.org/abs/https://doi.org/10.1080/19648189.2012.666900.

[70] Nicot, F., Daouadji, A., Hadda, N., Jrad, M. & Darve, F. Granular media failure along triaxial proportional strain paths. *European Journal of Environmental and Civil Engineering* **17** (9), 777–790 (2013). URL https://doi.org/10.1080/19648189.2013.819301. https://doi.org/10.1080/19648189.2013.819301, https://arxiv.org/abs/https://doi.org/10.1080/19648189.2013.819301.

[71] Magnanimo, V., La Ragione, L., Jenkins, J. T., Wang, P. & Makse, H. A. Characterizing the shear and bulk moduli of an idealized granular material. *EPL (Europhysics Letters)* **81** (3), 34006 (2008).

[72] Mcnamara, S., Crassous, J. & Amon, A. Eshelby inclusions in granular matter: theory and simulations. *Physical Review E* **94** (2), 022907 (2016).

[73] La Ragione, L., Prantil, V. & Jenkins, J. A micromechanical prediction of localization in a granular material. *Journal of the Mechanics and Physics of Solids* **83**, 146–159 (2015).
[74] Tordesillas, A., Walker, D. M. & Lin, Q. Force cycles and force chains. *Physical Review E* **81** (1), 011302 (2010).

[75] Tordesillas, A., Lin, Q., Zhang, J., Behringer, R. & Shi, J. Structural stability and jamming of self-organized cluster conformations in dense granular materials. *Journal of the Mechanics and Physics of Solids* **59** (2), 265–296 (2011).

[76] Zhu, H., Nicot, F. & Darve, F. Meso-structure organization in two-dimensional granular materials along biaxial loading path. *International Journal of Solids and Structures* **96**, 25–37 (2016).