Antoine Wautier
email: antoine.wautier@irstea.fr

Stéphane Bonelli

François Nicot

Micro-inertia origin of instabilities in granular materials

Keywords: Stability and bifurcation, Second order work, Force chain, Unjamming, Bending, Granular material

Within the framework of the second order work theory, the onset of instabilities is explored numerically in loose granular materials through three dimensional DEM simulations. Stress controlled directional analysis are performed in Rendulic's plane and a particular attention is paid to transient evolutions at the microscale. Thanks to a micromechanical analysis, the onset and development of transient mechanical instabilities is explored. It is shown that these instabilities result from the unjamming and bending of a few force chains associated with a local burst of kinetic energy. This burst of kinetic energy propagates to the whole sample and provokes a generalized unjamming of force chains. As force chains buckle, a phase transition from a quasistatic to an inertial regime is observed. At the macroscopic scale, this results in a transient softening and a loss of controllability. After the collapse of existing force chains, the development of plastic strain is eventually stopped as new stable force chains are built.

INTRODUCTION

Granular materials are in the meantime both very simple and complex materials. Indeed, they can be described as simply as a set of grains interacting together through contact laws. The discrete nature of the assembly, however results in highly non-linear and non-regular constitutive behaviors in the sense that the direction of the strain rate depends on the loading direction and on the previous stress history. The overall response is thus inherently dynamic and relies widely on particles rearrangements. These dynamic effects are most of the time accounted for at the macroscale through velocity fluctuations in both quasi-static and inertial regime [START_REF] Roland | Towards a constitutive law for the unsteady contact stress in granular media[END_REF][START_REF] Johan | Quasistatic to inertial transition in granular materials and the role of fluctuations[END_REF][START_REF] Ryan | Friction in inertial granular flows: competition between dilation and grain-scale dissipation rates[END_REF], however when instabilities or failure are considered, this statistical quantity cannot account for sudden transitions from quasi-static to inertial regime [START_REF] François | Bifurcation in granular materials: An attempt for a unified framework[END_REF][START_REF] François | Inertia effects as a possible missing link between micro and macro second-order work in granular media[END_REF][START_REF] François | Failure in rate-independent granular materials as a bifurcation toward a dynamic regime[END_REF][START_REF] Richard | Diffuse instabilities with transition to localization in loose granular materials[END_REF]. In this case, only local instantaneous velocities are relevant to describe the rapid increase in kinetic energy. Based on the concept of the second-order work criterion introduced in the middle of the twentieth century by Hill [START_REF] Hill | A general theory of uniqueness and stability in elastic-plastic solids[END_REF], recent papers have demonstrated the ability of this criterion to anticipate the occurrence of material instabilities characterized by an outburst of kinetic energy [START_REF] François | Bifurcation in granular materials: An attempt for a unified framework[END_REF][START_REF] Ali | Diffuse failure in geomaterials: experiments, theory and modelling[END_REF][START_REF] François | Failure in rate-independent granular materials as a bifurcation toward a dynamic regime[END_REF][START_REF] Richard | Diffuse instabilities with transition to localization in loose granular materials[END_REF][START_REF] François | Second-order work criterion: from material point to boundary value problems[END_REF]. It has even been shown that the vanishing of the second order work is a necessary condition [START_REF] François | Bifurcation in granular materials: An attempt for a unified framework[END_REF][START_REF] Ali | Diffuse failure in geomaterials: experiments, theory and modelling[END_REF][START_REF] Richard | Diffuse instabilities with transition to localization in loose granular materials[END_REF] for the existence of any types of instabilities characterized for instance by localized failure [START_REF] Rudnicki | Conditions for the localization of deformation in pressure-sensitive dilatant materials[END_REF][START_REF] Nova | Controllability of the incremental response of soil specimens subjected to arbitrary loading programmes[END_REF] or by a classical Mohr-Coulomb plastic failure. This lead to the concept of bifurcation domain which consists in the set of states for which the second order work is negative. A bifurcation point corresponds to a potentially unstable state for which loading programs leading to the material failure exist. However not every loading program will necessarily trigger off the underlying instabilities.

Within the framework of small strain continuum mechanics, Hill's instability criterion applied to a material point states that for a given equilibrium ( , ) reached after a given loading history, if there exists at least one stress increment d , associated with a strain response d such that 2 = d ∶ d < 0, the material point is unstable. The physical meaning of Hill's criterion corresponds to a situation in which the deformation of the mechanical system can be pursued without any input of energy from the observer. From the writing of the energy balance, it can be shown that the variation of kinetic energy of a given mechanical system originally in an equilibrium state is a second order term d 2 equal to the difference between the external work resulting from the boundary conditions and the volume integral of the second order work [START_REF] François | A micro-mechanical investigation of bifurcation in granular materials[END_REF][START_REF] François | Bifurcation in granular materials: An attempt for a unified framework[END_REF][START_REF] François | Failure in rate-independent granular materials as a bifurcation toward a dynamic regime[END_REF]. Let Ω denotes a mechanical system of boundary Ω. Let be the force per surface unit applied on Ω and the displacement on Ω. Then the variation of kinetic energy from an equilibrium state to any new infinitely close state reads

d 2 = |Ω| ext 2 - 2 = ∫ Ω ⋅ d -∫ Ω 2 d (1)
where This fundamental equation written here in an updated Lagrangian formalism states that the external work is always larger than the integral of the second order work of the system. This has two important consequences:

-In case of a quasi-static evolution (d 2 = 0), the mean external second order work is equal to the mean internal second order work. As the external second order work is often easier to evaluate compared to its internal counterpart, numerical or experimental assessment of the second order work is carried out by computing the external second order work between two infinitely close equilibrium states.

-In the case of a nil external second order work, a negative second order work associated with an equilibrium position (at the first order d = 0) will automatically results in an increase in the kinetic energy of the system as d 2 > 0. As a result a transition from a quasi-static regime to an inertial one should be expected.

For granular materials, it has been proved that the second order work (at the scale of the continuum) can be computed from micro-quantities such as branch vectors and contact forces in the quasi-static case [START_REF] François | A micro-mechanical investigation of bifurcation in granular materials[END_REF][START_REF] François | From microscopic to macroscopic second-order work in granular assemblies[END_REF][START_REF] Nejib | Micromechanical analysis of second order work in granular media[END_REF][START_REF] François | Second-order work analysis for granular materials using a multiscale approach[END_REF]. However, to the best of our knowledge, the vanishing of the second order work has never been explicitly linked to physical processes occurring at the microscale. Thanks to the combine use of photoelastic grains and numerical simulations, the overall behavior of granular materials was shown to rely on a limited number of grains organized in mesostructures called force chains [START_REF] Drescher | Photoelastic verification of a mechanical model for the flow of a granular material[END_REF][START_REF] Liu | Force fluctuations in bead packs[END_REF][START_REF] Farhang | Bimodal character of stress transmission in granular packings[END_REF][START_REF] Bardenhagen | Numerical study of stress distribution in sheared granular material in two dimensions[END_REF][START_REF] Bernard | Micromechanics of granular materials[END_REF]. Following the definition as well as the algorithm given in Peters et al. [START_REF] Peters | Characterization of force chains in granular material[END_REF], recent results highlighted the importance of these mesoscale structures with respect to the overall mechanical properties of the granular materials [START_REF] Huaxiang | On a common critical state in localized and diffuse failure modes[END_REF][START_REF] Antoine | Scale separation between grain detachment and grain transport in granular media subjected to an internal flow[END_REF][START_REF] Antoinette | Force cycles and force chains[END_REF]. In particular, it has been shown, that the strength of a granular material results from its ability to build relatively long force chains and its ability to constantly rearrange existing ones to cope with any change in the boundary conditions. In addition, as shown by Zhu et al. [START_REF] Huaxiang | On a common critical state in localized and diffuse failure modes[END_REF][START_REF] Huaxiang | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF] and Tordesillas et al. [START_REF] Antoinette | Force cycles and force chains[END_REF] in 2D for instance, the particles not included in these chains forms cycles that have an important stabilizing role. The aim of this paper is to investigate numerically the elementary mechanisms taking place in loose granular materials when a loading program leading to the vanishing of the second order work has been detected. The results presented in this paper may also be useful to provide some insight on the standard numerical procedure used to assess the mechanical stability of a granular material thanks to the use of the second order criterion. This is done thanks to discrete element simulations on a representative elementary volume (REV) of a granular assembly. What is the micro-signature of an instability? What is the sequence of micromechanical mechanisms leading to the vanishing of the second order work at the continuum scale? These are the main questions that the manuscript should address. This paper is organized as follows. In section 2, a representative elementary volume of an idealized loose granular material is generated and submitted to a drained triaxial test. For different stress states, its mechanical stability is assessed in section 3 following a classical stress controlled directional analysis procedure with use of the second order work criterion. A particular care is paid to (i) the influence of the pre-stabilizing step of the procedure and (ii) the magnitude of the stress probe used in the directional analysis. In the last section of the paper, a micromechanical analysis of the onset and development of instabilities is performed for a particular loading increment corresponding to an observed vanishing of the second order work. The physical mechanisms at stakes are eventually identified.

SAMPLE PREPARATION AND TRIAXIAL TESTING

Numerical modeling

The micromechanical analysis performed throughout this paper considers non cohesive granular materials modeled as polydisperse assemblies of spheres. The interaction between two particles is modeled by the classical elasto-frictional contact law proposed by Cundal et al. [START_REF] Peter | A discrete numerical model for granular assemblies[END_REF] and is illustrated in Figure 1 . Two spherical particles are said to be in contact if they overlap. Based on the direction of the vector joining the particle centers, a normal force is defined as proportional to the overlapping distance Δ between the two particles where the contact normal stiffness is proportional to the material's Young modulus and to the harmonic average of the two particles radii 1 and 2 (see Fig. 1 ). In addition to the normal force, a tangential force is introduced. This tangential force is proportional to the relative tangential displacement Δ between the two particles where the horizontal contact stiffness is a fraction of its normal counterpart (see Fig. 1 ). Δ is defined in an incremental form as the relative tangential displacement of the contact point between a reference configuration (when the contact is first established) and the current configuration. The last parameter of the implemented contact law is the internal friction angle , which defines the largest accessible ratio ∕ according to the Mohr Coulomb theory (see Fig. 1 ). The input parameters used in this elasto-frictional contact law are reported in Table 1 . They are chosen equal to those of Hadda et al. [START_REF] Nejib | Micromechanical analysis of second order work in granular media[END_REF] and Wautier et al. [START_REF] Antoine | Scale separation between grain detachment and grain transport in granular media subjected to an internal flow[END_REF]. After computing all inter-particle contact forces, the induced particles displacements are integrated based on Newton's second law of motion using the DEM open source code YADE (28).

Sample preparation

A cubic assembly of spheres is generated randomly with a uniform radius distribution between min and max = 3.5 min . After generating a cloud of 10,000 non-overlapping spheres surrounded by six bounding planes defining a cube, the particles are inflated and allowed to rearrange according to the radius expansion technique. This process is stopped when the confining pressure applied on the bounding planes reaches 20 kPa and the unbalanced force ( unb ) of the system decreases below 10 -5 . By definition unb is equal to the mean summary force on particles (or equivalently to the mean mass time acceleration of each particle) divided by the mean contact force magnitude. As a result unb has no units and provide an intrinsic measure of how close to equilibrium a mechanical system is. During this process, the inter-particle friction angle is maintained to its value of 35 • in order to prepare a loose granular material with a void index of = 0.73. The resulting sample may be visualized in Figure 2 . Associated with the generated cubic grain assembly, a Cartesian coordinate system ( , , ) is defined such that the axis directions coincide with the edges of the cube.

Triaxial testing

From the obtained equilibrium state, a drained triaxial loading is then imposed to the sample in the form of a two step procedure. First, the confining pressure ( 0 ) is increased from 20 kPa to 100 kPa by allowing the bounding walls to move. Once a new equilibrium state is reached ( unb < 10 -5 ), a vertical compression strain rate ̇ = 0.01 s -1 is applied up to 6 % of deformation while keeping the same lateral confining pressure 0 . This strain rate is chosen similar to the one used in previous numerical studies [START_REF] Nejib | Micromechanical analysis of second order work in granular media[END_REF][START_REF] François | Failure in rate-independent granular materials as a bifurcation toward a dynamic regime[END_REF] and is supposed to be sufficiently small so that the loading can be considered as quasi-static. Throughout this paper the classical soil mechanics conventions are adopted with compressions counted positive. The homogeneous Cauchy stress tensor is defined at the REV scale from the forces applied on the bounding walls. The stress ratio is then introduced as the ratio between the deviatoric stress and the mean pressure

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ = = √ 3 2 dev ∶ dev dev = - = 1 3 Tr( ) . ( 2
)
where stands for the identity tensor and ":" stands for the double dot contraction product. Likewise, a macroscopic homogeneous strain tensor is defined from the bounding walls displacements. In order to be consistent with the second order work definition 2 = ∶ , compaction is counted positive. The volumetric strain is then simply defined as

= Tr( ) (3)
In Figure 3 , the stress ratio and volumetric strain responses are shown. A typical contractive behavior is observed as both and monotonously increase with the vertical compaction ( > 0). This behavior is typical of non compacted granular material (recent gravel deposits or stocks of industrial materials for instance). After a sharp increase, the stress ratio slowly increases on average from after less than 1 % of vertical strain (the critical state is not completely reached after 6 % of vertical strain). This second phase of the mechanical response is characterized by a "bumpy" curve resulting from important grain rearrangement processes occurring from time to time. As a result, as soon as these sudden drops in are observed the mechanical response of the sample becomes inherently dynamic and the quasi-static hypothesis becomes questionable over these short periods. In order to study the mechanical stability of the particular sample under different mechanical stress states, 9 samples are saved during the triaxial loading for ∈ {0.01, 0.20, 0.30, 0.35, 0.45, 0.50, 0.55, 0.60, 0.65}. The corresponding states are marked with diamonds in Figure 3 .

MACROSCOPIC ASSESSMENT OF BIFURCATION POINTS

In this section, the mechanical stability of the 9 mechanical states saved during the triaxial loading presented in the previous section is assessed thanks to a classical stress controlled directional analysis procedure with use of the second order work criterion [START_REF] Luc | Analysis of failure occurrence from direct simulations[END_REF][START_REF] François | From microscopic to macroscopic second-order work in granular assemblies[END_REF][START_REF] François | Bifurcation in granular materials: An attempt for a unified framework[END_REF]. This numerical procedure is composed of two steps which are reviewed in detail hereafter. First, a pre-stabilization is required to define a reference equilibrium state. Then, an incremental stress probe is imposed at the scale of the whole sample which is accompanied by a transient mechanical response until a second equilibrium state is reached. In addition to the directional analysis of the vanishing of the second order work criterion (section 3.2), a particular care is paid to the influence of the pre-stabilizing step (section 3.1), to the transient response observed between the two equilibrium states (section 3.3) and to the influence of the magnitude of the stress probe (section 3.4).

Pre-stabilization step

Even if the vertical compression rate is set at a reasonably low value, the mechanical response recorded during the triaxial loading is inherently dynamic as it relies on particle rearrangements and the integration of Newton's second law of motion. As a result, a pre-stabilization step is required prior to study the mechanical stability of saved samples thanks to the second order work criterion reviewed in the introduction. This pre-stabilization step is achieved by keeping the lateral pressure to the constant value 0 and by imposing a vertical stress corresponding to the stress ratio considered. The equilibrium state is assumed to be reached once unb < 10 -5 . Even if the stress state is kept constant, delayed deformations are observed during this pre-stabilization step. The final stressstrain states reached at the end of the pre-stabilization step are shown for each stress ratio in Figure 3 by blue circles. For ≤ 0.45 and = 0.55, negligible delayed deformations are observed while for ∈ {0.50, 0.60, 0.65} observed deformations are larger than 1 %. In addition, even when no important delayed deformation is observed, some localized bursts of kinetic energy may occur as illustrated in Figure 4 for = 0.45. These two observations prove the occurrence of local particle reorganizations which will affect the mechanical stability as the state variables obtained at the end are different from the initial ones as shown in particular for and in Figure 3 . This impact is even more important if the stress ratio is larger than = 0.45 as the quasistatic regime becomes questionable. The mechanical response of the sample is then governed by sudden rearrangements in the microstructure which yield fluctuations in the stress-strain response. These rearrangements are not instantaneous and continue to occur during the pre-stabilization step which result at the scale of the sample in delayed permanent deformations. A parallel can be drawn here with the concept of "delayed plasticity" introduced by Di Prisco ( 30) and also mentioned in Froiio et al. [START_REF] Francesco | Incremental response of a model granular material by stress probing with DEM simulations[END_REF].

Directional analysis

As introduced by ( 32) and widely used in the literature, directional analysis is a convenient framework to track the existence of stress increments leading to the vanishing of the second order work [START_REF] Jean-Pierre | Numerical simulations of the incremental responses of idealized granular materials[END_REF][START_REF] Luc | Analysis of failure occurrence from direct simulations[END_REF][START_REF] François | From microscopic to macroscopic second-order work in granular assemblies[END_REF][START_REF] François | Bifurcation in granular materials: An attempt for a unified framework[END_REF][START_REF] Barthélémy | Directional plastic flow and fabric dependencies in granular materials[END_REF]. In this study, we restrict our stability analysis to stress increments lying in the plane of axisymmetry such that d = d (Rendulic's plane). In this plane ( √ 2d , d ), a stress increment is fully described by its polar coordinates || || and such that

⎧ ⎪ ⎨ ⎪ ⎩ √ 2d = || || cos d = || || sin || || = √ d 2 + 2 d 2 . ( 4
)
In practice, finite stress increments of || || = 5 kPa are imposed in the form of a stress loading rate of 142 kPa.s -1 (corresponding to 10 000 numerical time increments) followed by a stabilization phase letting the system evolves toward a new equilibrium position as unb < 10 -5 . The size of the stress increment is similar to those classically used in the literature [START_REF] Calvetti | A numerical investigation of the incremental behavior of granular soils[END_REF][START_REF] Yuji | On the incremental nonlinearity observed in a numerical model for granular media[END_REF][START_REF] Luc | Analysis of failure occurrence from direct simulations[END_REF]. At the end, for each imposed stress increment , a strain increment is obtained. A normalized second order work is then defined at the material point scale (REV scale) as

norm 2 = ∶ || || || || (5)
Provided that the strain increment depends linearly with the applied stress increment for each direction (this will be investigated in section 3.4), the normalized second order work introduced in (5) depends only on and not || ||. For the considered stress increment || || = 5 kPa, the normalized second order work envelopes for the 9 pre-stabilized samples are shown in Figure 5 .

A circular representation is used as the second order work is positive outside the red dashed circle and negative inside.

In Figure 5 , the normalized second order work envelopes on the left correspond to the quasi-static regime ( ≤ 0.45) while the ones on the right corresponds to the dynamic regime of the triaxial test ( > 0.45). For small stress ratios, no vanishing of the second order work is observed. As the stress ratio increases, the second order work gradually decreases for incremental stress directions around 210 • and an instability cone ( values such that 2 ≤ 0) is eventually obtained for = 0.45. According to the second order work criterion the sample enters the bifurcation domain while approaching the onset of the dynamic regime during the triaxial loading as there exist incremental loading programs associated with negative 2 . Then, for higher stress ratio, the second order work suddenly becomes positive again for all the incremental stress directions (except for = 0.55 for which an instability cone is still visible). The non vanishing of 2 for ∈ {0.5, 0.6, 0.65} should be linked to the large delayed deformations observed for these stress state in Figure 3 . The grain rearrangements result in a stable configuration for which no stress increments lead to the vanishing of the second order work (at least for loading directions within Rendulic's plane).

Macroscopic analysis of the onset and development of instabilities

In the previous subsection, the second order work has been systematically computed between two equilibrium states. In this subsection, transient evolutions induced by macroscopic stress probes are analyzed at the scale of the REV. In Figure 6 , the incremental stress and strain time responses are presented for a stress ratio = 0.45 and for the two loading directions = 30.5 • (stable) and = 210.5 • (unstable). For the considered stable direction ( = 30.5 • ) the evolution of the incremental stress follows perfectly the prescribed loading program in the form of a ramp over the first 10 000 numerical time iterations ( < 0.035 s). By contrast the measured stress increments corresponding to the unstable direction = 210.5 • fails to follow prescribed ramp and a deviation is observed for ∈ [0.015, 0.04] s. Indeed, contrary to stable directions a transient loss of controllability (in the sense of Nova [START_REF] Nova | Controllability of the incremental response of soil specimens subjected to arbitrary loading programmes[END_REF]) is observed for unstable directions. This has to be linked to the development of large incremental strains up to a few percent whereas in the case of stable directions, the incremental strain is limited to approximately 10 -5 (see inset graph in Figure 6 ). The considered unstable direction ( = 210.5 • ) corresponds to a physical configuration in which the sample is slightly deconfined simultaneously in the vertical and horizontal directions (d < 0 and d = d < 0). For a stable material, this loading program should result in an increase in the volume of the sample which is only the case here for < 0.01 s (before large strains develop). A sudden collapse of the sample in the vertical direction is observed and a densification of the sample is obtained as the horizontal dilatancy does not counterbalance the vertical contraction in Figure 6 (d = ∑ d > 0). In Figure 7 , the transient evolution of the second order work is plotted for the two loading directions ∈ {30.5 • , 210.5 • } and for the same stress ratio = 0.45. It should be noted here that as the sample does not follow a quasi-static evolution, the transient second order work shown in Figure 7 is indeed the external second order work which is an upper bound for the internal second order work as d 2 ≥ 0 in equation (1) [START_REF] François | Second-order work criterion: from material point to boundary value problems[END_REF]. Initially positive, the second order work vanishes after = 0.014 s when the loss of controllability is observed. Then ext 2 decreases, goes through a minimum and eventually stabilizes after = 0.05 s around a negative value. This non monotonic evolution should be underlined as for some stress ratios and some stress loading directions the final increase of the second order 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 t (s) > 2 according to equation ( 1)) before getting back to its original position as soon as the softening regime is contained.

Influence of the stress increment on the onset of instabilities

As already pointed out in the literature concerning the numerical assessment of the mechanical stability of granular materials [START_REF] Francesco | Incremental response of a model granular material by stress probing with DEM simulations[END_REF][START_REF] Luc | Analysis of failure occurrence from direct simulations[END_REF], the stress increment || || used in the directional analysis cannot be infinitely small as assumed theoretically. This is usually assumed to be linked to the discrete nature of granular materials as the material point definition used at the basis of continuum mechanics is not rigorously fulfilled. The particular choice in || || is however arbitrary. In this study, the specific influence of || || on the normalized second order work envelopes shown in Figure 5 , has been explored for the stress ratio = 0.45. In Figure 9 , the normalized second order work envelopes are shown for six stress increments || || ∈ {0.5, 1, 2.5, 5, 7, 10} kPa. For each stress increment, the testing procedure remains the same as in section 3.2. In particular, the same stress loading rate || ̇ || = 142 kPa.s -1 is used. If no change is visible for stable directions, the magnitude of the applied stress increment has a noticeable impact on the width of the instability cone which can even disappear for small stress increments. Indeed, the onset of instabilities requires the applied perturbation to be large enough. Once all the underlying instabilities have been triggered off, the normalized second order work does not depend on || ̇ || anymore. In the present case, a minimal value between 2.5 kPa and 5 kPa is required. In order to explore this dependence on the stress increment at the microscale, it is interesting to define for each contact sliding index as

= tan - || || || || tan (6)
where is the inter-particle friction angle, and are respectively the tangential and normal contact forces as defined in section 2.1. This sliding index is normalized between 0 and 1, = 0 meaning that the Mohr-Coulomb criterion of the contact law (see Figure 1 ) is reached. In Figure 10 , probability density functions (pdf) associated with this sliding index are shown for the pre-stabilized sample corresponding to = 0.45. The solid line correspond to the reference pdf before the application of any stress increment while the dashed lines correspond to pdf reached after a stress increment || || = 0.5 kPa or || || = 5 kPa. Two loading directions are considered as = 30.5 • (stable direction) or = 210.5 • (unstable direction). In Figure 10 , the initial pdf presents a "S" shape with a local maximum for ≃ 0.75 and a maximum for = 0. Indeed, in the initial configuration, a significant number of contacts have reached or are close to the Morh-Coulomb limit. For both directions, the application of a stress increment || || = 0.5 kPa hardly change the contact distribution whereas a larger stress increment of || || = 5 kPa has a noticeable impact on the pdf as fewer contacts remain close to sliding. It should be noted that

FIGURE 11

Remaining fraction of contacts with a sliding index smaller than 0.1 with respect to the applied incremental stress magnitude for a stress ratio = 0.45 and an unstable direction = 210.5 • (solid line) and a stable direction = 30.5 • (dashed line).

the shape of the pdf remains mostly unchanged for > 0.1. As a result, a finite stress increment is required to have an impact on the microstructure geometry which is required to observe the vanishing of the second order work [START_REF] Luc | Analysis of failure occurrence from direct simulations[END_REF][START_REF] Francesco | Incremental response of a model granular material by stress probing with DEM simulations[END_REF]. The comparison between the two loading directions for || || = 5 kPa shows that for the unstable direction = 210.5 • , the probability density falls down to 0 for = 0 which is not the case for = 30.5 • . Indeed, an incremental load in a stable direction will provoke a partial reorganisation of the microstructure through gain and loss of contacts, whereas an incremental load in an unstable direction will result in a complete microstructure reorganization. In both cases many contacts close to sliding are either lost or released and for = 210.5 • no contacts remain close to sliding anymore. In order to highlight the specific influence of the incremental stress || || on contact reorganizations, the fraction of contacts remaining close to sliding ( < 0.1) was investigated for ∈ {30. || remaining fraction of contacts close to sliding drops to a low value around 30 % and remains stable. These observations are consistent with the comments made previously and can be interpreted through a micromechanical approach of stability based on jamming analysis [START_REF] Bagi | On the concept of jammed configurations from a structural mechanics perspective[END_REF]. Indeed, an incremental loading along a stable direction will result in incremental modifications of the microstructure whereas an incremental loading along a unstable direction will result in a generalized unjamming/re-jamming process, provided that || || > || * ||. As a result, for unstable directions, there exists a threshold value responsible for the triggering of the observed instability. As soon as the threshold is reached, the size of the stress increment does not have any influence on the number of contacts close to sliding anymore. It should also be noted that this value is found to depend on the loading rate and the direction of loading (not shown here). The existence of a threshold value is consistent with DEM framework in which particle overlapping is allowed and controlled through a normal stiffness . As a result, contact points are not limited to single points and the loss of contact requires a finite perturbation. A complementary study to the present work (not shown here) has been carried out to highlight the dependence of || * || to . An increase in reduces grain overlap and a smaller threshold is observed.

MICROMECHANICAL ANALYSIS OF THE ONSET AND DEVELOPMENT OF INSTABILITIES

In the previous section, the mechanical stability of granular samples has been assessed at the scale of the REV and the onset of instabilities has been explained as the ability of an incremental load to trigger off microstructure reorganizations. If the mechanical state of the considered sample is in the bifurcation domain (or very close to it) and loaded along an unstable direction, the incremental loading induces generalized microstructure reorganizations which result in a macroscopic transient softening responsible for a loss of controllability. In the end, a new equilibrium is reached which is characterized by a contact population relatively far from sliding. The purpose of this section is to provide a microscale investigation of the physical processes leading to the vanishing of the second order work during the transient loss of controllability phase observed macroscopically in Figure 6 . In all this section, the particular stress state = 0.45 and the unstable loading direction = 210.5 • are considered. A particular attention is paid to the chronology of events leading to the softening of the granular assembly.

Spatial distribution of kinetic energy

As previously highlighted in the literature [START_REF] Claudio | Experimental analysis and theoretical interpretation of triaxial load controlled loose sand specimen collapses[END_REF][START_REF] Luc | Analysis of failure occurrence from direct simulations[END_REF][START_REF] François | Failure in rate-independent granular materials as a bifurcation toward a dynamic regime[END_REF][START_REF] François | Bifurcation in granular materials: An attempt for a unified framework[END_REF], the main ingredient enabling microstructure reorganizations in granular materials is the particles' kinetic energy. As a result, it is of particular interest to track the time evolution of the kinetic energy of the individual particles while applying a stress increment. In Figure 12 , snapshots of the considered sample are shown for different time steps for two sufficiently large stress increments || || ∈ {1, 5} kPa (see Figure 11 ). The particles are colored according to their kinetic energy, and the most energetic ones are highlighted. An arbitrary threshold of * = 10 -8 J is chosen corresponding to the most energetic particles in the initial state.

In both cases, a localized burst of kinetic energy appears at the same spot and approximately at the same time (once || || has reached its targeted value in the 1 kPa case and during the transient increase of || || in the 5 kPa case). Then the local burst of kinetic energy propagates to the whole sample. In the case where || || is below the threshold value || * || identified in section 3.4, no burst of kinetic energy is visible (not shown here). Indeed, the observed threshold value in direction = 210.5 • corresponds to the minimal perturbation required to trigger off the burst of kinetic energy shown in Figure 12 . Once initiated, the burst of kinetic energy can propagate to the whole sample. The sample is unstable as a whole and a suitable perturbation will trigger off the underlying instability leading to a generalized microstructure rearrangement. This contrasts with local bursts of kinetic energy observed during the pre-stabilization step in Figure 4 . The observed vanishing of the second order work is thus a material property and not a structural one as the microstructure modifications do not stay localized in some regions of the sample with a length scale similar to the one of the whole sample. This is the reason why we can still consider this instability as a material instability. The sample is unstable as a whole and a small perturbation will trigger off the underlying instability. A structural instability would have stayed confined in a restricted zone of our sample with a length scale similar to the one of the whole sample. In order to describe the onset and the propagation of the burst of kinetic energy, the time evolution of the mean kinetic energy per particle is plot in Figure 13 for the whole sample and for the small control volume visible in Figure 12 . In Figure 13 , an early increase in the mean kinetic energy per particle occurs simultaneously in the control volume for ||d || = 1 kPa and ||d || = 5 kPa for ≃ 0.005 s. This corresponds to the initiation of the burst of the kinetic energy shown in Figure 12 . After ≃ 0.015 s the kinetic energy per particle follows the same evolution in the control volume and for the whole sample. This marks the end of the propagation of the burst of the kinetic energy. For ||d || = 1 kPa, the transient incremental loading is stopped before the end of the kinetic burst propagation (blue dotted line). As a result, a small plateau is observed in the general increase in the kinetic energy corresponding to the time needed for the instabilities to affect the whole sample. Such a plateau is not observed for ||d || = 5 kPa as kinetic energy is continuously provided through evolving boundary conditions (the end of the transient loading is reached even after the mean kinetic energy per particle reaches its maximum).

Force chains definition

Now commonly used in the literature, the force chain concept provides a relevant mesoscopic scale to account for the macroscopic mechanical behavior of granular materials [START_REF] Huaxiang | On a common critical state in localized and diffuse failure modes[END_REF][START_REF] Zhang Lingran | The role of force chains in granular materials: from statics to dynamics[END_REF][START_REF] Antoine | Scale separation between grain detachment and grain transport in granular media subjected to an internal flow[END_REF][START_REF] Antoinette | Force cycles and force chains[END_REF]. In particular, it has been shown, that the strength of a granular material results from its ability to build relatively long force chains and to constantly rearrange the existing force chains to cope with any change in the boundary conditions [START_REF] Antoinette | Structural stability and jamming of self-organized cluster conformations in dense granular materials[END_REF][START_REF] Huaxiang | On a common critical state in localized and diffuse failure modes[END_REF][START_REF] Antoine | Scale separation between grain detachment and grain transport in granular media subjected to an internal flow[END_REF]. Details on the algorithm used to identify the chained particles of a sample can be found in Peters et al. [START_REF] Peters | Characterization of force chains in granular material[END_REF]. The definition of a force chain used throughout this manuscript is briefly reviewed:

-The particles belonging to a force chain have a higher principal stress than the mean particle principal stress. -The principal stress direction of chained particles is aligned with the geometrical contact direction (less than 45 • deviation).

-A force chain contains at least three contacting particles.

Based on this definition, the particles responsible for the stress transmission can be identified in the initial state and tracked while the incremental stress is applied. By comparing the current set of chained particles to the initial one, force chain renewal can be quantified as born current chained particles were not identified as such initially and died initially chained particles are no longer identified as such. A visualization of the force chains identified in the initial and final configurations is shown in Figure 14 as well as the time evolution of the total number of chained particles ( tot ) and the time evolution of the new and former ones.

In this figure the microstructure evolution induced by the applied incremental stress is visible. Between the initial and final states roughly 50 % of the chained particles are renewed while the total number of chained particles remains quasi-constant (a slight decrease from 2821 to 2661 chained particles is observed). In addition, in the final state the vertical anisotropy of the force chains seems visually to be less pronounced than initially.

Chained particles lifespan and life expectancy

In order to detect birth and death of force chains, the notions of lifespan and life expectancy are introduced similarly as in [START_REF] Huaxiang | On a common critical state in localized and diffuse failure modes[END_REF]. This analysis is not carried out directly on force chains but on groups of three chained particles. For a given time and a given group of 3 chained particles , the lifespan ( , ) is defined as

( , ) = -birth - , ( , ) ∈ [0, 1] (7)
where = 0 is the time at which the incremental stress is applied and birth is the time at which the three particles composing were first identified as chained particles and from which they remain continuously identified as such until the current time . ( , ) = 1 means that the group of three chained particles have been existing since the initial configuration before the incremental stress is applied. On the contrary, ( , ) = 0 means that has just appeared. Similarly, the life expectancy ( , ) of computed at time is defined as

( , ) = death - - , ( , ) ∈ [0, 1] (8)
where = 0.087 s is the time at which a new equilibrium ( unb < 10 -5 ) is reached after the application of the incremental stress and death is the time which the three particles composing were last identified as particles starting from time . ( , ) = 1 means that the group of three chained particles exists until the final stabilization of the sample. On the contrary, ( , ) = 0 means that is about to disappear. In Figure 15 , the probability density functions and are shown for different time ∈ [ , ] in order to identify whether the current force chains are rather young or old and whether they have a long life ahead or are about to disappear. In the beginning of the incremental loading, the force chains are stable as all groups of three chained particles are old with a non zero life expectancy. This statement remains true when the burst of kinetic energy identified in the previous subsection initiates and propagates to the whole sample ( ∈ [0.005, 0.013] s). From the moment that the kinetic energy reaches its maximum level ( = 0.0388 s), groups of three particles start collapsing replaced by young force chains with a short life expectancy. Once the kinetic energy of the sample start decreasing (from = 0.0634 s) some time after the targeted stress increment reaches its final value ( = 0.035 s), force chains become stable again as an aging population of force chains with a long life expectancy is encountered. The main conclusions to be retained from this analysis are that:

-the initiation of the burst of kinetic energy is not due to the destruction of force chains; -the loss of controllability appears prior to the generalized collapse of force chains; -the force chains reorganization appears simultaneously with the maximal level of kinetic energy; -the final stabilization of the sample is reached as soon as long lasting force chains appear.

Localized force chains bending

As the loss of controllability occurs prior to the destruction of force chains, their collapse should not be regarded as the triggering microscale mechanism responsible for the onset of the burst of kinetic energy [START_REF] Antoinette | Force cycles and force chains[END_REF]. Force chains are elongated column like structures loaded in compression. As a result, their effective failure is very likely to be related to the onset of bending, before destruction of any force chains. Given a group of three chained particles for which two contact directions form an angle ∈ [0, ], the bending rate ̇ is simply the time derivative of . A strictly positive value ̇ > 0 characterize the bending of . In Figure 16 , a zoom of the time evolution of the mean kinetic energy and the mean bending rate are shown for the whole sample and a small control volume around the location of the burst of kinetic energy (see Figure 12 ). On the first graph, the evolution of the difference between the actual incremental stress components and their prescribed values is shown. In this figure, both the loss of controllability, the onset of the burst of kinetic energy and the increase in the bending rate in the control volume shown in Figure 12 occur simultaneously for ≃ 0.008 s, approximately 0.007 s prior to a general increase in mean kinetic energy per particle and mean bending rate for the whole sample (black solid curves in Figure 16 ). Indeed, the localized bending of a few force chains seems to be sufficient to generate a loss of controllability at the scale of the REV. As recently shown for 2D granular assemblies, force chain loss of stability results from the opening of contacting grain cycles [START_REF] Antoinette | Force cycles and force chains[END_REF][START_REF] Antoinette | Structural stability and jamming of self-organized cluster conformations in dense granular materials[END_REF][START_REF] Huaxiang | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]. In order to investigate this feature in 3D, the time evolutions of the number of contacts between two chained particles ( ), two non-chained particles ( ) and a non-chained and a chained particle ( ) are represented in Figure 17 . As for Figure 16 , the time evolutions are given for the whole sample and for the same control volume around the location of the burst of kinetic energy (see Figure 12 ). The onset and the propagation of the burst of kinetic energy are indicated by two vertical solid lines. In Figure 17 , the number of non-chained/chained contacts ( ) decreases in the control volume even prior to the onset of the burst of kinetic energy while and remain more or less constant up to the end of the propagation of the burst of kinetic energy. This observation suggests that the localized bending highlighted in Figure 16 results from an early unjamming of force chains within the used control volume (and consequently a local dilatancy). The evolution of , and for the whole sample confirms this mechanism with a generalized drop in observed during the propagation of the burst of kinetic energy. This failure mechanism is consistent with the 2D opening of the three particles loops around force chains mentioned in Tordesillas et al. [START_REF] Antoinette | Structural stability and jamming of self-organized cluster conformations in dense granular materials[END_REF] or Zhu et al. [START_REF] Huaxiang | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]. In the end, following the generalized collapse of force chains, the initially loose sample gets denser and the number of contacts between non-chained particles increases. The decrease in follows qualitatively the trend mentioned in section 4.2 for the number of chained particles. However, as the number of chained particles is reduced by roughly 160, drops by approximately 260. This might result from changes in the topology of force chains (less branching or shorter chains for instance). Complementary to Figure 17 , the evolution of the mean coordination number of the particles belonging to force chains ( Fch ) is shown in Figure 18 for ||d || ∈ {1, 5} kPa. In Figure 18 , the evolution of the mean kinetic energy for the whole sample is shown in dashed lines while the onset and propagation of the burst of kinetic energy are still indicated by two vertical solid lines. As soon as the incremental load starts, Fch drops which is a signature for local dilatancy around force chains. This decrease in Fch appears before any noticeable increase in kinetic energy and its decreasing rate accelerates as the local burst of kinetic energy appears. These observations are consistent with the failure mechanism identified previously. In addition, it should be noted that drops in Fch occur simultaneously for || || = 1 kPa and || || = 5 kPa (both above the threshold value identified in section 3.4) so as the onset of the burst of kinetic energy. In the case || || = 1 kPa, the drop is observed after the end of the transient incremental loading while for || || = 5 kPa, the drop occurs during the transient increase of || || from 0 to 5 kPa (see Figure 13 ). As a result, the failure mechanism seems to be triggered off as soon as the incremental stress overcomes the threshold value || * || = 0.9 kPa (see section 3.4). Then, the failure is not instantaneous and develops in an intrinsic timespan which is independent of the size of the stress increment. The concept of delayed instability observed experimentally in soils by Di Prisco et al. [START_REF] Claudio | Experimental analysis and theoretical interpretation of triaxial load controlled loose sand specimen collapses[END_REF] is recovered here at the microscale on a shorter timescale.

FIGURE 16

Absolute error between the actual stress components and their targeted values during the incremental stress loading (top), mean kinetic energy per particle for the whole sample and close to initiation of the kinetic burst (bottom left), mean bending rate of the groups of three chained particles for the whole sample and around the kinetic burst (bottom right). The sample considered is characterized by a stress ratio = 0.45 and is loaded along an unstable direction = 210.5 • with a stress increment ||d || = 5 kPa.

CONCLUSION AND OUTLOOK

Thanks to the combine use of the second order work theory at the macroscale and the force chain concept at the microscale, physical mechanisms responsible for the onset and development of instabilities in loose granular material at the scale of the continuum have been identified at the microscale. It was shown that an incremental stress loading leading to the vanishing of the second order work provokes the unjamming of force chains. This results in their bending accompanied by an increase in the kinetic energy. At the macroscale, force chains bending is visible through a softening behavior and a loss of controllability as the actual incremental stress rotates toward the direction of the instability cone. Initially localized, this phenomenon propagates to the whole sample in the form of a burst of kinetic energy. Once the kinetic energy of the whole sample reaches a sufficient value, existing force chains collapse and important microstructure reorganizations are observed resulting in a macroscopic densification. This densification is eventually stopped as soon as new stable force chains are built. Macroscopic softening ends and the incremental stress rotates back to its prescribed direction. In addition to this physical interpretation of the vanishing of the second order work, a particular attention has been paid to the pre-stabilization step required to define a reference equilibrium state as well as to the magnitude of the stress increments required ). The evolution is given for the whole sample (left) and for the control volume shown in Figure 12 (right). The onset and the propagation of the burst of kinetic energy is delimited by two vertical lines. The sample considered is characterized by a stress ratio = 0.45 and is loaded along an unstable direction = 210. to trigger off instabilities. Even for quasi-static evolutions, it was shown that localized micro-inertia crisis occur during this pre-stabilization step. As a consequence, threshold values exist for stress increments below which the sign of the normalized second order work depends on the magnitude of stress increments. These observations are hardly ever reported in the literature and could provide clues to carry out systematic numerical assessment of the bifurcation domain of granular materials. As outlook, the micromechanical approach presented in this paper could provide valuable tools to further explore the influence of some specific microstructure modifications onto the mechanical stability of granular materials. For instance, the specific role of rattlers (i.e. particles with no contacts in the absence of gravity) can be assessed. Even if they are useless to withstand the initial stress state, they are mobilized as soon as force chains collapse. As a result they are expected to play a stabilizing role. In addition, more subtle microstructure evolutions may be explored as those induced by suffusion, a specific type of internal erosion processes. Suffusion occurs in saturated granular materials submitted to an internal flow. It has a noticeable impact on free particles that are either flushed out or clogged. The resulting impact on the mechanical stability remains today an open question.

1 |Ω|∫ Ω 2 d

 12 to as the mean external second order work and 2 = is the mean (internal) second order work.

FIGURE 1

 1 FIGURE 1 Elasto-frictional contact law used in DEM simulations.

FIGURE 2

 2 FIGURE 2 Sample visualization in the initial state.

FIGURE 3

 3 FIGURE 3Drained triaxial test response. The diamond points correspond to the mechanical states considered for the stability analysis. The circles corresponds to the stress-strain state reached after the pre-stabilization step (see subsection3.1)

FIGURE 4

 4 FIGURE 4 Two typical visualization of observed local bursts of kinetic energy during the pre-stabilizing step ( = 0.45). Particles are colored according to their kinetic energy level. Particles with > 10 -8 J are highlighted.

FIGURE 5

 5 FIGURE 5 Circular second order work for different values in quasi-static regime (left) and in dynamic regime (right). 2 = 0 is indicated by a red dashed circle.

5 FIGURE 6

 56 FIGURE 6 Macroscopic strain (left) and stress (right) response for the two loading directions = 30.5 • (dashed lines) and = 210.5 • (solid lines). The vertical stress (d ) and strain (d ) components are plotted in blue while the horizontal ones are plotted in red (d and d ) and green (d and d ).

FIGURE 7 2 = 2 1 FIGURE 8 2 = 0

 7221820 FIGURE 7External normalized second order works for the two loading directions = 30.5 • (dashed line) and = 210.5 • (solid line). The vanishing of 2 is indicated by the horizontal red line.

FIGURE 9 ddFIGURE 10

 910 FIGURE 9Circular second order work for = 0.45 and different stress increment magnitudes. 2 = 0 is indicated by a red dashed circle.

5 •

 5 ; 210.5 • } and = 0.45 for 14 values of || || in Figure 11 . In Figure 11 , stable and unstable directions show completely different results. For = 30.5 • , the number of contacts remaining close to sliding smoothly decrease with || || whereas a threshold value || * || = 0.9 kPa is observed for = 210.5 • . For || || < || * ||, the final number of contacts ( < 0.1) close to sliding increases monotonously by comparison with the initial one ( < 0.1). Then, for || || > || *

FIGURE 12 FIGURE 13

 1213 FIGURE 12 Onset and propagation of the burst of kinetic energy leading to the loss of controllability for || || = 1 kPa (a) and || || = 5 kPa (b). The sample considered is characterized by a stress ratio = 0.45 and is loaded in an unstable direction = 210.5 • . Particles with > 10 -8 J are highlighted. A control volume around the location of the kinetic burst is shown in black.

FIGURE 14

 14 FIGURE[START_REF] François | From microscopic to macroscopic second-order work in granular assemblies[END_REF] Visualization of force chains identified in the initial (a) and final (b) configurations for a stress ratio = 0.45 in an unstable direction = 210.5 • for || || = 5 kPa. Force chains are composed of contacting particles of same color. (c) Time evolution of the number of chained particles (black), the number of new (green) and former (red) chained particles compared with the initial state.

  .0775 s, close to the second equilibrium state

FIGURE 15

 15 FIGURE[START_REF] Nejib | Micromechanical analysis of second order work in granular media[END_REF] Probability density functions associated to the lifespan (blue) and life expectancy (red) of the groups of three chained particles at different moments of the incremental loading. The sample considered is characterized by a stress ratio = 0.45 and is loaded in an unstable direction = 210.5 • with a stress increment ||d || = 5 kPa. A high value of around 1 corresponds to the presence of old force chains, a high value of around 1 corresponds to the presence long lasting force chains.

FIGURE 18

 18 FIGURE 17 Time evolution of the number of contacts between two chained particles (), two non-chained particles ( ) and a non-chained and a chained particle (). The evolution is given for the whole sample (left) and for the control volume shown in Figure12(right). The onset and the propagation of the burst of kinetic energy is delimited by two vertical lines. The sample considered is characterized by a stress ratio = 0.45 and is loaded along an unstable direction = 210.5 • with ||d || = 5 kPa.

TABLE 1

 1 Mechanical parameters used in the elasto-frictional contact law implemented in YADE.

	Parameters	Value
	Density	3,000 kg.m -3
	Young Modulus ( )	356 MPa
	Stiffness ratio ( )	0.42
	Inter-particle friction angle ( )

  Time evolution of the number of contacts between two chained particles ( ), two non-chained particles ( ) and a non-chained and a chained particle (

		5400			Whole sample		125	Control volume around kinetic burst
		5000 5200						120
	N nn , N nc , N cc	4200 4400 4600 4800				N cc N nc N nn	N nn , N nc , N cc	105 110 115
		3800 4000						100
		0.00 3600	0.02	0.04 t (s)	0.06	0.08		0.000 95	0.004	0.008 t (s)	0.012	0.016

FIGURE 17