# Mesoscopic scale instability in particulate materials

François Nicot, Guillaume Veylon, Huaxiang Zhu, Jean Lerbet, Félix Darve

► To cite this version:

François Nicot, Guillaume Veylon, Huaxiang Zhu, Jean Lerbet, Félix Darve. Mesoscopic scale in-stability in particulate materials. Journal of Engineering Mechanics - ASCE, 2016, 142 (8), pp.1-12. 10.1061/(ASCE)EM.1943-7889.0001100. hal-01352023

HAL Id: hal-01352023

https://hal.science/hal-01352023

Submitted on 6 Jan 2020

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
# Mesoscopic Scale Instability in Particulate Materials

François Nicot<fnref n="1" />; Guillaume Veylon, Ph.D.<fnref n="2" />; Zhu Huaxiang, Ph.D.<fnref n="3" />; Jean Lerbet<fnref n="4" />; and Félix Darve<fnref n="5" />

**Abstract:** This manuscript investigates some instability features in granular materials by considering an elementary grain arrangement on the intermediate scale. Although force chains have long been recognized as playing a basic role in the strength of granular specimens, the collaborative contribution of grain loops (grain arrangement) has been highlighted more recently. As a result, the stability of grain loops is expected to strongly govern the stability of the whole assembly. This paper shows that such elementary patterns can be destabilized even though the contact law between granules is elastic. This behavior stems from the nature of the kinematical model describing the geometrical interaction between neighboring grains.

**Keywords:** Stress; Instability; Collapse; Failure; Second-order work; Granular materials; Microstructure; Grain loops.

## Introduction

It has been well known for half a century that force chains appear in granular assemblies when loaded (Dantu 1957). These force chains are made up of more or less aligned contacting grains. They carry the largest forces in the direction of the major principal stress. More recently (Tordesillas et al. 2012), other meso-structures have been characterized by considering contacting grains forming closed loops [in two-dimensional (2D) conditions], which have been called "force loops". These force loops can comprise three, four, five, or even more grains in contact, but the number of such loops decreases rapidly as the number of grains involved decreases (Tordesillas and Muthuswamy 2009; Tordesillas et al. 2012; Kruyt 2012; Zhu et al. 2016).

Furthermore, microstructural analyses of stability in granular media (Nicot and Darve 2006, 2007; Nicot et al. 2007b, 2012b) have shown the close relation between macroscopic second-order work and the sum of microscopic discrete second-order works. The first is defined from the macroscopic incremental stress and strain tensors applied to the boundary of the body or of the representative elementary volume at the macro level. The second are obtained from the incremental interaction forces and relative displacements between grains at the micro level. Moreover, several papers have been published showing that the second-order work criterion is a general stability criterion for rate-independent materials, subjected to quasi-static loading conditions leading to divergence instability (Laouafa and Darve 2002; Darve et al. 2004; Nicot and Darve 2006, 2007; Nicot et al. 2007a, b; Daouadji et al. 2011). This paper will not come back on these basic aspects of instability in granular media. The loss of stability by divergence appears to be related

to the loss of definite positiveness of the constitutive matrix (or equivalently of its symmetric part) for homogeneous elastoplastic problems and of the stiffness matrix for (1) elastic discrete systems (Challamel et al. 2008; Nicot et al. 2012a; Lerbet et al. 2013) or (2) elastoplastic boundary value problems simulated by the Finite Element Method (Laouafa et al. 2011).

Focusing on the microstructural scale, the first preliminary analyses (Hadda et al. 2013) of granular instability seem to show failure mechanisms corresponding to force chain breaking by buckling attributable to the microstructural failure of force loops laterally supporting the force chains. This scenario is currently being fully confirmed (see for example, Sibille et al. 2015) by discrete element computations with the discrete element method (DEM) code called "Yade" (Smilauer et al. 2010), demonstrating that force loop and force chain instability plays a fundamental role in granular failure.

Therefore, it is useful to independently investigate the stability of force loops and force chains as the basic meso-structures, as observed in granular media. This is precisely the purpose of this paper: to propose such stability analyses for a four-grain loop referred to as "the diamond pattern". Indeed, taking into account the fact that three grain loops appear in DEM modelling as a basically stable loop in a granular assembly, the four-grain loop may be viewed as the simplest loop that may be unstable depending on the previously discussed second-order work criterion. As is usual in discrete element methods, the interaction forces between grains are based on elastic linear stiffnesses in the normal and tangential directions with respect to the intergranular contact plane. Moreover, the tangential force is limited by a coulombian friction, leading to an elastic-plastic behavior in the tangential direction.

In this context, it will be shown that an explicit stiffness matrix can be exhibited in both cases of the purely elastic and the elastic-plastic diamond meso-structures. The stiffness matrix, whose expression is analytically established, is not symmetric even in the case of a purely elastic intergranular behavior. This nonsymmetry appears to result from geometrical effects in the case of solely elastic components, and it is induced both by material and geometrical aspects in the case of elastic-plastic interaction forces (see also Kuhn and Chang 2006). Thus, because of this nonsymmetry, a stability analysis using the second-order work criterion is pertinent. As obtained for large granular assemblies (Sibille et al. 2007, 2008), a bifurcation domain is exhibited in the macroscopic force

<footnote n="1">¹Professor, IRSTEA, Geomechanics Group, ETNA, Grenoble F38402, France (corresponding author). E-mail: francois.nicot@irstea.fr</footnote>

<footnote n="2">²IRSTEA, OHAX, Aix-en-Provence, France.</footnote>

<footnote n="3">³IRSTEA, Geomechanics Group, ETNA, Grenoble F38402, France.</footnote>

<footnote n="4">⁴Professor, IBISC, Université d'Evry Val d'Essonne, Evry, France.</footnote>

<footnote n="5">⁵Professor, UJF-INPG-CNRS, Laboratoire Sols Solides Structures Risques, Grenoble, France.</footnote>
plane with a lower bound given by the curve related to the vanishing values of the determinant of the symmetric part of the stiffness matrix. In this bifurcation domain, instability cones are observed which correspond to the “isotropic cones” of linear algebra theory for nonpositive definite matrices. A linear perturbation approach which is classical in a continuum mechanics framework (e.g., Benallal and Comi 2003) is not proposed here, but rather a direct geometric instability analysis to properly take into account the essential discrete nature of a granular force loop.

Consequently, the intriguing macroscopic granular instabilities appear to be present at the mesoscale in the force loops observed experimentally in 2D and numerically by DEM. Thus this paper will show that these meso-structures constitute the elemental place where instability arises, in fine leading to macroscopic failure.

## Contact Model between Spheres

### Constitutive Formulation

The classical description of sphere interaction assumes that the spheres are perfectly rigid and can overlap, resulting in contact forces. Let us consider two spheres, $S_1$ and $S_2$, with the same radius $r_g$, in contact (with a possible overlapping). A local frame $(\mathbf{n}, \mathbf{t})$ can be attached to this configuration, where $\mathbf{n}$ is given by the line joining the two centers $G_1$ and $G_2$, and $\mathbf{t}$ is normal to $\mathbf{n}$ (Fig. 1).

The relative motion of $S_2$ with respect to $S_1$ can be described by analyzing the incremental displacement of a given Point I belonging to the interfacial region between the two spheres. By denoting $\mathbf{u}_M$ the displacement of any Point M with respect to the reference frame R ($x_1, x_2$), it follows that

$$ \delta \mathbf{u}_I = \delta \mathbf{u}_{I \in S_2 / S_1} = \delta \mathbf{u}_{I \in S_2 / R} - \delta \mathbf{u}_{I \in S_1 / R} \quad (1) $$

Because

$$ \delta \mathbf{u}_{I \in S_2 / R} = \delta \mathbf{u}_{G_2 / R} + \mathbf{IG}_2 \times \delta \mathbf{\omega}_2 \quad (2a) $$

$$ \delta \mathbf{u}_{I \in S_1 / R} = \delta \mathbf{u}_{G_1 / R} + \mathbf{IG}_1 \times \delta \mathbf{\omega}_1 \quad (2b) $$

where the symbol $\times$ stands as the cross product. Eq. (1) gives the standard kinematic relation

$$ \delta \mathbf{u}_I = \delta \mathbf{u}_{G_2 / R} - \delta \mathbf{u}_{G_1 / R} + \mathbf{IG}_2 \times \delta \mathbf{\omega}_2 - \mathbf{IG}_1 \times \delta \mathbf{\omega}_1 \quad (3) $$

The first term $\delta \mathbf{u}_{G_2 / R} - \delta \mathbf{u}_{G_1 / R}$ corresponds to the incremental change in the vector $\mathbf{G}_1 \mathbf{G}_2$, whereas the second term

<figure><img src="image_2.png" /><figcaption>Fig. 1. Contact between two spheres: geometrical setting</figcaption></figure>

$\mathbf{IG}_2 \times \delta \mathbf{\omega}_2 - \mathbf{IG}_1 \times \delta \mathbf{\omega}_1$ takes the proper rotation of spheres around their center of mass into account. Noting $\mathbf{G}_1 \mathbf{G}_2 = l \mathbf{n}$, it follows that

$$ \delta(\mathbf{G}_1 \mathbf{G}_2) = \delta l \mathbf{n} - l \delta \alpha \mathbf{t} \quad (4) $$

where $\alpha = \text{rotation of the vector } \mathbf{G}_1 \mathbf{G}_2$ with respect to the frame $(\mathbf{n}, \mathbf{t})$ attached to the sphere $S_1$, resulting from the motion of $S_2$ with respect to $S_1$ (Fig. 2).

If the authors ignore hereafter the proper rotation of spheres around their center of mass, the relative incremental displacement of $S_2$ with respect to $S_1$, therefore, reads

$$ \delta \mathbf{u} = \delta l \mathbf{n} - l \delta \alpha \mathbf{t} \quad (5) $$

where $\delta u_n = \delta l = \text{incremental normal displacement}$; and $\delta u_t = -l \delta \alpha = \text{incremental tangential displacement}$. Eq. (5) is independent of the choice of the arbitrary Point I, justifying the notation $\delta \mathbf{u}$.

In discrete element modelling and in absence of particle rotation around their center of mass, the increments of both normal $\Delta u_n$ and tangential $\Delta u_t$ displacements are computed from the current positions of the centers of the spheres (Fig. 2). Noting $\mathbf{G}_1 \mathbf{G}_2' = l' \mathbf{n}$, the following expressions can be derived with respect to the frame $(\tilde{n}, \tilde{t})$

$$ \Delta u_n = l' \cos \alpha - l \quad (6a) $$

$$ \Delta u_t = -l' \sin \alpha \quad (6b) $$

For small increments ($\alpha = \Delta \alpha \ll 1$), Eq. (6) gives in a first-order approximation

$$ \Delta u_n \approx l' - l = \Delta l \quad (7a) $$

$$ \Delta u_t \approx -l \Delta \alpha \quad (7b) $$

which corresponds to the incremental formulation in Eq. (5). The kinematical description developed previously, leading to Eq. (5), is, therefore, thought to be general enough and representative of the broad class of the discrete element method (enabling particle overlapping).

It should be stressed that this kinematical analysis is based on the displacement of a material point (Point I in the initial configuration, Fig. 3). This point is attached to sphere $S_2$, and is moved with sphere $S_2$ with respect to sphere $S_1$ (Point $I'$ in the final configuration, Fig. 3). Both vectors $\mathbf{G}_2 \mathbf{G}_2'$ and $\mathbf{II}'$ are equal.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Relative motion between two contacting spheres: geometrical model and notations</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Material description of the relative motion between two contacting spheres without proper rotation</figcaption></figure>

A contact law can, therefore, be built by relating a contact force to the displacement of this material point. More precisely, this contact law can be built by relating the normal (N) and tangential (T) components of the contact force acting between the two spheres to the normal ($u_n$) and tangential ($u_t$) components of the relative displacement.

The most popular model for describing the local behavior is an elastic-plastic mechanical model, including a Mohr-Coulomb criterion. This model can be expressed under the following rate formalism, which introduces a normal elastic stiffness $k_n$ and a tangential elastic stiffness $k_t$, both constant, and a local friction angle $\varphi_g$:

$$ \dot{N} = -k_n \dot{u}_n = -k_n \dot{l} \quad (N = 0 \text{ if } l > 2r_g) \quad (8a) $$

The standard soil mechanics sign convention is adopted throughout the text: forces are counted positive in compression.

In the elastic regime ($|T| < \tan \varphi_g N$, or $|T| = \tan \varphi_g N$ with $\dot{T} < \tan \varphi_g \dot{N}$),

$$ \dot{T} = -k_t \dot{u}_t = k_t l \dot{\alpha} \quad (8b) $$

In the plastic regime ($|T| = \tan \varphi_g N$ and $|\dot{T}| = \tan \varphi_g |\dot{N}|$),

$$ \dot{T} = -\zeta_n \zeta_t \tan \varphi_g k_n \dot{u}_n = -\zeta_n \zeta_t \tan \varphi_g k_n \dot{l} \quad (8c) $$

where $\zeta_n$ = sign of $\dot{u}_n$; and $\zeta_t$ = sign of $\dot{u}_t$.

In absence of the proper rotation of spheres around their center of mass, this model is the one commonly used in most numerical discrete element models.

Coming back to the definition of the $u_n$ and $u_t$, $u_n$ is clearly an exact differential, which is not the case for $u_t$. The relation $\delta u_t = l\delta\alpha$ is not integrable. An explicit expression $u_t(l, \alpha)$ satisfying the differential equation $\delta u_t = l\delta\alpha$ cannot be derived.

The absence of the exact character of the differential $u_t$ is at the origin of some peculiar features. In the elastic regime, Eqs. (8a) and (8b) hold. These equations are associated with a nonconservative behavior. The behavior is clearly reversible (there is no dissipation). However, the external energy necessary for the system to evolve from a given state to another state depends on the path followed. Eqs. (8a) and (8b) (no sliding occurs along the tangential direction of contact), therefore, model a hypoelastic behavior (Truesdell 1955; Ericksen 1958; Bernstein 1960).

To illustrate the dependence of the external energy with respect to the loading path, the following example can be considered (Fig. 4). Two spheres (with the same radius) are initially in contact

<figure><img src="image_5.png" /><figcaption>Fig. 4. Illustration of the path dependency of the contact law: geometrical model and notations</figcaption></figure>

(Configuration 1), with $OA = l$. Sphere 1 is fixed, whereas Sphere 2 is moved to a final position defined by the position C of its center of mass (Configuration 2). Two loading paths can be considered:

* Path 1 (ABC) is composed of a purely shearing motion (AB), followed by a purely normal compression (BC).

* Path 2 (ADC) is composed of a purely normal compression (AD), followed by a purely shearing motion (DC), with $BC = AD = \Delta$.

Assuming an elastic behavior at the contact between the two spheres, the external energy associated with Path 1 is given by

$$ E_1 = (k_t l^2 \alpha^2 + k_n \Delta^2)/2 \quad (9) $$

Likewise, the external energy associated with Path 2 is given by

$$ E_2 = [k_n \Delta^2 + k_t (l - \Delta)^2 \alpha^2]/2 \quad (10) $$

When the tangential displacement $\Delta$ is not zero, Eqs. (9) and (10) show that $E_1 \neq E_2$. Thus, the external work required to transform the system from the initial Configuration 1 to the final Configuration 2 depends on the loading path followed, even though no internal dissipation is introduced in the system. This very basic example illustrates a hypoelastic behavior.

The purpose of the next section is to examine some constitutive consequences of this local nonconservativity, on the scale of a small particle arrangement. In particular, it will be shown that a nonsymmetric tangent stiffness matrix can be derived, even in the elastic regime. As reviewed in the introduction, this nonsymmetry, stemming from the nonconservativity of the contact law, will induce the existence of a proper bifurcation domain, enabling the occurrence of various failure modes.

## Diamond Pattern

### Constitutive Formulation

This paper considers a regular assembly of four spherical particles in contact, experiencing a diamond pattern as depicted in Fig. 5.
<figure><img src="image_6.png" /><figcaption>Fig. 5. The diamond pattern: geometrical model and notations</figcaption></figure>

The spheres have the same radius $r$. The elementary assembly is assumed to be loaded by a system of two centered perpendicular forces, aligned along the principal loading directions: an axial force, $F_1$, and a lateral force, $F_2$. The geometry of the assembly is described by indicates of the two parameters $L_1$ and $L_2$. For symmetry reasons, the distance between the centers of any two adjoining spheres is the same and is denoted $l$. At any stage of the loading, the following geometrical relations hold:

$$L_1 = 2l \cos \alpha \qquad (11a)$$

$$L_2 = 2l \sin \alpha \qquad (11b)$$

where $\alpha$ denotes the angle between the straight line joining the centroids of any spheres in contact (for example, centroids $G_1$ and $G_2$ in Fig. 5).

Before any loading, $l_o = 2r$. The initial configuration is, therefore, properly described by the angle $\alpha_o$. It is also convenient to introduce the kinematic variables $\Delta_1 = L_{1o} - L_1$ and $\Delta_2 = L_{2o} - L_2$, where $L_{1o} = 2l_o \cos \alpha_o$ and $L_{2o} = 2l_o \sin \alpha_o$.

The symmetry of the problem ensures that no torque is applied to the particles. Hence, the motion of each particle is only composed of a translational component. The particles have no rotation. Applying an external loading ($F_1$, $F_2$) makes the geometry of the assembly change. The change in geometry is related to the relative motion between the spheres. The normal component of the relative displacement is denoted $u_n$, whereas the tangential counterpart is denoted $u_t$. This relative displacement directs a contact force between adjoining grains, where $N$ = normal component; and $T$ = tangential component. At any stage of the loading, the following balance equations hold:

$$F_1 = 2(N \cos \alpha + T \sin \alpha) \qquad (12a)$$

$$F_2 = 2(N \sin \alpha - T \cos \alpha) \qquad (12b)$$

Balance Eq. (12) can be written under a rate form, relating the kinematic variables ($\dot{\Delta}_1$, $\dot{\Delta}_2$) to the static variables ($\ddot{F}_1$, $\ddot{F}_2$). Differentiating Eq. (12) gives

$$\dot{F}_1 = 2(\dot{N} \cos \alpha + \dot{T} \sin \alpha) - F_2 \dot{\alpha} \qquad (13a)$$

$$\dot{F}_2 = 2(\dot{N} \sin \alpha - \dot{T} \cos \alpha) + F_1 \dot{\alpha} \qquad (13b)$$

The local behavior, described properly using the elastic-plastic mechanical model presented in the previous section, can be introduced in Eq. (13), by distinguishing the elastic and plastic regimes. Elastic regime

$$\dot{F}_1 = 2(-k_n \cos \alpha \dot{u}_n - k_t \sin \alpha \dot{u}_t) - F_2 \dot{\alpha} \qquad (14a)$$

$$\dot{F}_2 = 2(-k_n \sin \alpha \dot{u}_n + k_t \cos \alpha \dot{u}_t) + F_1 \dot{\alpha} \qquad (14b)$$

Accounting for Eq. (11) yields

$$\dot{L}_1 = 2l \cos \alpha - 2l \sin \alpha \dot{\alpha} \qquad (15a)$$

$$\dot{L}_2 = 2l \sin \alpha + 2l \cos \alpha \dot{\alpha} \qquad (15b)$$

Thus, as $\dot{\Delta}_1 = -\dot{L}_1$ and $\dot{\Delta}_2 = -\dot{L}_2$, and recalling that $\dot{u}_n = l$ and $\dot{u}_t = -l\dot{\alpha}$

$$2\dot{u}_n = -\cos \alpha \dot{\Delta}_1 - \sin \alpha \dot{\Delta}_2 \qquad (16a)$$

$$2\dot{u}_t = -\sin \alpha \dot{\Delta}_1 + \cos \alpha \dot{\Delta}_2 \qquad (16b)$$

Combining this with Eq. (14) gives

$$\begin{aligned} \dot{F}_1 = & \left( k_n \cos^2 \alpha + k_t \sin^2 \alpha - F_2 \frac{\sin \alpha}{\sqrt{L_1^2 + L_2^2}} \right) \dot{\Delta}_1 \\ & + \left[ (k_n - k_t) \cos \alpha \sin \alpha + F_2 \frac{\cos \alpha}{\sqrt{L_1^2 + L_2^2}} \right] \dot{\Delta}_2 \end{aligned} \qquad (17a)$$

$$\begin{aligned} \dot{F}_2 = & \left[ (k_n - k_t) \cos \alpha \sin \alpha + F_1 \frac{\sin \alpha}{\sqrt{L_1^2 + L_2^2}} \right] \dot{\Delta}_1 \\ & + \left( k_t \cos^2 \alpha + k_n \sin^2 \alpha - F_1 \frac{\cos \alpha}{\sqrt{L_1^2 + L_2^2}} \right) \dot{\Delta}_2 \end{aligned} \qquad (17b)$$

where $\cos \alpha = L_1/\sqrt{L_1^2 + L_2^2}$ and $\sin \alpha = L_2/\sqrt{L_1^2 + L_2^2}$. Elastoplastic regime:

$$\dot{F}_1 = 2(-k_n \cos \alpha \dot{u}_n - \zeta_n \zeta_t \sin \alpha \tan \varphi_g k_n \dot{u}_n) - F_2 \dot{\alpha} \qquad (18a)$$

$$\dot{F}_2 = 2(\zeta_n \zeta_t \cos \alpha \tan \varphi_g k_n \dot{u}_n - k_n \sin \alpha \dot{u}_n) + F_1 \dot{\alpha} \qquad (18b)$$

and by virtue of Eq. (16), it follows that:

$$\begin{aligned} \dot{F}_1 = & \left[ k_n \cos \alpha (\cos \alpha + \zeta_n \zeta_t \sin \alpha \tan \varphi_g) - F_2 \frac{\sin \alpha}{\sqrt{L_1^2 + L_2^2}} \right] \dot{\Delta}_1 \\ & + \left[ k_n \sin \alpha (\cos \alpha + \zeta_n \zeta_t \sin \alpha \tan \varphi_g) \right. \\ & \left. \qquad + F_2 \frac{\cos \alpha}{\sqrt{L_1^2 + L_2^2}} \right] \dot{\Delta}_2 \end{aligned} \qquad (19a)$$
$$ \dot{F}_2 = \left[ k_n \cos \alpha (-\zeta_n \zeta_t \cos \alpha \tan \varphi_g + \sin \alpha) + F_1 \frac{\sin \alpha}{\sqrt{L_1^2 + L_2^2}} \right] \dot{\Delta}_1 + \left[ k_n \sin \alpha (-\zeta_n \zeta_t \cos \alpha \tan \varphi_g + \sin \alpha) - F_1 \frac{\cos \alpha}{\sqrt{L_1^2 + L_2^2}} \right] \dot{\Delta}_2 \quad (19b) $$

Eqs. (17) and (19) can also be expressed as $\dot{\mathbf{F}} = \mathbf{K}\dot{\mathbf{\Delta}}$, where $\mathbf{K}$ is the tangent stiffness matrix

$$ \tilde{\mathbf{K}} = k_n \left[ \begin{array}{cc} \cos^2\alpha + \Lambda_1\sin^2\alpha - \sin\alpha \frac{F_2}{k_n\sqrt{L_1^2 + L_2^2}} & (1-\Lambda_2)\cos\alpha\sin\alpha + \cos\alpha \frac{F_2}{k_n\sqrt{L_1^2 + L_2^2}} \\ (1-\Lambda_1)\cos\alpha\sin\alpha + \sin\alpha \frac{F_1}{k_n\sqrt{L_1^2 + L_2^2}} & \sin^2\alpha + \Lambda_2\cos^2\alpha - \cos\alpha \frac{F_1}{k_n\sqrt{L_1^2 + L_2^2}} \end{array} \right] \quad (20) $$

In the elastic regime, $\Lambda_1 = k_t/k_n$ and $\Lambda_2 = k_t/k_n$

In the plastic regime, $\Lambda_1 = \zeta_n\zeta_t(\tan\phi_g/\tan\alpha)$ and $\Lambda_2 = -\zeta_n\zeta_t(\tan\phi_g/\tan\alpha)$

Noting $\mathbf{S} = \begin{bmatrix} S_1 \\ S_2 \end{bmatrix}$ and $\mathbf{E} = \begin{bmatrix} E_1 \\ E_2 \end{bmatrix}$, with $S_1 = F_1/L_{2o}$, $S_2 = F_2/L_{1o}$, $E_1 = (L_{1o} - L_1)/L_{1o}$, and $E_2 = (L_{2o} - L_2)/L_{2o}$, the constitutive equation $\dot{\mathbf{F}} = \tilde{\mathbf{K}}\dot{\mathbf{\Delta}}$ can be expressed in a Lagrangian formalism as $\dot{\mathbf{S}} = \hat{\mathbf{K}}\dot{\mathbf{E}}$. Then

$$ L_{2o} \dot{S}_1 = \tilde{K}_{11} L_{1o} \dot{E}_1 + \tilde{K}_{12} L_{2o} \dot{E}_2 \quad \text{and} \quad L_{1o} \dot{S}_2 = \tilde{K}_{21} L_{1o} \dot{E}_1 + \tilde{K}_{22} L_{2o} \dot{E}_2 \quad (21) $$

Thus

$$ \hat{K}_{1j} = \tilde{K}_{1j} L_{jo} / L_{2o} \quad \text{and} \quad \hat{K}_{2j} = \tilde{K}_{2j} L_{jo} / L_{1o} \quad (\text{with no summation on repeated indices}) \quad (22) $$

which gives

$$ \hat{\mathbf{K}} = \mathbf{A}\tilde{\mathbf{K}}\mathbf{B} \quad (23) $$

where $\mathbf{A} = \begin{bmatrix} 1/\tan\alpha_o & 0 \\ 0 & 1 \end{bmatrix}$, $\mathbf{B} = \begin{bmatrix} 1 & 0 \\ 0 & \tan\alpha_o \end{bmatrix}$ and $\tan\alpha_o = L_{2o}/L_{1o}$.

As $\hat{\mathbf{K}} = \begin{bmatrix} 1/\tan\alpha_o & 0 \\ 0 & 1 \end{bmatrix} \tilde{\mathbf{K}} \begin{bmatrix} 1 & 0 \\ 0 & \tan\alpha_o \end{bmatrix}$, both matrices $\hat{\mathbf{K}}$ and $\tilde{\mathbf{K}}$ have the same spectral properties; they are associated with the same endomorphism (Appendix I). Thus, the problem can be investigated by considering equation $\dot{\mathbf{F}} = \tilde{\mathbf{K}}\dot{\mathbf{\Delta}}$ or equation $\dot{\mathbf{S}} = \hat{\mathbf{K}}\dot{\mathbf{E}}$ independently, with

$$ \hat{\mathbf{K}} = \cos^2\alpha k_n \left[ \begin{array}{cc} \left(1 + \Lambda_1 k_t \tan^2\alpha - \tan\alpha \frac{S_2/k_n}{1-E_1}\right) / \tan\alpha_o & (1-\Lambda_2 k_t) \tan\alpha + \frac{S_2/k_n}{1-E_1} \\ (1-\Lambda_1 k_t) \tan\alpha + \tan^2\alpha \frac{S_1/k_n}{1-E_2} & \left(\tan^2\alpha + \Lambda_2 k_t - \tan\alpha \frac{S_1/k_n}{1-E_2}\right) \tan\alpha_o \end{array} \right] \quad (24) $$

Finally, the constitutive equations can be expressed as a dimensionless relation

$$ \dot{s} = \mathbf{K}\dot{\mathbf{E}} \quad (25) $$

By setting $\mathbf{s} = \mathbf{S}/k_n$ and $\mathbf{K} = k_n\hat{\mathbf{K}}$, with

$$ \mathbf{K} = \cos^2\alpha \left[ \begin{array}{cc} \left(1 + \Lambda_1 \tan^2\alpha - \tan\alpha \frac{s_2}{1-E_1}\right) / \tan\alpha_o & (1-\Lambda_2)\tan\alpha + \frac{s_2}{1-E_1} \\ (1-\Lambda_1)\tan\alpha + \tan^2\alpha \frac{s_1}{1-E_2} & \left(\tan^2\alpha + \Lambda_2 - \tan\alpha \frac{s_1}{1-E_2}\right)\tan\alpha_o \end{array} \right] \quad (26) $$

The only geometrical parameters involved in the tangent stiffness matrix are angles $\alpha_o$ and $\alpha$, with $\tan\alpha_o = L_{2o}/L_{1o}$ and $\tan\alpha = \tan\alpha_o(1-E_2)/(1-E_1)$. The nonlinearity stems from the terms $\alpha$, $s_k$, and $E_k$ in the components $K_{ij}$.

## Material and Geometry Effects

In both situations (elastic regime and plastic regime), Eq. (26) reveals that the tangent stiffness matrix is nonsymmetric. However, the origin of the nonsymmetry is not the same in the elastic and

plastic regimes. Indeed, $\mathbf{K}$ can be split into two terms, $\mathbf{K} = \mathbf{K}^{\text{mat}} + \mathbf{K}^{\text{geo}}$, where $\mathbf{K}^{\text{mat}}$ represents the material origin of the constitutive behavior, involving only the elastic or the plastic behavior on the contact scale between the spheres

$$ \mathbf{K}^{\text{mat}} = \cos^2\alpha \left[ \begin{array}{cc} (1 + \Lambda_1 \tan^2\alpha)/\tan\alpha_o & (1-\Lambda_2)\tan\alpha \\ (1-\Lambda_1)\tan\alpha & (\tan^2\alpha + \Lambda_2)\tan\alpha_o \end{array} \right] \quad (27) $$

$\mathbf{K}^{\text{geo}}$ accounts for the geometrical origin of the constitutive behavior, with
$$ \mathbf{K}^{\text{geo}} = \cos^2\alpha \left[ \begin{matrix} -\frac{\tan \alpha}{\tan \alpha_o} \frac{s_2}{1-E_1} & \frac{s_2}{1-E_1} \\ \tan^2\alpha \frac{s_1}{1-E_2} & -\tan \alpha_o \tan \alpha \frac{s_1}{1-E_2} \end{matrix} \right] \quad (28) $$

Indeed, unlike $\mathbf{K}^{\text{geo}}$, $\mathbf{K}^{\text{mat}}$ contains the material parameters $\Lambda_1$ and $\Lambda_2$. $\mathbf{K}^{\text{geo}}$ includes (in addition of the stress terms $s_1$ and $s_2$) the geometrical parameter $\alpha$.

When the local behavior is plastic, the matrix $\mathbf{K}^{\text{mat}}$ is clearly nonsymmetric

$$ K_{12}^{\text{mat}} - K_{21}^{\text{mat}} = \zeta_n \zeta_t \tan \varphi_g \quad (29) $$

As expected, $\mathbf{K}^{\text{mat}}$ is symmetric when the local behavior is elastic.

As shown in Eq. (28), $\mathbf{K}^{\text{geo}}$ is nonsymmetric, except in the situation where

$$ \tan^2\alpha \frac{s_1}{1-E_2} = \frac{s_2}{1-E_1} \quad (30) $$

It can readily be shown that Eq. (30) is equivalent to the relation

$$ \cos \alpha F_2 = \sin \alpha F_1 \quad (31) $$

By virtue of Eqs. (12), Eq. (31) leads to $T = 0$. Thus, when the tangential force is zero between contacting particles (no shearing takes place), the matrix $\mathbf{K}^{\text{geo}}$ is symmetric in the elastic regime.

Thus, the nonsymmetry of the tangent stiffness matrix has a double origin: material and geometric. When the local behavior is elastic, the matrix $\mathbf{K}^{\text{mat}}$ is symmetric, but the matrix $\mathbf{K}^{\text{geo}}$ is non-symmetric as soon as the tangential force $T$ is nonzero, making the whole tangent stiffness matrix nonsymmetric.

This is the consequence of the nature of the contact law along the tangential direction. When the tangential behavior is activated ($T \neq 0$), the nonconservativeness of the contact law along the tangential direction makes the tangent stiffness matrix nonsymmetric. This property arises in discrete materials, made up of an assembly of contacting or connecting bodies. Because relative motions can take place between the constitutive bodies, such nonconservative geometrical effects may occur. This result, well established in structural mechanics (Challamel et al. 2008; Nicot et al. 2012a, b; Lerbet et al. 2013), is shown to also be valid in any granular material.

To check this result numerically, the response of the system is investigated along a biaxial loading path in the next section.

## Numerical Investigation

Starting from an unloaded initial geometrical configuration ($\alpha_o$) in which particles are barely touching, an isotropic compression loading ($s_1 = s_2 = s_o$) is first applied. Then, the lateral force is kept constant, whereas the axial loading is monotonously increased

$$ \dot{s}_2 = 0 \quad (32a) $$

$$ \dot{E}_1 = \text{const (positive)} \quad (32b) $$

The simulation is run using first the parameters reported in Table 1.

Because a purely isotropic initial state is considered ($\alpha_o = \pi/4$), a low value was affected to the friction angle ($\varphi_g = 5$ deg) to activate the plastic regime during the biaxial loading. The evolution of the normalized axial stress ($s_1$) as a function of the axial strain ($E_1 = \Delta_1/L_{1o}$) is reported in Fig. 6. The curve is composed of three parts. The first part, up to Point A, corresponds to the isotropic

<table><caption>**Table 1.** Numerical Simulation: Constitutive Parameters and Initial Conditions</caption><thead><tr><th>Parameters (unit)</th><th>Values</th></tr></thead><tbody><tr><td>α<sub>o</sub> (degrees)</td><td>45</td></tr><tr><td>s<sub>o</sub></td><td>0.3</td></tr><tr><td>k<sub>t</sub>/k<sub>n</sub></td><td>0.2</td></tr><tr><td>φ<sub>g</sub> (degrees)</td><td>5</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Biaxial test: evolution of the axial stress</figcaption></figure>

confinement. Both the second part, from Point A to Point B, and the third path beyond Point B correspond to the biaxial loading. The lateral stress $S_2$ is constant ($s_2 = 0.3$). During the second part, the local behavior is elastic, and the nonlinearity character of the response stems from the change in the geometry. Then, after Point B, the local behavior is plastic. Immediately after Point A, the axial stress monotonously decreases (softening regime). Even though the decrease is more pronounced in the plastic regime, it starts in the elastic regime. This is also confirmed by the analysis of the (normalized) determinant of the tangent stiffness matrix, as reported in Fig. 7.

As expected, the determinant is first positive during the confining stage until Point I. Then it becomes and stays negative (end of confining stage, and biaxial loading stage), whatever the nature of the local constitutive regime. This clearly proves that during the confining stage, the material becomes unstable (from Point I). The discontinuity that appears at Point B corresponds to the transition from an elastic regime to an elastoplastic regime, inducing a discontinuous change in the components of the stiffness matrix.

It must be recalled that the axial stress rate $\dot{s}_1$, along the drained biaxial path, is given by

$$ \dot{s}_1 = \frac{K_{11}K_{22} - K_{12}K_{21}}{K_{22}} \dot{E}_1 = \frac{\det \mathbf{K}}{K_{22}} \dot{E}_1 \quad (33) $$

Thus, at Point A, when the loading path turns from the isotropic compression loading to the drained biaxial loading, $\det \mathbf{K} < 0$. Therefore, according to Eq. (33), the slope of the curve $E_1 - s_1$ is negative at Point A, indicating that a softening regime takes place.

More interestingly, $\det \mathbf{K}$ can be computed analytically. According to Eq. (26), the following relation can be derived:
$$ \det \mathbf{K} = \sin^2 \Lambda_1 + \cos^2 \Lambda_2 - \cos \alpha \sin \alpha \left( \frac{s_2}{1-E_1} + \frac{s_1}{1-E_2} \right) + (\Lambda_1 - \Lambda_2) \cos \alpha \sin \alpha \left( \cos^2 \alpha \frac{s_2}{1-E_1} - \sin^2 \alpha \frac{s_1}{1-E_2} \right) \quad (34) $$

Taking the definition of $\Lambda_1$ and $\Lambda_2$ into account, in the elastic regime, one has

$$ \det \mathbf{K} = k_t/k_n - \cos \alpha \sin \alpha \left( \frac{s_2}{1-E_1} + \frac{s_1}{1-E_2} \right) \quad (35) $$

Likewise, in the plastic regime

$$ \det \mathbf{K} = \zeta_n \zeta_t \tan \varphi_g \left( \cos^2 \alpha \frac{s_2}{1-E_1} - \sin^2 \alpha \frac{s_1}{1-E_2} \right) - \cos \alpha \sin \alpha \left( \frac{s_2}{1-E_1} + \frac{s_1}{1-E_2} \right) \quad (36) $$

Restricting the analysis to the elastic case, it appears from Eq. (35) that the vanishing of $\det \mathbf{K}$ occurs when

$$ \cos \alpha \sin \alpha \left( \frac{s_2}{1-E_1} + \frac{s_1}{1-E_2} \right) = k_t/k_n \quad (37) $$

One can particularize this analysis by considering the isotropic compression case. Then, $E_1 = E_2 = E$. Noting that $\alpha = \alpha_o$, Eq. (26) gives

$$ \mathbf{K} = \cos^2 \alpha_o \begin{bmatrix} \left(1/\tan^2 \alpha_o + k_t/k_n\right) \tan \alpha_o - \frac{s_2}{1-E} & \left(1-k_t/k_n\right) \tan \alpha_o + \frac{s_2}{1-E} \\ \left(1-k_t/k_n\right) \tan \alpha_o + \tan^2 \alpha_o \frac{s_1}{1-E} & \left(\tan^2 \alpha_o + k_t/k_n\right) \tan \alpha_o - \tan^2 \alpha_o \frac{s_1}{1-E} \end{bmatrix} \quad (38) $$

Thus, the following can be readily obtained:

$$ \dot{s}_1 = \dot{E} / \tan \alpha_o \quad \text{or} \quad s_1 = E / \tan \alpha_o \quad (39a) $$

$$ \dot{s}_2 = \tan \alpha_o \dot{E} \quad \text{or} \quad s_2 = \tan \alpha_o E \quad (39b) $$

In these conditions, Eq. (38) can be rewritten as

$$ \mathbf{K} = \cos^2 \alpha_o \begin{bmatrix} \left(1/\tan^2 \alpha_o + k_t/k_n - \frac{E}{1-E}\right) \tan \alpha_o & \left(1-k_t/k_n + \frac{E}{1-E}\right) \tan \alpha_o \\ \left(1-k_t/k_n + \frac{E}{1-E}\right) \tan \alpha_o & \left(\tan^2 \alpha_o + k_t/k_n - \frac{E}{1-E}\right) \tan \alpha_o \end{bmatrix} \quad (40) $$

K is, therefore, symmetric along the isotropic loading path.

Likewise, by virtue of Eq. (39), Eq. (37) gives

$$ \frac{E_*}{1 - E_*} = k_t/k_n \quad \text{or} \quad E_* = \frac{k_t/k_n}{1 + k_t/k_n} \quad (41) $$

Eq. (41) specifies the value $E_*$ of the strain $E$ at which $\det \mathbf{K}$ vanishes along an isotropic compression path. At this state

$$ \mathbf{K}_* = \cos \alpha_o \sin \alpha_o \begin{bmatrix} 1/\tan^2 \alpha_o & 1 \\ 1 & \tan^2 \alpha_o \end{bmatrix} \quad (42) $$

The vanishing of $\det \mathbf{K}$ can be obtained even at small strains if the ratio $k_t/k_n$ is small with respect to 1. For example, if $k_t/k_n = 0.5$, $E_* = 0.33$; the vanishing of the determinant is observed at relatively large strains. On the contrary, if $k_t/k_n = 0.2$, $E_* = 0.17$ (as it is perfectly verified in Fig. 7), the vanishing of the determinant occurs at moderate strains. Remarkably, these values are independent of the initial geometrical configuration (angle $\alpha_o$). However, it must be emphasized that for $\alpha_o = \pi/4$, as considered in the previous numerical simulation (Table 1), the stiffness matrix takes the simplified form

$$ \mathbf{K} = 1/2 \begin{bmatrix} 1 + k_t/k_n - \frac{E}{1-E} & 1 - k_t/k_n + \frac{E}{1-E} \\ 1 - k_t/k_n + \frac{E}{1-E} & 1 + k_t/k_n - \frac{E}{1-E} \end{bmatrix} \quad (43) $$

During the isotropic compression loading, the strain increases gradually, until reaching the value $E_* = (k_t/k_n)/(1+k_t/k_n)$ at which $\det \mathbf{K}$ vanishes. As $\det \mathbf{K} = k_t/k_n - [E/(1-E)]$, $\det \mathbf{K}$ decreases from $k_t/k_n$ to zero, and finally takes negative values for strains larger than $E_*$. Thus, although not intuitive, the occurrence of instability (as depicted by the softening regime associated with the drained biaxial loading path) requires a sufficient level of deformation ($E_*$) to be reached at the end of the isotropic confining process. In other words, a sufficient level of normal overlapping (critical normal overlapping) is required at the contacts between grains. If this critical overlapping is not reached, no softening will appear during the drained biaxial loading in the elastic regime.

When $E = E_*$, the stiffness matrix $\mathbf{K}_*$ degenerates into the spherical matrix $1/2\mathbf{J}$, $\mathbf{J}$ being the unit matrix, as all terms are equal to $1/2$. So, strictly speaking, the determinant of $\mathbf{K}$, when the tangent stiffness matrix is not degenerated into a spherical matrix,
<figure><img src="image_8.png" /><figcaption>Fig. 7. Biaxial test: evolution of the determinant of the tangent stiffness matrix</figcaption></figure>

<table><caption>**Table 2.** Refined Numerical Simulation: Constitutive Parameters and Initial Conditions</caption><thead><tr><th>Parameters (unit)</th><th>Values</th></tr></thead><tbody><tr><td>α<sub>o</sub> (degrees)</td><td>35</td></tr><tr><td>s<sub>o</sub></td><td>0.35</td></tr><tr><td>k<sub>t</sub>/k<sub>n</sub></td><td>0.24</td></tr><tr><td>φ<sub>g</sub> (degrees)</td><td>14</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Biaxial test: evolution of the axial stress</figcaption></figure>

is never zero. Before the state $E_*$ is reached, all eigenvalues of **K** are strictly positive ($\lambda_1 = 1$ and $\lambda_2 = k_t/k_n - E/(1-E)$ with $E/(1-E) < k_t/k_n$), as $\det K$ is. For strains larger than $E_*$, **K** admits one negative eigenvalue ($\lambda_2 = k_t/k_n - E/(1-E)$ with $E/(1-E) > k_t/k_n$), making the determinant negative as well.

As shown in Eq. (42), this peculiar situation cannot occur when $\alpha_o \neq \pi/4$. The stiffness matrix $K_*$ does not degenerate into a spherical matrix and admits two eigenvalues, one being zero. This case is considered in the following numerical simulation, using the parameters reported in Table 2.

The evolution of the normalized axial stress as a function of the axial strain is reported in Fig. 8. As for the purely isotropic case, the curve is composed of three parts. The first part, up to Point A, corresponds to the isotropic confinement. After Point A, the axial stress monotonously decreases (softening regime), even in the elastic regime. The softening regime is less pronounced than in the purely isotropic case $\alpha_o = \pi/4$.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Biaxial test: evolution of the determinant of the tangent stiffness matrix</figcaption></figure>

Similar to what was obtained in the previous purely isotropic case, the determinant of the tangent stiffness matrix, as seen in Fig. 9, is first positive during the confining stage until Point I. Then it becomes, and remains, negative (end of the confining stage and biaxial loading stage), whatever the nature of the local constitutive regime. As expected, the vanishing of the determinant occurs at a strain state $E_*$, which depends only on the ratio $k_t/k_n$: $E_* = (k_t/k_n)/(1 + k_t/k_n) = 0.19$.

## Stability Analysis

Because the tangent stiffness operator **K** can be nonsymmetric, even when the local behavior is elastic, material instability can occur. Such instabilities are properly detected by the vanishing of the second-order work (Hill 1958; Darve et al. 2004; Nicot and Darve 2007; Nicot et al. 2012a). The second-order work involves the loading parameters ($\dot{E}_1$, $\dot{E}_2$) and ($\dot{s}_1$, $\dot{s}_2$), related through the constitutive Eq. (26), as follows

$$W_2 = \dot{s}_1 \dot{E}_1 + \dot{s}_2 \dot{E}_2 \quad (44)$$

The mechanical state considered is reputed to be unstable when at least one loading direction exists along which the second-order work takes a negative value. Considering the numerical simulation performed in the previous section (Table 2), a directional analysis is carried out at three different states: points I, A, and B (Figs. 8 and 9). Kinematic probes ($\dot{E}_1$, $\dot{E}_2$) are imposed along all the
<figure><img src="image_11.png" /><figcaption>Fig. 10. Evolution of the eigenvalues of both the tangent stiffness matrix and its symmetric part over the loading path</figcaption></figure>

directions of the strain rate space, and the stress rate response ($\dot{s}_1, \dot{s}_2$) is computed. Then, the normalized second-order work is determined

$$w_2 = \frac{\dot{s}_1 \dot{E}_1 + \dot{s}_2 \dot{E}_2}{\|\dot{s}\| \|\dot{E}\|} \quad (45)$$

For convenience, a polar plot is chosen (Laouafa and Darve 2002), made up of points of coordinates [$\cos \alpha_E(r + w_2)$, $\sin \alpha_E(r + w_2)$], where $\tan \alpha_E = \dot{E}_1 / \dot{E}_2$ and $r = 2$. When the second-order work is negative along a given direction, the corresponding point of the diagram is located inside the circle of radius $r = 2$.

In fact, the second-order work is a quadratic form associated with the symmetric part $K_s$ of $K$ (Nicot and Darve 2009). According to the Bromwich theorem (Ishaq 1955; Iordache and Willam 1998), stating that the smallest eigenvalue of the symmetric part $M_s$ of any square matrix $M$ is lower than any real part of the eigenvalues of $M$ (the inequality is strict when $M$ is nonsymmetric), it follows that the determinant of $M_s$ always vanishes before that of $M$. Hence, when the determinant of $K$ is zero, at least one eigenvalue of $K_s$ is strictly negative. As shown in the previous section, $K$ is symmetric until Point A. After Point A, the eigenvalues of $K$ and $K_s$ are very close and both determinants ($\det K$ and $\det K_s$) as well, until the plastic regime is activated, as seen in Fig. 10. Thus, until Point B, the sign of the second-order work is directly related to the sign of the eigenvalues of $K$.

As shown in Fig. 11, when $\det K$ vanishes (Point I), the instability cone reduces to a single direction. When the loading is pursued (Points A and B), $\det K$ takes strictly negative values. The instability cone opens. The evolution of the instability cone opening over the loading path can be analyzed in Fig. 12. The instability cone opening is represented by hatched lines in Fig. 12. The amplitude of the cones follows the evolution of $|\det K|$. When $\det K$ is negative, the amplitude of the cone increases (resp. decreases) when $|\det K|$ also increases (resp. decreases). A strong discontinuity takes place at Point B, when the local behavior becomes plastic. In addition, the evolution of $\alpha_E$ over the loading path is reported in Fig. 12. It can be observed that the loading direction is not included within the instability cones, after Point I, until Point A is reached. Then, at Point A, the loading path is changed (from an isotropic compression to a biaxial loading with a constant lateral force), and the loading direction belongs to the instability cone. As a result,

<figure><img src="image_12.png" /><figcaption>Fig. 11. Polar diagram of the normalized second-order work. The different diagrams correspond to the different loading stages (Points I, A and B—see Figs. 8 and 9)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Evolution of the instability cone over the loading path</figcaption></figure>
the axial stress $s_1$ decreases after Point A (softening regime), as shown in Fig. 8.

In conclusion, even though the local behavior between particles is purely elastic, an unstable behavior can be observed on the specimen scale. Instability appears long before the plastic criterion is met on the local (contact) scale. It is clearly material on the specimen scale, detected by the second-order work criterion. The local origin of this instability is related to the contact law used to describe the evolution of the contact forces as a function of the relative displacements between spheres. Because of the hypoelastic character of the law before the plastic criterion is met, instability may arise, detected by the vanishing of the second-order work.

## Closure Discussion

The diamond pattern considered in the previous section can be regarded as one basic generic arrangement taking place within a 2D granular assembly. It is now well recognized that a granular material, once loaded, is composed of two distinct, interacting phases: a strong phase is constituted by the force chains gathering along quasi-linear patterns, with the contacting grains transmitting substantial contact forces. This was confirmed by several authors using discrete element simulations (Radjai et al. 1999) or considering experimental tests (Oda and Konishi 1974). These force chains, whose stability is ensured by the surrounding weak phase, are known to be responsible for the strength of the medium. Correlatively, the collapse of these chains directs a limitation of the strength of the medium, as observed at the peak (or near the plateau) of the drained triaxial test (Tordesillas et al. 2010, 2011, 2012). The weak phase is made up of the rest of the medium, mixing different grain loops in 2D. Unlike force chains, the deformability of the grain loops is important attributable to the relative motion between grains. Therefore, as a first approximation, it is believed that the deformability of the granular assembly is ensured by the grain loops within the weak phase, whereas the strength of the medium stems from the force chains. These two phases act in close collaboration, the weak phase (grain loops) directing a stabilization effect for the force chains. This confining effect is even more efficient because the grain loops constituting the weak phase are stable.

The diamond-like pattern considered in this manuscript represents a grain loop. The same approach can be carried out with other grain loops containing more grains (pentagonal or hexagonal patterns), even though the occurrence frequency of grain cycles within a granular assembly decreases when the number of grains increases (Kruyt 2012). Basically, the same conclusions would have been drawn. According to the loading conditions applied, the axial force reaches a peak and follows a descending branch. At the peak, the system is not able to sustain a higher axial force. If an additional axial force is applied to the system, an abrupt collapse is expected to occur, stemming from an imbalance between the external force and the internal contact forces between grains. It is worth emphasizing that this evolution occurs even though the local behavior between grains is purely elastic, shedding light on the geometrical effects related to the local kinematical description at contacts (section 2).

These geometrical effects are particularly visible in both Eqs. (12) and (13). The expressions of $\tilde{F}_1$ and $\tilde{F}_2$ in Eq. (13) are composed of two parts. A first part involves the constitutive behavior on the contact scale: terms $2(\dot{N} \cos \alpha + \dot{T} \sin \alpha)$ and $2(\dot{N} \sin \alpha - \dot{T} \cos \alpha)$. These terms give rise to the symmetric operator $\mathbf{K}^{\text{mat}}$ [Eq. (27)]. A second part involves a purely geometrical effect: terms $F_2\dot{\alpha}$ and $F_1\dot{\alpha}$. These terms, together with the

definition of the tangent incremental displacements given in Eq. (5), are directly related to the nonsymmetric operator $\mathbf{K}^{\text{geo}}$, see Eq. (28).

Eq. (12) can be rewritten as

$$ \begin{bmatrix} F_1 \\ F_2 \end{bmatrix} = 2 \begin{bmatrix} \cos \alpha & -\sin \alpha \\ \sin \alpha & \cos \alpha \end{bmatrix} \begin{bmatrix} N \\ -T \end{bmatrix} \quad (46) $$

where $\mathbf{R} = \begin{bmatrix} \cos \alpha & -\sin \alpha \\ \sin \alpha & \cos \alpha \end{bmatrix} = \text{rotation matrix of angle } \alpha$. Any change in $\mathbf{R}$, directed by a change in the orientation $\alpha$ of the branch vectors joining the centroids of the spheres, results in a change in the external forces ($F_1, F_2$). Thus, the evolution of ($F_1, F_2$) also stems from the change in the micro-structure. This feature can be generalized from the stress averaging Love-Weber formula (Love 1927; Weber 1966; Christoffersen et al. 1981; Mehrabadi et al. 1982). The average stress $\sigma$ acting within a granular assembly of volume $V$ containing $N_c$ contacts is given by

$$ \sigma_{ij} = \frac{1}{V} \sum_{c=1}^{N_c} f_i^c l_j^c \quad (47) $$

where $\mathbf{l}^c = \text{branch vector relating the centers of contacting particles; and } \mathbf{f}^c = \text{inter-particle forces. Ignoring the change in the volume, the evolution of } \sigma$ is governed by the change in the contact forces $\mathbf{f}^c$ (depending on the local behavior) and the change in the branch vectors $\mathbf{l}^c$ (the distribution of the branch vectors relates to the micro-structure). This is a generalization of the feature investigated in this manuscript on the elementary pattern scale.

The emergence of a nonconservative behavior on the grain arrangement scale, characterized by the nonsymmetry of the tangent stiffness operator, makes the diamond pattern a heuristic example of a hypoelastic structure. This hypoelastic behavior observed on the grain arrangement scale derives from the mathematical structure of the kinematical model used to describe the interaction between the particles. It should again be stressed that this model corresponds to that widely used in the discrete element models, in which particles can overlap.

The nonsymmetry of the tangent stiffness operator was shown to give rise to a proper bifurcation domain (even in the elastic regime). This domain gathers all the mechanical states where the second-order work can take negative values along given loading directions. Thus, for the mechanical states within the bifurcation domain, instability modes can occur according to the loading conditions (this corresponds, for example, to the softening regime observed in Figs. 6 and 8). These findings are perfectly in accordance with the recent investigation on elastic structures (Ziegler column), in which it is shown that some loading conditions can destabilize an elastic structure, especially when kinematical constraints are prescribed (Challamel et al. 2008; Nicot et al. 2012a; Lerbet et al. 2013).

## Appendix I. Spectral Properties of Matrices $\hat{\mathbf{K}}$ and $\tilde{\mathbf{K}}$

Matrices $\tilde{\mathbf{K}}$ and $\hat{\mathbf{K}}$ are related by the equation

$$ \hat{\mathbf{K}} = \mathbf{A}\tilde{\mathbf{K}}\mathbf{B} \quad (48) $$

where $\mathbf{A} = \begin{bmatrix} 1/\tan \alpha_o & 0 \\ 0 & 1 \end{bmatrix}$, $\mathbf{B} = \begin{bmatrix} 1 & 0 \\ 0 & \tan \alpha_o \end{bmatrix}$ and $\tan \alpha_o = L_{2o}/L_{1o}$.

The authors observe that $\mathbf{A} = \mathbf{R}_A\mathbf{R}_A$ and $\mathbf{B} = \mathbf{R}_B\mathbf{R}_B$, where matrices $\mathbf{R}_A$ and $\mathbf{R}_B$ are diagonal and, therefore, commute $(\mathbf{R}_A\mathbf{R}_B = \mathbf{R}_B\mathbf{R}_A)$. The quadratic form $\hat{q}$ associated with matrix $\hat{\mathbf{K}}$ reads
$$ \hat{q}(\mathbf{X}) = {}'\mathbf{X}\hat{\mathbf{K}}\mathbf{X} \qquad (49) $$

Using Eq. (48), and taking into account that all matrices $\mathbf{R}_A$ and $\mathbf{R}_B$ and their inverses $\mathbf{R}_A^{-1}$ and $\mathbf{R}_B^{-1}$ are symmetric and commute, Eq. (49) can be transformed into

$$ \hat{q}(\mathbf{X}) = {}'\mathbf{X}'\mathbf{R}_B'\mathbf{R}_A\mathbf{R}_B^{-1}\mathbf{R}_A\tilde{\mathbf{K}}\mathbf{R}_A^{-1}\mathbf{R}_B\mathbf{R}_A\mathbf{R}_B\mathbf{X} \quad (50) $$

and thus

$$ \hat{q}(\mathbf{X}) = {}'(\mathbf{R}_A\mathbf{R}_B\mathbf{X})(\mathbf{R}_A^{-1}\mathbf{R}_B)^{-1}\tilde{\mathbf{K}}(\mathbf{R}_A^{-1}\mathbf{R}_B)(\mathbf{R}_A\mathbf{R}_B\mathbf{X}) \quad (51) $$

Noting $\mathbf{Y} = \mathbf{R}_A\mathbf{R}_B\mathbf{X}$, $\mathbf{P} = \mathbf{R}_A^{-1}\mathbf{R}_B$ and $\mathbf{H} = \mathbf{P}^{-1}\tilde{\mathbf{K}}\mathbf{P}$ yields

$$ \hat{q}(\mathbf{X}) = {}'\mathbf{Y}\mathbf{H}\mathbf{Y} = q_H(\mathbf{Y}) \quad (52) $$

where $q_H$ = quadratic form associated with matrix $\mathbf{H}$.

As $\mathbf{Y} = \mathbf{R}_A\mathbf{R}_B\mathbf{X}$, where the matrix $\mathbf{R}_A\mathbf{R}_B$ is symmetric, definite and positive, the quadratic forms $q_H$ and $\hat{q}$ are equivalent. Matrices $\tilde{\mathbf{K}}$ and $\mathbf{H}$, therefore, have the same spectral properties. Moreover, as $\mathbf{H} = \mathbf{P}^{-1}\tilde{\mathbf{K}}\mathbf{P}$, matrices $\tilde{\mathbf{K}}$ and $\mathbf{H}$ are similar. They are associated with the same endomorphism and have the same eigenvalues.

In conclusion, both matrices $\tilde{\mathbf{K}}$ and $\mathbf{K}$ have the same spectral properties; they are associated with the same endomorphism.

## Appendix II. Dimensionless Constitutive Relation

In both situations (elastic regime and plastic regime), the tangent stiffness matrix $\tilde{\mathbf{K}}$ can be split into two terms, $\tilde{\mathbf{K}} = \tilde{\mathbf{K}}^{\text{mat}} + \tilde{\mathbf{K}}^{\text{geo}}$, where $\tilde{\mathbf{K}}^{\text{mat}}$ represents the material origin of the constitutive behavior, involving only the elastic or the plastic behavior on the contact scale between the spheres, and $\tilde{\mathbf{K}}^{\text{geo}}$ accounts for the geometrical origin of the constitutive behavior

$$ \tilde{\mathbf{K}}^{\text{mat}} = k_n \begin{bmatrix} \cos^2\alpha + \Lambda_1\sin^2\alpha & (1-\Lambda_2)\cos\alpha\sin\alpha \\ (1-\Lambda_1)\cos\alpha\sin\alpha & \sin^2\alpha + \Lambda_2\cos^2\alpha \end{bmatrix} \quad (53) $$

with $\Lambda_1 = k_t/k_n$ and $\Lambda_2 = k_t/k_n$ in the plastic regime, and $\Lambda_1 = \zeta_n\zeta_t(\tan\phi_g/\tan\alpha)$ and $\Lambda_2 = -\zeta_n\zeta_t\tan\phi_g\tan\alpha$ in the plastic regime

$$ \tilde{\mathbf{K}}^{\text{geo}} = \frac{1}{\sqrt{L_1^2 + L_2^2}} \begin{bmatrix} -F_2 \sin \alpha & F_2 \cos \alpha \\ F_1 \sin \alpha & -F_1 \cos \alpha \end{bmatrix} \quad (54) $$

Formally, noting $S_1 = F_1/L_{2o}$, $S_2 = F_2/L_{1o}$, $E_1 = (L_{1o} - L_1)/L_{1o}$ and $E_2 = (L_{2o} - L_2)/L_{2o}$, the governing Eqs. (17) and (19) can be expressed as

$$ \dot{\mathbf{S}} = (\hat{\mathbf{K}}^{\text{mat}} + \hat{\mathbf{K}}^{\text{geo}})\dot{\mathbf{E}} \quad (55) $$

with $\mathbf{S} = \begin{bmatrix} S_1 \\ S_2 \end{bmatrix}$ and $\mathbf{E} = \begin{bmatrix} E_1 \\ E_2 \end{bmatrix}$

$$ \hat{\mathbf{K}}^{\text{mat}} = \cos^2\alpha k_n \begin{bmatrix} (1 + \Lambda_1 k_t \tan^2\alpha) / \tan\alpha_o & (1 - \Lambda_2 k_t) \tan\alpha \\ (1 - \Lambda_1 k_t) \tan\alpha & (\tan^2\alpha + \Lambda_2 k_t) \tan\alpha_o \end{bmatrix} \quad (56) $$

and

$$ \tilde{\mathbf{K}}^{\text{geo}} = \cos^2\alpha k_n \left[ \begin{array}{cc} -\tan\alpha \frac{S_2}{1-E_1} & \tan\alpha_o \\ \tan^2\alpha \frac{S_1/k_n}{1-E_2} & -\tan\alpha \frac{S_1}{1-E_2} \tan\alpha_o \end{array} \right] \quad (57) $$

The relation $\mathbf{K} = \mathbf{K}^{\text{mat}} + \mathbf{K}^{\text{geo}}$, together with Eqs. (56) and (57), gives Eq. (25).

## Acknowledgments

The authors would like to express their sincere thanks to the French Research Network MeGe (Multiscale and multi-physics couplings in geoenvironmental mechanics GDR CNRS 3176/2340, 2008-2015) for having funded this work.

## References

Benallal, A., and Comi, C. (2003). "Perturbation growth and localization in fluid-saturated inelastic porous media under quasi-static loadings." *J. Mech. Phys. Solids*, 51(5), 851-899.

Bernstein, B. (1960). "Hypoelasticity and elasticity." *Arch. Ration. Mech. Anal.*, 6(1), 89-104.

Challamel, N., Nicot, F., Lerbet, J., and Darve, F. (2008). "On the stability of non-conservative elastic systems under mixed perturbations." *Eur. J. Environ. Civ. Eng.*, 13(3), 347-367.

Christoffersen, J., Mehrabadi, M. M., and Nemat-Nasser, S. (1981). "A micro-mechanical description of granular material behavior." *J. Appl. Mech.*, 48(2), 339-344.

Dantu, P. (1957). "Contribution à l'étude mécanique et géométrique des milieux pulvérulents." *Proc., 4th Int. Conf. Soils Mechanical Foundation Engineering*, London, 144-148.

Daouadji, A., et al. (2011). "Diffuse failure in geomaterials, experiments, theory and modelling." *Int. J. Numer. Anal. Methods Geomech.*, 35(16), 1731-1773.

Darve, F., Servant, G., Laouafa, F., and Khoa, H. D. V. (2004). "Failure in geomaterials, continuous and discrete analyses." *Comp. Methods Appl. Mech. Eng.*, 193(27-29), 3057-3085.

Ericksen, J. L. (1958). "Hypoelastic potentials." *Q. J. Mech. Appl. Math.*, 11(1), 67-72.

Hadda, H., Nicot, F., Bourrier, F., Sibille, L., Radjai, F., and Darve, F. (2013). "Micromechanical analysis of second order work in granular media." *Granular Matter*, 15(2), 221-235.

Hill, R. (1958). "A general theory of uniqueness and stability in elastic-plastic solids." *J. Mech. Phys. Solids*, 6(3), 236-249.

Iordache, M. M., and Willam, K. J. (1998). "Localized failure analysis in elastoplastic Cosserat continua." *Comput. Methods Appl. Mech. Eng.*, 151(3-4), 559-586.

Ishaq, M. (1955). "Sur les spectres des matrices." *Séminaire Dubreuil, Al-gèbre et théorie des nombres*, 9, 1-14 (in French).

Kruyt, N. P. (2012). "Micromechanical study of fabric evolution in quasi-static deformation of granular materials." *Mech. Mater.*, 44, 120-129.

Kuhn, M. R., and Chang, C. S. (2006). "Stability, bifurcation, and softening in discrete systems: A conceptual approach for granular materials." *Int. J. Solids Struct.*, 43(20), 6026-6051.

Laouafa, F., and Darve, F. (2002). "Modelling of slope failure by a material instability mechanism." *Comput. Geotech.*, 29(4), 301-325.

Laouafa, F., Prunier, F., Daouadji, A., Al Gali, H., and Darve, F. (2011). "Stability in geomechanics, experimental and numerical analyses." *Int. J. Numer. Anal. Methods Geomech.*, 35(2), 112-139.

Lerbet, J., Kirillov, O., Aldowaji, M., Challamel, N., Nicot, F., and Darve, F. (2013). "Additional constraints may soften a non-conservative structural system: Buckling and vibration analysis." *Int. J. Solids Struct.*, 50(2), 363-370.
Love, A. E. H. (1927). *A treatise of mathematical theory of elasticity*, Cambridge University Press, Cambridge, U.K.

Mehrabadi, M. M., Oda, M., and Nemat-Nasser, S. (1982). "On statistical description of stress and fabric in granular materials." *Int. J. Numer. Anal. Methods Geomech.*, 6(1), 95-108.

Nicot, F., Challamel, N., Lerbet, J., Prunier, F., and Darve, F. (2012a). "Some insights into structure instability and the second order work criterion." *Int. J. Solids Struct.*, 49(1), 132-142.

Nicot, F., and Darve, F. (2006). "Micro-mechanical investigation of material instability in granular assemblies." *Int. J. Solids Struct.*, 43(11-12), 3569-3595.

Nicot, F., and Darve, F. (2007). "A micro-mechanical investigation of bifurcation in granular materials." *Int. J. Solids Struct.*, 44(20), 6630-6652.

Nicot, F., Darve, F., and Khoa, H. D. V. (2007a). "Bifurcation and second-order work in geomaterials." *Int. J. Numer. Anal. Methods Geomech.*, 31(8), 1007-1032.

Nicot, F., Hadda, N., Bourrier, F., Sibille, L., Wan, R., and Darve, F. (2012b). "Inertia effects as a possible missing link between micro and macro second-order work in granular media." *Int. J. Solids Struct.*, 49(10), 1252-1258.

Nicot, F., Sibille, L., Donzé, F., and Darve, F. (2007b). "From microscopic to macroscopic second-order works in granular assemblies." *Mech. Mater.*, 39(7), 664-684.

Oda, M., and Konishi, J. (1974). "Microscopic deformation mechanism of granular material in simple shear." *Soils Found.*, 14(4), 15-32.

Radjai, F., Roux, S., and Moreau, J. J. (1999). "Contact forces in a granular packing." *Chaos*, 9(3), 544-550.

Sibille, L., Donzé, F., Nicot, F., Chareyre, B., and Darve, F. (2008). "Bifurcation detection and catastrophic failure." *Acta Geotecnica*, 3(1), 14-24.

Sibille, L., Nicot, F., Donze, F. V., and Darve, F. (2007). "Material instabilities in granular assemblies from fundamentally different models." *Int. J. Numer. Anal. Methods Geomech.*, 31(3), 457-481.

Sibille, L., Tordesillas, A., Nicot, F., Hadda, N., and Darve, F. (2015). "Granular plasticity, a contribution from discrete mechanics." *J. Mech. Phys. Solids*, 75, 119-139.

Šmilauer, V., et al. (2010). "Yade reference documentation." (http://yade-dem.org/doc/).

Tordesillas, A., Lin, Q., Zhang, J., Behringer, R. P., and Shi, J. (2011). "Structural stability and jamming of self-organized cluster conformations in dense granular materials." *J. Mech. Phys. Solids*, 59(2), 265-296.

Tordesillas, A., and Muthuswamy, M. (2009). "On the modeling of confined buckling of force chains." *J. Mech. Phys. Solids*, 57(4), 706-727.

Tordesillas, A., Pucilowski, S., Sibille, L., Nicot, F., and Darve, F. (2012). "Multiscale characterization of diffuse granular failure." *Philos. Mag.*, 92(36), 4547-4587.

Tordesillas, A., Walker, D. M., and Lin, Q. (2010). "Force cycles and force chains." *Phys. Rev. E*, 81(1), 011302.

Truesdell, C. (1955). "Hypoelasticity." *J. Ration. Mech. Anal.*, 4, 83-133.

Weber, J. (1966). "Recherches concernant les contraintes intergranulaires dans les milieux pulvérulents." *Bulletin de Liaison des Ponts-et-Chaussées*, 20, 1-20.

Zhu, H., Nicot, F., and Darve, F. (2016). "Meso-structure evolution in a 2D granular material during biaxial loading." *Granular Matter*, in press.