Mesoscopic scale instability in particulate materials

: 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.


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 andDarve 2006, 2007;Nicot et al. 2007bNicot et al. , 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 andDarve 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 elasticplastic 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(Sibille et al. , 2008)), a bifurcation domain is exhibited in the macroscopic force 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.

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 (n, t) can be attached to this configuration, where n is given by the line joining the two centers G 1 and G 2 , and t is normal to 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 u M the displacement of any Point M with respect to the reference frame R (x 1 , x 2 ), it follows that Because where the symbol × stands as the cross product.Eq. ( 1) gives the standard kinematic relation The first term δu G 2 =R − δu G 1 =R corresponds to the incremental change in the vector G 1 G 2 , whereas the second term where α = rotation of the vector G 1 G 2 with respect to the frame (n, 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 where δu n ¼ δl = incremental normal displacement; and δu t ¼ −lδα = incremental tangential displacement.Eq. ( 5) is independent of the choice of the arbitrary Point I, justifying the notation δu.
In discrete element modelling and in absence of particle rotation around their center of mass, the increments of both normal Δu n and tangential Δu t displacements are computed from the current positions of the centers of the spheres (Fig. 2).Noting G 1 G 0 2 ¼ l 0 n, the following expressions can be derived with respect to the frame (ñ, t) For small increments (α ¼ Δα ≪ 1), Eq. ( 6) gives in a firstorder approximation 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 0 in the final configuration, Fig. 3).Both vectors G 2 G 0 2 and II 0 are equal.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 φ g : The standard soil mechanics sign convention is adopted throughout the text: forces are counted positive in compression.
In the elastic regime (jTj < tan φ g N, or jTj ¼ tan φ g N with Ṫ < tan φ g Ṅ), In the plastic regime (jTj ¼ tan φ g N and j Ṫj ¼ tan φ g j Ṅj), where ζ n = sign of un ; and ζ t = sign of ut .
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 δu t ¼ lδα is not integrable.An explicit expression u t ðl; αÞ satisfying the differential equation δu t ¼ lδα 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 (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 ¼ Δ.
Assuming an elastic behavior at the contact between the two spheres, the external energy associated with Path 1 is given by Likewise, the external energy associated with Path 2 is given by When the tangential displacement Δ is not zero, Eqs. ( 9) and (10) show that E 1 ≠ 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.

Constitutive Formulation
This paper considers a regular assembly of four spherical particles in contact, experiencing a diamond pattern as depicted in Fig. 5.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: where α 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 α o .It is also convenient to introduce the kinematic variables Δ 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: Balance Eq. ( 12) can be written under a rate form, relating the kinematic variables ( Δ1 , Δ2 ) to the static variables ( Ḟ1 , Ḟ2 ).Differentiating Eq. ( 12) gives 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 Accounting for Eq. ( 11) yields Thus, as Δ1 ¼ − L1 and Δ2 ¼ − L2 , and recalling that un ¼ l and ut ¼ −l α Combining this with Eq. ( 14) gives and by virtue of Eq. ( 16), it follows that: Fig. 5.The diamond pattern: geometrical model and notations Eqs. ( 17) and ( 19) can also be expressed as Ḟ ¼ K Δ, where K is the tangent stiffness matrix In the elastic regime, which gives , both matrices K and K have the same spectral properties; they are associated with the same endomorphism (Appendix I).Thus, the problem can be investigated by considering equation Finally, the constitutive equations can be expressed as a dimensionless relation By setting s ¼ S=k n and K ¼ k n K, with The only geometrical parameters involved in the tangent stiffness matrix are angles α o and α, with tan The nonlinearity stems from the terms α, 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, K can be split into two terms, K ¼ K mat þ K geo , where K mat represents the material origin of the constitutive behavior, involving only the elastic or the plastic behavior on the contact scale between the spheres K geo accounts for the geometrical origin of the constitutive behavior, with Indeed, unlike K geo , K mat contains the material parameters Λ 1 and Λ 2 .K geo includes (in addition of the stress terms s 1 and s 2 ) the geometrical parameter α.
When the local behavior is plastic, the matrix K mat is clearly nonsymmetric As expected, K mat is symmetric when the local behavior is elastic.
As shown in Eq. ( 28), K geo is nonsymmetric, except in the situation where It can readily be shown that Eq. ( 30) is equivalent to the relation 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 K 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 K mat is symmetric, but the matrix K geo is nonsymmetric 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 ≠ 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 (α 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, whereasthe axial loading is monotonously increased The simulation is run using first the parameters reported in Table 1.
Because a purely isotropic initial state is considered (α o ¼ π=4), a low value was affected to the friction angle (φ 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 ¼ Δ 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 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 ṡ1 , along the drained biaxial path, is given by Thus, at Point A, when the loading path turns from the isotropic compression loading to the drained biaxial loading, det 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 K can be computed analytically.According to Eq. ( 26), the following relation can be derived: Fig. 6.Biaxial test: evolution of the axial stress Taking the definition of Λ 1 and Λ 2 into account, in the elastic regime, one has Likewise, in the plastic regime Restricting the analysis to the elastic case, it appears from Eq. ( 35) that the vanishing of det K occurs when One can particularize this analysis by considering the isotropic compression case.Then, 26) gives Thus, the following can be readily obtained: In these conditions, Eq. ( 38) can be rewritten as K is, therefore, symmetric along the isotropic loading path.
Likewise, by virtue of Eq. ( 39), Eq. ( 37) gives Eq. ( 41) specifies the value E Ã of the strain E at which det K vanishes along an isotropic compression path.At this state The vanishing of det 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 α o ).However, it must be emphasized that for α o ¼ π=4, as considered in the previous numerical simulation (Table 1), the stiffness matrix takes the simplified form During the isotropic compression loading, the strain increases gradually, until reaching the value 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 K Ã degenerates into the spherical matrix 1=2J, J being the unit matrix, as all terms are equal to 1=2.So, strictly speaking, the determinant of K, when the tangent stiffness matrix is not degenerated into a spherical matrix, is never zero.Before the state E Ã is reached, all eigenvalues of K are strictly positive (λ 1 ¼ 1 and λ 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 (λ 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 α o ≠ π=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 α o ¼ π=4.
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 :

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 ( Ė1 , Ė2 ) and (ṡ 1 , ṡ2 ), related through the constitutive Eq. ( 26), as follows 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 ( Ė1 , Ė2 ) are imposed along all the   directions of the strain rate space, and the stress rate response (ṡ 1 , ṡ2 ) is computed.Then, the normalized second-order work is determined For convenience, a polar plot is chosen (Laouafa and Darve 2002), made up of points of coordinates ½cos α 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 j det Kj.When det K is negative, the amplitude of the cone increases (resp.decreases) when j det Kj also increases (resp.decreases).A strong discontinuity takes place at Point B, when the local behavior becomes plastic.In addition, the evolution of α 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, 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(Tordesillas et al. , 2011(Tordesillas et al. , 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 Ḟ1 and Ḟ2 in Eq. ( 13) are composed of two parts.A first part involves the constitutive behavior on the contact scale: terms 2ð Ṅ cos α þ Ṫ sin αÞ and 2ð Ṅ sin α − Ṫ cos αÞ.These terms give rise to the symmetric operator K mat [Eq.( 27)].A second part involves a purely geometrical effect: terms F 2 α and F 1 α.These terms, together with the definition of the tangent incremental displacements given in Eq. ( 5), are directly related to the nonsymmetric operator K geo , see Eq. ( 28).
Eq. ( 12) can be rewritten as where R ¼ h cos α − sin α sin α cos α i = rotation matrix of angle α.Any change in R, directed by a change in the orientation α 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 σ acting within a granular assembly of volume V containing N c contacts is given by where l c = branch vector relating the centers of contacting particles; and f c = inter-particle forces.Ignoring the change in the volume, the evolution of σ is governed by the change in the contact forces f c (depending on the local behavior) and the change in the branch vectors 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 secondorder 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).
Using Eq. ( 48), and taking into account that all matrices R A and R B and their inverses R −1 A and R −1 B are symmetric and commute, Eq. ( 49) can be transformed into and thus where q H = quadratic form associated with matrix H.As Y ¼ R A R B X, where the matrix R A R B is symmetric, definite and positive, the quadratic forms q H and q are equivalent.Matrices K and H, therefore, have the same spectral properties.Moreover, as H ¼ P −1 KP, matrices K and H are similar.They are associated with the same endomorphism and have the same eigenvalues.
In conclusion, both matrices K and 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 K can be split into two terms, K ¼ Kmat þ Kgeo , where Kmat represents the material origin of the constitutive behavior, involving only the elastic or the plastic behavior on the contact scale between the spheres, and Kgeo accounts for the geometrical origin of the constitutive behavior with Λ 1 ¼ k t =k n and Λ 2 ¼ k t =k n in the plastic regime, and Λ 1 ¼ ζ n ζ t ðtan ϕ g = tan αÞ and Λ 2 ¼ −ζ n ζ t tan φ g tan α in the plastic regime 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 and The relation K ¼ K mat þ K geo , together with Eqs. ( 56) and ( 57), gives Eq. ( 25).

Fig. 3 .Fig. 4 .
Fig. 3. Material description of the relative motion between two contacting spheres without proper rotation

Fig. 7 .
Fig. 7. Biaxial test: evolution of the determinant of the tangent stiffness matrix

Fig. 10 .
Fig. 10.Evolution of the eigenvalues of both the tangent stiffness matrix and its symmetric part over the loading path

Table 1 .
Numerical Simulation: Constitutive Parameters and Initial Conditions

Table 2 .
Refined Numerical Simulation: Constitutive Parameters and Initial Conditions