Keywords: Granular materials, Critical state, Proportional loading path, Dynamical system, Microstructure, Multiscale modeling, DEM, Micromechanics

my friends Lingling Xuan, Xinming Xing, Dou Quan, Yanning Hu, Bo Ren and Qiangqiang Sun for enjoyable and unforgettable time with them after work during the three years.

In the end, my sincere gratitude to my family for their unwavering love and support all the time and deep thanks to Jiangfeng with whom I can share all my emotions, thoughts and dreams.

Abstract

Multiscale modeling of granular materials has received an increasing attention in recent years based on notable progress in instrumentation technologies and computing power. It will probably play a prominent role in the future of constitutive modeling as an alternative method to phenomenological models in large-scale engineering problems, especially those concerned by heterogeneity and strain localization, e.g., landslides. In multiscale modeling, the main challenge is to simulate continuum-scale behavior using information obtained at finer scales. As for the continuum-scale behavior in granular materials, critical state is one of the most important features that manifests in drained triaxial tests by the capacity of a geomaterial to accommodate shear deformation under constant stresses and without any volume change.

A fresh perspective is presented in this thesis for the standard concept of critical state (CS) in granular materials by suggesting that CS can be defined using a single proportional strain test. A comparison between proportional strain tests and biaxial tests simulated with discrete element methods (DEM) shows that the CS line characterized by the stress state (effective pressure and deviatoric stress), void ratio and fabric indexes can act as an attractor. The mechanical responses and fabric metrics evolve along dilatant proportional strain loading paths according to similar values after the strain level has become large enough to wipe out the initial material memory in the homogeneous domains considered in this analysis, i.e., the shear band domain in dense samples subjected to strain localization and the whole sample domain in loose samples. This suggests that the contact network of a granular material subjected to proportional strain loading evolves while preserving its ability to withstand shearing without volume change. Therefore, the CS concept can be found in a wide class of loading paths which shows that CS acts as a general attractor irrespective of the loading path considered.

Once critical state is reached, a stationary regime takes place in a biaxial test. The macroscopic state remains constant while the microstructure is constantly rearranging. In this thesis, we examine Critical State in granular materials from a dynamical systems theory (DST) perspective. The aim is to enrich the classical view of the critical state regime in granular materials, and more broadly, improve fundamental understanding of the underlying mesoscale mechanisms responsible for macroscopic stationary states in complex systems. The lifespan and life expectancy of specific cluster conformations comprising particles in force chains and grain loops, are tracked and quantified. Results suggest that these conformational clusters reorganize at similar rates in the critical state regime, depending on strain magnitudes and confining pressures. We quantified this rate of reorganization i ii and found that the material memory rapidly fades, with an entirely new generation of force chains and grain loops replacing the old ones within a few percent strain.

Based on these extensive investigations on the critical state, mesostructural transformations are found which bring vanishing and generating mesostructures along the loading path as a central ingredient that should be included in multiscale constitutive modeling. A deactivation/reactivation procedure acting on the local mesoscale is proposed in a specific micromechanical-based constitutive model (the H-model). Preliminary results show that this is a promising direction to improve the model to predict relatively large strains and to make stationary regimes emerge naturally without drawing from any empirical law at the macroscopic scale.

Résumé

La modélisation multi-échelles des matériaux granulaires a fait l'objet d'une attention croissante ces dernières années, grâce aux progrès notables de des techniques d'instrumentation et à l'augmentation des puissances de calcul. Il est probable qu'elle jouera un rôle essentiel à l'avenir en tant que méthode alternative aux modèles phénoménologiques dans les problèmes d'ingénierie à grandes échelles, en particulier ceux qui concernent l'hétérogénéité et la localisation des déformations, par exemple, dans la modélisation des glissements de terrain. Dans la modélisation multi-échelle, le principal défi est de simuler le comportement à l'échelle du continuum en utilisant au maximum les informations obtenues aux échelles plus fines dans le système, plutôt que de recourir à des modèles constitutifs empiriques. En ce qui concerne le comportement à l'échelle du continuum dans les matériaux granulaires, l'état critique est l'une des caractéristiques les plus importants, que l'on définit généralement dans des essais triaxiaux drainés comme la capacité pour un géomatériau à encaisser de la déformation de cisaillement à état de contrainte constant et sans changement de volume.

Une approche originale est présentée dans cette thèse concernant concept classique d'état critique dans les matériaux granulaires en suggérant que l'EC peut être défini par l'utilisation d'un seul essai en déformation proportionnelle. Une comparaison entre des essais en déformation proportionnelle et des essais biaxiaux simulés par la méthode aux éléments discrets (DEM) montre que la ligne d'état critique caractérisée par l'état de contrainte (pression effective et contrainte déviatorique), l'indice des vides et le déviateur du tenseur de texture peut agit comme un attracteur. Toutes les réponses mécaniques et les indicateurs de microstructure évoluent lors des chemins de chargement dilatés à déformation proportionnelle selon des valeurs similaires après que le niveau de déformation soit suffisamment important pour effacer la mémoire initiale du matériau dans les domaines homogènes considérés dans cette analyse, c'est-à-dire la zone de la bande de cisaillement dans les échantillons denses soumis au phénomène de localisation de la déformation et l'échantillon tout entier dans les échantillons lâches. Cela suggère que la microstructure d'un matériau granulaire soumis à un chargement de déformation proportionnelle évolue avec une capacité constante à résister au cisaillement sans changement de volume. Par conséquent, le concept d'état critique peut être retrouvé dans une large classe de chemins de chargement, puisqu'il agit comme un attracteur, quel que soit le chemin de chargement considéré.

Une fois l'état critique atteint, un régime stationnaire peut se mettre en place dans un essai biaxial. L'état macroscopique reste inchangé alors que la microstructure se réarrange en permanence. Dans ce travail de thèse, nous examinons ces régimes stationnaires iii iv dans les matériaux granulaires du point de vue de la théorie des systèmes dynamiques (DST). L'objectif est d'enrichir la vision classique de l'état critique dans les matériaux granulaires et, plus largement, d'améliorer la compréhension fondamentale des mécanismes mésoscopiques sous-jacents responsables des états stationnaires macroscopiques dans les systèmes complexes. La durée de vie et l'espérance de vie de conformations spécifiques de clusters, comprenant quelques particules comme les chaînes de force et les cycles de grains, sont suivies et quantifiées. Les résultats suggèrent que ces structures se réorganisent à des vitesses similaires à l'état critique, proportionnellement aux amplitudes de déformation et aux pressions de confinement considérées. Nous avons quantifié ces taux de réorganisation et constaté que la mémoire du matériau s'estompe rapidement, une génération entièrement de nouvelles chaînes de force et cycles de grains remplaçant l'ancienne en quelques pourcentages de déformation.

Sur la base de ces études micromécaniques détaillées sur l'état critique, on constate que les transformations mésostructurales correspondant à la disparition et la génération de mésostructures au cours du chargement est un ingrédient central qui devrait être inclus dans les modèles multi-échelles du comportement constitutif des matériaux granulaires. Une procédure de désactivation/réactivation de mésostructures est proposée dans un modèle micromécanique particulier : le modèle H. Les résultats préliminaires montrent que cette procédure constitue une direction prometteuse pour améliorer le modèle et faire émerger naturellement un régime stationnaire sans s'appuyer sur une loi empirique à l'échelle macroscopique.

Mots clés : Matériaux granulaires, État critique, Chemin de chargement proportionnel, Systèmes dynamiques, Microstructure, Modélisation multi-échelle, DEM, Micromécanique.

List of Figures

2.1 Typical deviatoric stress and volumetric strain responses in a dense (left) and a loose (right) granular assembly under a biaxial loading path. . . . .

2.2 images of a system of photoelastic disks for an isotropically compressed state (left) and a sheared state (right) [START_REF] Luding | Contact force measurements and stress-induced anisotropy in granular materials[END_REF] . . . . .

2.3

Avatar specimen under isotropic stress (left) and at the end of triaxial loading (right) [START_REF] Kawamoto | All you need is shape: predicting shear banding in sand with LS-DEM[END_REF] 3.5 Principal stresses in the whole sample (σ 1 and σ 2 ) and within the shear band (σ * 1 and σ * 2 ), and relative orientation of principal stress (θ * σ 1 -θ σ 1 ). The 100kPa-dense sample under biaxial loading path is adopted. . . . . . .

3.6

Stress paths of biaxial tests and proportional strain tests in p -q plane.

The start and end points are marked by triangles and squares, respectively. The critical stress ratio curve and the maximum stress ratio curve are drawn according to the critical states and the maximum value from biaxial tests, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.7

Stress paths in p -q plane along six proportional strain tests and one biaxial test. Dilatant (λ = -1.2), undrained (λ = -1) and contracting (λ = -0.5) proportional strain paths are conducted with the samples labelled 100kPadense and 100kPa-loose. A biaxial loading path is performed in the sample 100kPa-dense. The square domain is zoomed up on the right hand side. . .

3.16

Stress-strain responses along combined proportional strain and biaxial loading paths. The switch states and ultimate states are marked by dots and squares, respectively. The corresponding stress paths are shown in Fig. 3.15

3.17 The relation between imposed dilatancy from proportional strain test and unconstrained dilatancy if switched to a biaxial condition according to results in Fig. 3.15 and Fig. 3.16. Note that soil mechanics convention is adopted with positive compression. . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.18

The relation between imposed dilatancy/contractancy in proportional strain tests and unconstrained dilatancy/contractancy expected in biaxial tests in the plane of axial strain vs. volumetric strain (left) and the possible stress paths in the plane of deviatoric stress vs. mean stress (right). The major solid lines represent the imposed dilatancy/contractancy, namely the proportional strain loading paths; the fins denote the incremental unconstrained dilatancy/contractancy if the loading is switched to biaxial conditions. . . 4.1 Evolution of deviatoric stress q and volumetric strain ε v with respect to the axial strain ε 2 along biaxial tests. Note that soil mechanics convention is adopted with positive compression and positive contraction. . . . . . . . . 4.5 Lifespan and life expectancy PDFs (probability density functions) of chained particles in two samples: (a) 100kPa-dense sample where only the spatial domain within the shear band is considered and (b) 100kPa-loose sample with the whole domain considered. For each sample, two reference conformations belonging to the stationary regime are adopted, labeled as M and N. The deviatoric stress curve (dot-dashed line) and the volumetric strain (dashed line) curves are recalled. . . . . . . . . . . . . . . . . . . . . . . . 66 4.6 The averaged PDFs of chained particles with respect to the incremental strain evolution (∆ε 2 ), over six microstructure conformations at different strain levels belonging to the critical regime (ε 2 = 0.0610, 0.0708, 0.0806, 0.0904, 0.1002 and 0.1100 for the 100kPa-dense sample and ε 2 = 0.0605, 0.0708, 0.0802, 0.0906, 0.1005 and 0.1106 for the 100kPa-loose sample). . . 67 4.7 The strain magnitudes within the shear band domain (100kPa-dense*), whole dense sample (100kPa-dense) and the whole loose sample (100kPaloose). The strain magnitude ratio between 100kPa-dense* and 100kPadense is around 2.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.8 PDFs for chained particles from 100kPa-loose sample and from 100kPa-dense sample stretched by the strain magnitude ratio 2.8. . . . . . . . . . . . . . 68 4.9 Lifespan and life expectancy PDFs (probability density functions) of grain loops in two samples with respect tho the axial strain ε 2 : (a) 100kPa-dense sample where only the spatial domain within the shear band is considered and (b) 100kPa-loose sample with the whole domain considered. For each sample, two reference conformations belonging to the stationary regime are adopted. The deviatoric stress curve (dot-dashed line) and the volumetric strain (dashed line) curves are recalled in the subfigure of Loop 4. . . . . . 70

List of Figures xiii

4.10

The averaged PDFs of grain loops with respect to the incremental strain evolution (∆ε 2 ), over six microstructure conformations at different strain levels belonging to the critical regime (ε 2 = 0.0610, 0.0708, 0.0806, 0.0904, 0.1002 and 0.1100 for the 100kPa-dense sample and ε 2 = 0.0605, 0.0708, 0.0802, 0.0906, 0.1005 and 0.1106 for the 100kPa-loose sample). . . . . . . 4.16 Two characteristic strains (ε c1 and ε c2 ), the characteristic probability (P c ) and the ratio ε c1 ε c1 with respect to the mean stress at the critical state. The critical mean stress under different confining pressures (40 kPa, 100 kPa and 400 kPa) are 59 kPa, 143 kPa and 559 kPa, respectively. . . . . . . . . . .

4.17

Motivation and objectives

Granular materials are widespread in nature, like sands, rocks, rices etc., and in industry, for example, medicine powders [START_REF] Jaeger | Granular solids, liquids, and gases[END_REF]. Despite their prevalence, they present complex behavior even in the lab. Even though significant efforts have been devoted into understanding the mechanisms that govern deformation in granular materials subjected to compression and shear based on experiments and simulations, we still have a blurred vision about how the interactions of contacts and grains led to complex behavior in the system of granular materials [START_REF] Papadopoulos | Evolution of network architecture in a granular material under compression[END_REF].

Multiscale modeling of granular materials has received an increasing attention in the recent years based on notable progress in instrumentation technology and computing power. Multiscale modeling describes processes that aim to simulate continuum-scale behavior using information obtained from computational models of finer scales in the system, rather than resorting to empirical constitutive models, provided that the fine-scale physics and structure are better understood than those at a coarser scale [START_REF] Fish | Mesoscopic and multiscale modelling in materials[END_REF]. Even though this modelling method cannot replace conventional phenomenological paradigm currently in engineering which describes material behaviour directly at a coarser scale relevant for analysis and design and easier to compare to available experimental data rather than on microscale physics, it is believed that in future multiscale modelling could play a critical role, especially for large deformation problems [START_REF] Liang | Multiscale modeling of large deformation in geomechanics[END_REF].

Chapter 1. General introduction

In the macroscopic behavior in granular materials, one significant and attractive observation is the existence of stationary state where a multitude of conformations (local arrangements of particle structures) and conformational transitions exist for which the macroscopic stress and porosity remain constant under continuous shearing. Such macroscopic states is named as critical state in geomechanics. Critical state (CS), initially observed in geomaterials (e.g. soils) under biaxial or triaxial loading, form the focus of Critical State Theory (CST) [START_REF] Casagrande | Characteristics of cohesionless soils affecting the stability of slopes and earth fills[END_REF][START_REF] Roscoe | On the yielding of soils[END_REF][START_REF] Been | Critical state of sands[END_REF] which stands as a cornerstone of constitutive theory for granular materials. Constitutive laws built on CST are however phenomenological and cannot account for the underlying conformational transitions which are responsible for the emergence of critical states. To this respect, only multiscale models may capture the microscale detailed balance [START_REF] Sun | Experimental evidence of detailed balance in granular systems[END_REF]) responsible for critical state. But, many aspects of these transitions remain poorly understood, and no multiscale model exists in which the minimal set of conformational transitions responsible for critical state are embedded. Indeed, novel structures, transitions and regimes are still being uncovered in granular materials under shear, mostly thanks to the use of Discrete Element Modelling (DEM), as well as the development of lab equipment [START_REF] Rocks | Hidden topological structure of flow network functionality[END_REF].

Thus, it is intriguing to revisit critical state not only to improve fundamental understanding of stationary regimes that comprise large numbers of interacting units in complex systems [START_REF] Gardner | Connectance of large dynamic (cybernetic) systems: critical values for stability[END_REF][START_REF] May | Will a large complex system be stable?[END_REF]Papadopoulos et al. 2018), but also to find a method that can make the critical state as an emerging property in constitutive models.

This thesis focuses on understanding the mechanical responses of granular materials under shearing, with an emphasis on the micromechanics of critical state and on the emergence of CS in multiscale modelling based on a micromechanical-based constitutive model, following several questions:

• Can we observe the critical state in different loading paths, especially the dilatant or contractant proportional strain loading along which the volume never reaches a constant value?

• What is the micromechanical mechanism under the stationary state?

• How can we make the critical state emerge in micromechanics-based constitutive models without introducing empirical laws? Chapter 2 summarizes some of the relevant background literature. The discussion mainly revolves shearing behavior in granular materials, starting from the classically continuummechanical understanding and the description of the main features that are the most important to be considered in the development of a constitutive model for granular materials, and ending at the more recent explanation from the micromechanical perspective.

Chapter 3 offers a fresh perspective on the classic concept of critical state (CS) in granular materials by suggesting that CS can be defined through the use of a single proportional strain test. A comparison between proportional strain tests and biaxial tests simulated with DEM is presented. The mechanical responses and fabric metrics evolve along dilatant proportional strain loading paths according to similar values after the strain level has become large enough to wipe out the material memory in the homogeneous domains considered in this analysis, i.e., the shear band area in dense samples and the whole area in loose samples. This suggests that the micro-structure of a granular material subjected to proportional strain loading evolves with a constant ability to withstand shearing without volume change. Therefore, the CS concept can be generalized to a wide class of loading paths which shows that CS acts as a general attractor irrespective of the loading path considered (Deng et al. 2021c).

Chapter 4 examines the critical state in granular materials from a dynamical systems theory perspective. The lifespan and life expectancy of specific cluster conformations are tracked and quantified. Results suggest that these conformational clusters reorganize at similar rates in the critical state regime, depending on strain magnitudes and confining pressures. This rate of reorganization is quantified and it is found that the material memory rapidly fades, with an entirely new generation of force chains and grain loops replacing the old within a few percent strain. The chapter extends the results of the paper [START_REF] Deng | Lifespan dynamics of cluster conformations in stationary regimes in granular materials (in revision)[END_REF].

Chapter 5 presents a novel method to make critical state emerge in a micromechanicalbased constitutive model, namely the H-model (Nicot and Darve 2011b) by introducing a deactivation/reactivation procedure in local scale based on the results from Chapter 4.

The chapter starts with a review on H-model and finishes with preliminary stress-strain results predicted from the H-model with the novel method.

In conclusion, Chapter 6 summarizes the most salient contributions and findings of this thesis. Future lines of research related to this work are put forward. A brief literature review on the micromechanical investigations of granular materials is presented from two different standpoints, namely geomechanics and statistical physics, in this chapter. From a conventional geomechanical problem, the shearing behavior in granular materials is reviewed before being put into perspective with the recent micromechanical interpretations to the conventional problem (Section 2.1). The other line of approach starts with an introduction of Discrete Element Method (DEM) and statistical physics, and lands with some of their applications (Section 2.2-2.3). There will be some crossover between the two lines. In Section 2.4, a short review on multiscale modeling is presented. It may be seen as a statistical constitutive modeling approach applied to geomaterials.

Shearing behavior of granular materials

Granular materials under shearing exhibit complex behavior because of the collective rearrangement of particles. It is the discrete nature of granular materials that makes their 6 Chapter 2. A brief literature review behavior highly nonlinear, dissipative, and intimately linked to the microstructure. In the particular field of geomechanics, the complex behavior has been conceptualized and detailed conventionally through standard features such as hardening/softening, dilatancy/contractancy, critical state, etc., as shown in Fig. 2.1 in the framework of continuous mechanics. In addition to the features at the macro scale, the shearing behavior can be demonstrated by indexes at the micro/meso scale with a pair or a cluster of grains or contacts and their collective descriptors. Enormous efforts have been devoted to observing and interpreting these phenomena from macro to micro, from experiments to simulations, especially in geotechnical engineering, in order to find either the mechanisms of the specific features or appropriate descriptors that can feed constitutive modeling. In short, the aim of multiscale modeling is to narrow the gap between micro-scale characteristics and macro-scale features. In this section, reviewing work starts from the continuum-mechanical understanding and the description of the main features that are the most important to be considered in the development of a constitutive model for granular materials, and ends at the more recent explanation from the micromechanical perspective, ranging from detailed observations to abstract theories.

Dilatancy and contractancy

Dilatancy and contractancy are one of the most important properties of granular materials, defined by volume changes (increase and decrease, respectively) induced by shear deformation. The main lineage of understanding of dilatancy in the framework of continuum mechanics can be seen clearly thanks to the review work in [START_REF] Rowe | The stress-dilatancy relation for static equilibrium of an assembly of particles in contact[END_REF][START_REF] Houlsby | How the dilatancy of soils affects their behaviour[END_REF][START_REF] Li | Dilatancy for cohesionless soils[END_REF][START_REF] Guo | Stress dilatancy and fabric dependencies on sand behavior[END_REF]and Kruyt and Rothenburg 2016 at different ages. The concept of "dilatancy" was initially highlighted by Reynolds in 1885 Reynolds 1885. Dilatancy/contractancy plays an central role in constitutive modeling as it controls the non associated character of the flow rule in granular materials.

Empirically, the relationship of dilatancy with stress ratio, density (void ratio) and material internal state has been established in a list of theories, based on a series of experiments and motivated by soil behavior modeling. Taylor 1948 pointed out a link between shear strength and dilatancy properties of sands based on a hypothesis about the energy dissipation in frictional soils. The dilatancy equation used in Cam-Clay model by Schofield and Wroth 1968 for soil behavior is a direct analogy of Taylor's work hypothesis. Alternatively, Rowe 1962 proposed the stress-dilatancy theory which basically quantifies the effect of geometrical interlocking of the particles on the stress state. Following Rowe's theory, [START_REF] Bolton | The strength and dilatancy of sands[END_REF] gave a very simple empirical fit to friction and dilation angles according to a particularly comprehensive review of the experimental data. The dependency of dilatancy on initial stresses and density has also been exhibited experimentally and theoretically in literature [START_REF] Rowe | The stress-dilatancy relation for static equilibrium of an assembly of particles in contact[END_REF][START_REF] Wroth | A stress-strain relationship for the shearing behaviour of a sand[END_REF][START_REF] Matsuoka | Dilatancy characteristics of soil[END_REF]. It is expected that the denser the sand, the more it will tend to expand, but at a given density, the angle of friction reduces slightly with increasing stress level.

A physical estimation of dilatancy was initially conducted by Rowe 1962 based on uniform disks or rods in triangular close-packed arrangements, then on regular packing hexagons [START_REF] Li | Dilatancy for cohesionless soils[END_REF], which provides the physical basis of the relationship between the dilatancy and stress ratio. More recently, dilatancy has been treated as the emergence of a collective property within disordered packings. Numerical studies [START_REF] Peyneau | Frictionless bead packs have macroscopic friction, but no dilatancy[END_REF][START_REF] Azema | Internal friction and absence of dilatancy of packings of frictionless polygons[END_REF] and experiments [START_REF] Clavaud | Revealing the frictional transition in shear-thickening suspensions[END_REF] have indicated the absence of dilatancy in frictionless systems, which is different to the previous idea about "rearrangement" that dilatancy arises purely from geometric exclusive effect of hard particles. A work from [START_REF] Babu | Is friction essential for dilatancy and shear jamming in granular matter?[END_REF] investigated this paradox by comparing the dilatancy effect between amorphous and lattice assemblies, and established conditions under which dilatancy emerges naturally in frictionless sphere assemblies through jamming-unjamming processes.

In addition, recent experiments and DEM simulations of specimens suggest that significant levels of dilatancy emerge from highly coordinated particles motions, particularly within shear bands [START_REF] Oda | Microstructure of shear bands and its relation to the mechanisms of dilatancy and failure of dense granular soils[END_REF][START_REF] Tordesillas | Stress-dilatancy and force chain evolution[END_REF]; a micromechanical study of dilatancy conducted by Kruyt and Rothenburg (Kruyt and Rothenburg 2016) identified two microstructural mechanisms: dilatancy due to the deformation of mesostructures and dilatancy owing to topological changes in the contact network mainly governed by the parameter contact anisotropy and coordination number, respectively.

Chapter 2. A brief literature review

There are a number of constitutive models for granular materials that incorporate physical ideas about the evolution of granular microstructure within the continuum mechanics framework and accurately describe dilation properties. [START_REF] Guo | Stress dilatancy and fabric dependencies on sand behavior[END_REF] bypassed the minimization of the energy ratio and highlighted the importance of stress dilatancy and its microstructural dependence upon the behavior of sand along various stress/strain paths in axial symmetry conditions. Following the work of [START_REF] Guo | Stress dilatancy and fabric dependencies on sand behavior[END_REF][START_REF] Tordesillas | Stress-dilatancy and force chain evolution[END_REF] Tshaikiwsky 2011 characterized the evolution of stress-dilatancy with respect to the evolution of mesostructures and proposed a model incorporating not only fabric but also fabric evolution. Additional work in this direction has been carried out by [START_REF] Li | Modeling of dilative shear failure[END_REF]Dafalias 2012 and[START_REF] Gao | A critical state sand plasticity model accounting for fabric evolution[END_REF], where the link between dilatancy and microstructure is embedded within the flow rule in the constitutive model.

Softening and hardening

Strain softening and hardening are concepts originally developed for metals in plasticity theory. During plastic deformation the loading surface is said to harden (the yield limit increases in the stress space) or soften (the yield limit decreases in the stress space). In triaxial or biaxial compression tests of granular materials, strain hardening and softening phenomena are exhibited as the increase and decrease in deviatoric stress, coming along with dilatancy or contractancy normally demonstrated by the incremental volumetric strain. Particularly, the strain softening is a complex behavior, characterized by a gradual loss of shear resistance with strain after a peak strength has been reached [START_REF] Prévost | Soil mechanics and plasticity analysis of strain softening[END_REF][START_REF] Read | Strain softening of rock, soil and concrete-A review article[END_REF]. Discussions on weather the strain softening is an intrinsic property of granular materials or a structure response have never stopped.

On the one hand, it has been viewed as an inherently material property of granular materials and routinely contained in constitutive models, since it is commonly reported in the experiments of dense materials [START_REF] Chu | Strain-softening behavior of granular soil in strain-path testing[END_REF][START_REF] Verdugo | The steady state of sandy soils[END_REF]. For example, Drucker, Gibson, and Henkel 1957 introduced the concept of soil as a workhardening material in the extended Mohr-Coulomb model. A link between current peak strength and a state parameter which is a combination of volumetric and mean effective stress information has been demonstrated to reproduce strain hardening in [START_REF] Wood | Strain softening and state parameter for sand modelling[END_REF]Belkheir 1994. Yao, Hou, and[START_REF] Yao | UH model: three-dimensional unified hardening model for overconsolidated clays[END_REF] proposed a unified hardening parameter that is independent of stress paths in order to incorporate strain hardening and softening based on the modified Cam-Clay model. This hardening parameter is used in constitutive models for sands [START_REF] Yao | UH model for sands[END_REF][START_REF] Yao | Unified hardening (UH) model for clays and sands[END_REF].

On the other hand, it has been argued that the observed strain softening in soils is not a material property but mainly the consequence of inhomogeneity in the deformation field of the specimens during the loading, as pointed out by Read and Hegemier 1984 after reviewing an extensive set of experiments to discuss the strain softening for rock, soil and concrete. It has been demonstrated that the strain softening is strongly affected by the size and shape of the sample prepared, and is normally accompanied by strain localization [START_REF] Fu | Fabric evolution within shear bands of granular materials and its relation to critical state theory[END_REF][START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]. Also, some researchers advocate that both material softening and structural softening exist but appear under different conditions [START_REF] Sterpi | An analysis of geotechnical problems involving strain softening effects[END_REF]Liu et al. 2020a). For example, softening along with strain localization in biaxial tests is regarded as structural softening, while that obtained in homogeneous samples as material softening.

It is difficult to close the debate within the framework of plasticity theory, as experiments with perfectly homogeneous material does not exist, especially in granular materials. Various modes of bifurcation are possible and do actually develop during the experiments, as pointed out by [START_REF] Hettler | Behaviour of dry sand tested in a large triaxial apparatus[END_REF] More recent micromechanical investigations have helped to some extent by further understanding of the evolution of fabric and stress, and thereby by providing new insight into strain hardening and softening (Li and Dafalias 2012;[START_REF] Tordesillas | Transition dynamics and magic-number-like behavior of frictional granular clusters[END_REF]Liu et al. 2020a). They have been characterized by some statistical physics descriptors related to the rearrangement of grains. For instance, Liu et al. 2020a proposed mesoscopically-based framework to interpret both hardening and softening mechanisms; Tordesillas 2007 explained the hardening and softening process by jamming-unjamming theory; to name it. Thus, it might be asked: since strain hardening partly induced by the self-reorganization of particles is recognized as a material property and involved into constitutive models, why not softening that is also induced by rearrangement of particles? The hardening is the process to be more well-organized and easier to be described than the softening process where the structures are much more chaotic and difficult to grasp with higher heterogeneity within the contact network. This could be a main difficulty in describing the softening process. In addition, it has been widely accepted that the inhomogeneity property should be considered carefully in constitutive modeling. Anisotropy builds quickly as it relies mostly on the contact network. If the anisotropy of contacts can be regarded as an inherent property in granular materials, it would be important to incorporate strain softening in constitutive models from a micromechanical viewpoint.

Critical State

Critical state (CS), initially proposed by [START_REF] Casagrande | Characteristics of cohesionless soils affecting the stability of slopes and earth fills[END_REF], is a physical phenomenon observed from a macromechanical perspective, in the form of a stationary state where stress and volume tend to be constant under continuous shear strain. Dense granular materials reach the critical state as a result of decreasing in dilatancy normally with shear band generation, while loose materials tend to reach the same state after decreasing in contractancy. When the steady state is reached, deformation continues without volume changes. The classic critical state theory (CST) developed by [START_REF] Roscoe | On the yielding of soils[END_REF]Schofield and Wroth 1968 is the basic principle behind critical state soil models [START_REF] Been | A state parameter for sands[END_REF][START_REF] Been | Critical state of sands[END_REF][START_REF] Yao | UH model: three-dimensional unified hardening model for overconsolidated clays[END_REF][START_REF] Li | Dissipation consistent fabric tensor definition from DEM to continuum for granular media[END_REF], which implies two conditions: a steady stress ratio and a steady void ratio. They are expressed analytically by

η = η c = (q/p) c e = e c = êc (p) (2.1)
with e c = êc (p) the critical void ratio which defines the Critical State Line (CSL) in the e -p plane.

More recently, it has been proven that not only do both stress and density converge to a steady state at the critical state, but that this convergence and stationarity also manifest at the micro-scale [START_REF] Drescher | Photoelastic verification of a mechanical model for the flow of a granular material[END_REF][START_REF] Fu | Fabric evolution within shear bands of granular materials and its relation to critical state theory[END_REF][START_REF] Kruyt | On micromechanical characteristics of the critical state of two-dimensional granular materials[END_REF]Fu and Dafalias 2015;Kruyt and Rothenburg 2016;[START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF][START_REF] Kawamoto | All you need is shape: predicting shear banding in sand with LS-DEM[END_REF]. For instance, [START_REF] Rothenburg | Critical state and evolution of coordination number in simulated granular materials[END_REF] developed a relation between the critical state and coordination number based on DEM simulations; Fu and Dafalias 2011 studied the evolution of fabrics at the critical state using DEM simulations; Kruyt and Rothenburg 2014 re-investigated the relationship between the continuum features at the macroscale and micromechanical quantities including coordination number and fabric anisotropy at both the particle and interparticle contact level and the grain loop scale; Kuhn 2016b demonstrated the stationarity at the critical state as well as the evolution of their disorder toward the critical state based on several characteristics at the microscale; Zhu et al. 2016 investigated different mechanical features of two kinds of meso-structures (force chains and grain loops) and proved that the two failure modes (localized and diffuse) are homological with respect to the concept of the critical state.

According to several micromechanical investigations, it can be shown that the difference between the current porosity and the CS porosity does not in itself determine the evolution of the micro-structure. The coaxiality between the plastic strain rate direction and the fabric anisotropy proved to be a relevant state variable in defining critical state (Li and Dafalias 2012;Theocharis et al. 2019). A third condition that quantifies the role of fabric anisotropy in terms of its intensity and its relative orientation with respect to loading direction should also be taken into account, in addition to the aforementioned two conditions of constant stress ratio and void ratio in the anisotropic critical state theory (ACST) (Li and Dafalias 2012;Theocharis et al. 2019). In ACST, a fabric-related tensor F is adopted. To incorporate micromechanical observations into a macroscopic continuum mechanics description, the norm F ≥ 0 and the unit-norm n F of F are considered with:

F = F n F , F = √ F : F, n F : n F = 1, trn F = 0 (2.2)
The Fabric Anisotropy Variable A is introduced by

A = F : n = F n F : n = F N (2.3)
where n represents a direction along which the loading is applied, say, the direction of plastic flow; thus, it can be related to the stress tensor for monotonic radial loading. N = n F : n is a measurement of the relative orientation of F and n. A tends toward 1 at critical state because both F and N tend toward 1. Then, ACST can be expressed as follows:

η = η c , e = e c = êc (p), A = A c = 1 (2.4)
2.2 Discrete element view of granular materials

Discrete Element Method

Discrete mechanics not only helps to prove or correct the transitional understanding and empirical theories of granular materials (as explained in the previous section), but also brings new knowledge enriching our interpretation of the complex behavior of granular materials under shearing. There are two types of tools used in micromechanical analysis: experiments and numerical simulations. As for the former, the technique of particle image velocimetry (PIV), computer tomography, magnetic resonance imaging, and confocal microscopy can provide rich microscopic information of the granular material, like the shape of particles as shown in Fig. 2.3 [START_REF] Kawamoto | All you need is shape: predicting shear banding in sand with LS-DEM[END_REF], local void ratio (), local strain field and the geometric information of interparticle contacts and so on; in terms of mechanical data, photo-elastic rods are the most popular granular assemblies used to study the force distribution and transmission within a deformed packing, as shown in Fig. 2.2. The latter refers to the discrete element method simulation is used in this PhD thesis to run numerical experiments with unlimited access to micro data. Discrete element method (DEM), developed by [START_REF] Cundall | A computer model for simulating progressive, large-scale movement in blocky rock system[END_REF][START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], has been widely accepted as a powerful method to simulate the complex mechanical behavior of granular materials at the Representative Elementary Volume (REV) scale. Note that REV proposed by [START_REF] Hill | Elastic properties of reinforced solids: some theoretical principles[END_REF] is the smallest volume providing representative mechanical quantities at macroscale, like stress and strain, in the framework of continuum mechanics. DEM simulations lay out abundant details at the microscopic level, from geometric data to force and energy information, which is relative difficult to be accessed in concrete experiments. DEM has been regarded as a reliable tool to find the connection between local physics and macro properties of granular materials in considerable investigations including those reviewed in the previous section. In this subsection, Discrete Element Method (DEM) is briefly reviewed.

All DEM simulations performed in this PhD thesis are based on the open-source DEM software YADE (Yet Another Dynamic Engine), which is developed with C++ and Python programming languages (Smilauer et al. 2015).

Computation loop

In DEM, as introduced in the pioneering work [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], the interaction between particles is viewed as a transient problem with states of equilibrium developing whenever the internal forces balance. The numerical scheme is explicit. The equilibrium contact forces and displacements of a stressed assembly of non-deformable discs are found through a series of calculations tracing the movements of the individual particles, and the movements are the result of the propagation through the medium of disturbances originating at the boundaries: this is therefore a dynamic process. The dynamic process is described numerically. It is assumed that the time step chosen should be small enough in order that disturbance cannot propagate from any disc further than its immediate neighbours during a single time step. The calculations performed in DEM alternate between the application of Newton's second law to the discs and a constitutive law at the contacts. Newton's second law gives the motion of a particle resulting from the forces acting on it. The constitutive law is used to compute contact forces from displacements. Considering a granular material at time t 0 , geometrical data describing the positions and size of grains and walls, and physical properties of the grains are given. To determine the contact forces and positions of grains at time t 0 + ∆t, there are four steps as follows:

• Detect the contacts: two spherical grains are in contact only if the distance between centers of two particles is less than the sum of their radii. If this condition is met, go to the next step.

• Compute the contact forces based on the contact law and relative displacement increments.

• Integrate the force and momentum equations for each grain, and calculate the acceleration of them according to Newton's second law.

• Update the positions of the grains by integrating the acceleration during ∆t.

After that, the particle positions are updated at this time step, and become the initial state of the next time step. Then, repeat the circulation under control.

Contact law

The contact law is one important modeling hypothesis in DEM simulations since it affects the macroscopic behavior significantly together with particle arrangement. In DEM, some contact laws use strains and stresses while others are expressed through displacement -force relation. The contact law presented here is the most common and classic one in DEM, originally proposed by [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] as the simplest non-cohesive elastic-frictional contact model. As shown in Fig. 2.4, when a new contact is established, it is assumed that there are two linear springs awakened: the one in the normal direction and the other in the tangential direction to the contact, characterized by the normal and tangential stiffness k n and k t respectively. A Coulomb-type friction law is incorporated by introducing the intergranular friction angle ϕ. Consequently, the normal contact force F n is defined by the overlapping distance u n between the two particles and the normal contact stiffness k n , expressed as Equation 2.5; the tangential contact force F t after being assessed by the friction law can be written as Equation 2.6. Note that the stiffness k n depends on a material stiffness E, and the radii of the two particles R 1 and R 2 , and k t is determined from k n and a stiffness ratio α.

F n = k n u n (2.5) F t = k t u t if |F t | < F n tanϕ tanϕF n otherwise (2.6) k n = E 2R 1 R 2 R 1 + R 2 (2.7)

Micromechanical analysis in granular materials

As assumed in DEM simulations, Granular materials can be regarded as assemblies of rigid particles, where particles interact through contact forces and small contact areas. Thus, the microscale level is that of particles and contacts [START_REF] Kruyt | On micromechanical characteristics of the critical state of two-dimensional granular materials[END_REF], as shown in the Fig. 2.4. Granular materials can be decomposed into weak network and strong network according to the force transmission, or grain loops tessellating the 2D assembly, and many other types of mesoscopic cluster. In summary, the structural object of granular materials can be grouped through different scales: (i) microscale at the particle and contact level, (ii) mesoscale involving a cluster of grains, (iii) macroscale of the global contact network.

Related to these diverse structures, variables can be constructed from various aspects (Kuhn 2016b). Information of the state (being) and transition (becoming) can be involved in to the following three perspectives: (i) configuration quantities, such as coordinating number, contact orientation; (ii) force quantities, such as particle or cluster stresses and contact forces; (iii) movement quantities, such as particle translations, local strains. The first two aspects are spacial variables associated with the status of granular assembly. The third set of attributes is space and time related, comprising transitional or rotational rate quantities that are driven by the bulk deformation and depending on the loading direction.

Particle/contact

The microscale level in granular materials is that of particles and contacts. The relative displacement and the contact force are information of interest [START_REF] Kruyt | On micromechanical characteristics of the critical state of two-dimensional granular materials[END_REF][START_REF] Imole | Micro-macro correlations and anisotropy in granular assemblies under uniaxial loading and unloading[END_REF]. The coordination number, the fabric tensor and the distribution function of orientations of contacts or particles etc. are the most popular statistics based on this information used in multiple analysis.

The coordination number is the average number of contacts per particle, and refers to the compact degree (closely linked to the density and porosity) of the material [START_REF] Kruyt | On micromechanical characteristics of the critical state of two-dimensional granular materials[END_REF]. It is defined as

Z = 2N c N p (2.8)
where N c and N p represent the number of contacts and particles, respectively. Fabric tensor collects the orientation information in contacts or particles. Different types of fabric tensors can be constructed with different vectors [START_REF] Oda | Deformation mechanism of sand in triaxial compression tests[END_REF][START_REF] Satake | Fabric tensor in granular materials[END_REF][START_REF] Fu | Fabric evolution within shear bands of granular materials and its relation to critical state theory[END_REF]. For example, a fabric tensor of contact normals n c i can be expressed by

F ij = 1 N c c∈C n c i n c j (2.9)
where the sum is over all contacts c in the set of contacts C. Similar quantities can be constructed with the direction of contact force for instance.

Force chains

Force chains, a concept of quasilinear particle clusters where stress is concentrated, have attracted widespread interest. There is abundant experimental and numerical evidence that the spacial distribution and temporal evolution of the force structures on the mesoscale govern macroscopic behavior. Particularly, the static state in the sample and the buckling (or bending) under loading of force chains play the primary role in carrying external load and energy dissipation, respectively [START_REF] Drescher | Photoelastic verification of a mechanical model for the flow of a granular material[END_REF][START_REF] Radjai | Force distributions in dense two-dimensional granular systems[END_REF][START_REF] Peters | Characterization of force chains in granular material[END_REF][START_REF] Tordesillas | Force chain buckling, unjamming transitions and shear banding in dense granular assemblies[END_REF][START_REF] Tordesillas | Buckling force chains in dense granular assemblies: physical and numerical experiments[END_REF][START_REF] Tordesillas | Force cycles and force chains[END_REF][START_REF] Saitoh | A master equation for the probability distribution functions of forces in soft particle packings[END_REF]Wautier, Bonelli, and Nicot 2018a).

An objective description of what constitutes a force chain initially proposed by [START_REF] Peters | Characterization of force chains in granular material[END_REF] is illustrated in Figure 2.5 and summarized as follows:

• The number of contacting particles in a force chain is not less than 3. • The major principal stress on each particle belonging to a force chain exceeds the global average value.

• The major principal stress direction is aligned with the contact direction, for example, with less than 45 • deviation.

The stress tensor on each particle is computed according to the summation running over all external forces of the particle, as shown in Equation 2.10 in the initial method [START_REF] Peters | Characterization of force chains in granular material[END_REF]. A modified approach proposed in Muthuswamy and Tordesillas 2006 adopted a force moment tensor replacing the stress tensor for each particle calculated as in Equation 2.11 to avoid preference for small particles.

σ ij = 1 V N c=1 f c i r c j (2.10)
where V is the volume of the particle, N is the number of contacts of the particle, f c i is the ith component of the force acting at the contact, and r c j the jth component of the radius vector from the center of the particle to the point of contact.

S ij = 1 r N c=1 f c i r c j (2.11)
where r is the particle radius.

In short, there are three thresholds needed to detect a force chain: magnitude of the principal stress (or force), direction of the principal stress (or force), and the minimum number of particles.

As for the force chain itself formed in a specimen, it is affected by multiple factors. Findings from [START_REF] Muthuswamy | How do interparticle contact friction, packing density and degree of polydispersity affect force propagation in particulate assemblies?[END_REF] indicate that increasing interparticle friction, packing density and degree of polydispersity promotes the formation of straighter chains and a greater degree of branching in the force chain network. Denser systems tend to favour shorter chains. Thus, straighter and shorter force chains, combined with a greater degree of branching in the force chain network, result in a macroscopically stronger granular material.

The failure mode of force chains was investigated numerically and experimentally by [START_REF] Tordesillas | Buckling force chains in dense granular assemblies: physical and numerical experiments[END_REF]. It is demonstrated that frictional rolling is the predominant mode of contact failure in a force chain buckling, and in the contacts with and within its laterally supporting neighbours.

Most importantly, connections of force chains of discrete medium to the formation of the most compressive principal stress, and of force chains buckling to energy dissipation, unjamming transition and global failure have also been looked into carefully in literature [START_REF] Oda | Microstructure of shear bands and its relation to the mechanisms of dilatancy and failure of dense granular soils[END_REF][START_REF] Radjai | Bimodal character of stress transmission in granular packings[END_REF][START_REF] Tordesillas | Force chain buckling, unjamming transitions and shear banding in dense granular assemblies[END_REF][START_REF] Tordesillas | Buckling force chains in dense granular assemblies: physical and numerical experiments[END_REF][START_REF] Tordesillas | Stress-dilatancy and force chain evolution[END_REF][START_REF] Luding | Contact force measurements and stress-induced anisotropy in granular materials[END_REF]. Force chain buckling has been regarded as the characteristic of unjamming events at the mesoscale, to be responsible to non-affine strain, and the mechanism for energy dissipation [START_REF] Tordesillas | Force chain buckling, unjamming transitions and shear banding in dense granular assemblies[END_REF]. It has been attributed the governing mechanism for shear banding to force chain buckling that is concentrated inside shear bands, and responsible for the high levels of dilatancy observed within these bands [START_REF] Tordesillas | Force chain buckling, unjamming transitions and shear banding in dense granular assemblies[END_REF][START_REF] Tordesillas | Buckling force chains in dense granular assemblies: physical and numerical experiments[END_REF][START_REF] Tordesillas | Stress-dilatancy and force chain evolution[END_REF].

Grain loops

Grains loops initially attracted attention for computing local strain [START_REF] Satake | A discrete-mechanical approach to granular materials[END_REF][START_REF] Kruyt | Micromechanical definition of the strain tensor for granular materials[END_REF][START_REF] Kuhn | The critical state of granular media: convergence, stationarity and disorder[END_REF][START_REF] Kruyt | Statics and kinematics of discrete Cosserat-type granular materials[END_REF]. The geometric structure of a 2D granular assembly can be described by a particle graph which entirely subdivides the assembly into polygonal sub-domains, namely grain loops (also known as mesoloops or mesodomains), as shown in Figure 2.6. Thus, grain loops, consisting of contacting grains, are an important type of mesostructure in 2D in granular materials in terms of local strains. Note that force cycles mentioned once in the previous subsection is an example of the grain loops defined here.

The stress and strain tensors of each grain loop are defined from the contact forces and particle movements by [START_REF] Nguyen | Analysis of structure and strain at the meso-scale in 2D granular materials[END_REF][START_REF] Nguyen | Local stress analysis in granular materials at a mesoscale[END_REF]. The elongation degree and elongation direction defined based on the loop fabric tensor have been used in [START_REF] Nguyen | Analysis of structure and strain at the meso-scale in 2D granular materials[END_REF] to investigate the relation between structure and strain at the mesoscale. It has been revealed that strain in grain loops significantly depends on the elongation degree and the orientation of loops. It was shown that contractancy takes place within the grain loops which are elongated along the minor principal compression direction whereas dilatancy takes place within the grain loops which are elongated along the major principal compression orientation. In addition, the structure and stress on grain loops have been connected by [START_REF] Nguyen | Local stress analysis in granular materials at a mesoscale[END_REF] It is demonstrated that the distribution of the local mean stress is uniform in any grain loop, while the distribution of the deviatoric stress is significantly dependent on the elongation direction and on the elongation degree of the grain loops. The local stress ratio is higher within the grain loops that are elongated in the global compression direction than that within the ones elongated in the global extension direction.

In addition, grain loops at the critical state where the specimen volume keeps constant under shearing have been investigated by [START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF][START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF][START_REF] Kruyt | On micromechanical characteristics of the critical state of two-dimensional granular materials[END_REF]. The results from [START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF][START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF] demonstrate that the percentage of each loop category (grouped according to the edge number) reaches a steady state that corresponds to the critical state regime at the specimen scale. Consistent to this finding, the results on the topological exchange within grain loops along a biaxial test show that difference between the new and disappearing loops follow a non-monotonic evolution, until fluctuating around zero at the critical state for the different categories of grain loops [START_REF] Liu | Sustainability of internal structures during shear band forming in 2D granular materials[END_REF]. At a specific state belonging to the critical regime, the distributions of loop orientation characterized by the loop fabric tensor are identical for each category independent of the coordination number [START_REF] Kruyt | On micromechanical characteristics of the critical state of two-dimensional granular materials[END_REF]. The interplay between force chains and force cycles (N cycles whose contacts each bear above average force) has been examined by [START_REF] Tordesillas | Force cycles and force chains[END_REF] It is unraveled that the three-force cycles play an primary role in reinforcing the force chains and indicating failure resulting from force chain buckling. The significant role of grain loops in supporting force chains in the framework of force networks is another reason for considering grain loops (normally referred to as force cycles instead) [START_REF] Tordesillas | Force cycles and force chains[END_REF].

Multiscale modeling of granular materials

A multiscale modeling of granular materials involves a trade-off between precision and simplicity in describing the microstructure. One should consider at which scale to simplify the physical information because of the existence of multiple spatial and temporal scales. It follows that studying structures at multiple length scales offers a fruitful pathway for the understanding of static and dynamic properties of granular materials. Physical information can be incorporated in different ways at different length scales.

The multiscale framework of coupling Finite Element Method (FEM) and Discrete Element Method (DEM) is promising method to deal with small scale problems. It takes advantage of the efficiency of FEM at solving boundary value problems at structural level and the capability of DEM to capture complex material behaviors [START_REF] Meier | On the multiscale computation of confined granular media[END_REF][START_REF] Guo | Multiscale characterization of the shear behavior of granular media[END_REF][START_REF] Shahin | A study of the influence of REV variability in double-scale FEM× DEM analysis[END_REF][START_REF] Desrues | From discrete to continuum modelling of boundary value problems in geomechanics: An integrated FEM-DEM approach[END_REF]. This approach can reproduce interesting experimental observations such as the inception of shear band under smooth symmetric boundary conditions, large dilation, critical state reached within the shear band, etc. Nevertheless, these DEM-FEM methods have a limitation when geotechnical-engineering problems are considered due to the fact that a great number of particles need to be contained in a boundary value problem (BVP) in order to reach local convergence. It makes these methods waste computational resources resulting in inefficient simulations. In addition, all the information contained in the DEM simulation is not used because of averaging, which makes the approach not very economical in terms of data saved in memory.

Micromechanically-based constitutive models should be a good alternative to couple with a numerical methods e.g., the finite element method (FEM) [START_REF] Xiong | A multiscale work-analysis approach for geotechnical structures[END_REF], the smooth particle hydrodynamics (SPH) method [START_REF] Xiong | A novel multi-scale large deformation approach for modelling of granular collapse[END_REF], and Finite Difference Method (FDM) (Wautier et al. 2021). There are different types of micromechanical constitutive models (Li and Dafalias 2012;[START_REF] Chang | An elasto-plastic model for granular materials with microstructural consideration[END_REF]Nicot and Darve 2011b). Some of them incorporate statistical physics into the transitional phenomenal laws (Li and Dafalias 2012), the other ones construct a model from geometrical features and try to recover the mechanical behavior at macroscopic scale explicitly, for example the H-model (Nicot and Darve 2011b).

It is now recognized that the relative displacement of any contact is not proportional to the macroscopic strain. In the H-model (Nicot and Darve 2011b), one of most important hypotheses is that this approach describes the local behavior through mesostructures instead of contacts. One main advantage of such an approach is to give access to rich microstructural information that is out of reach for standard phenomenological models in which the microstructure is simplistically embedded as internal variables, and is not described in an explicit way. Such approaches try to use a minimal set of data available in DEM simulations. It is more computationally affordable than DEM.

Chapter 3

On the attraction power of critical state Although the structure of a granular material appears to be simple at the microscopic scale, its behavior shows itself to be complex at the macroscopic scale, mostly due to the collective rearrangement of particles. Based on the gain or loss of contacts, granular materials will easily adjust to any change under loading conditions. Among the accessible microstructures, some have the ability to withstand constant shearing with no change of volume. A state with such microstructures is known as critical state (CS) and plays a leading role in the most popular constitutive relations for granular materials, as reviewed in the second chapter.

Research on CS is nearly exclusively based on conventional tests (triaxial tests in 3D

or biaxial tests in 2D) in which a constant lateral stress condition dσ lateral = 0 is imposed together with a constant strain rate in one given direction dε axial = constant [START_REF] Been | A state parameter for sands[END_REF][START_REF] Been | Critical state of sands[END_REF][START_REF] Fu | Fabric evolution within shear bands of granular materials and its relation to critical state theory[END_REF]Fu and Dafalias 24 Chapter 3. On the attraction power of critical state 2015; [START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]. This preferential selection is linked to the tendency of generating constitutive models to work with stress and to calculate the resulting strains, as well as the difficulty of imposing fully controlled strain paths in the laboratory [START_REF] Ibraim | Strain path controlled shear tests on an analogue granular material[END_REF][START_REF] Daouadji | Phase transformation states of loose and dense granular materials under proportional strain loading[END_REF]. For conventional loading, the p -q path is preset; samples are loaded by increasing the axial incremental strain under constant lateral confining pressure; volumetric strain and deviatoric stress responses are recorded according to the loading path. This test allows for the existence of stationary states in volume and stress (the so-called critical state). Proportional strain tests, on the other hand, have been much less studied whereby the loading is imposed by proportional strain rates in the axial and lateral directions with dε lateral = λdε axial and dε axial = constant. Under proportional strain loading condition, the volumetric strain is no longer a response but a loading variable and the stress path is not known beforehand, which makes the interpretation of such tests more complex than it is for triaxial or biaxial tests. Except for the particular case of the undrained triaxial test, whereby the volume is kept constant, the continuous change in volume prevents stationary states in volume and stress along proportional strain paths to be observed. Without the existence of a stationary state, investigating the evolution of the microstructures along a proportional strain loading path becomes a major challenge. In other words, tracking the relation between triaxial (or biaxial) and proportional strain loading paths in terms of mechanical responses and underlying microstructures is an open and stimulating topic.

In order to generalize the concept of CS, we focus on the relation between the material responses under dilatant proportional strain loading paths and biaxial loading paths in 2D. These two types of loading paths, as well as mixed biaxial/proportional strain paths, have been simulated by DEM. The details of the DEM simulation are shown in Section 3.1. Mechanical responses along biaxial tests and proportional strain paths are presented in Section 3.2. The relation between proportional strain tests and biaxial tests in the p -q -e space is analysed in Section 3.3. In Section 3.4, the fabric-related CS locus is studied. The mechanical responses along the mixed loading paths are presented in Section 3.5, as is the relation between the imposed dilatancy in proportional strain tests and the dilatancy in biaxial tests.

DEM simulation

The Discrete Element Method (DEM) [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] is a powerful numerical method to simulate the global behavior of a set of grains interacting through contact laws. It has been widely used to simulate the mechanical response of granular assemblies under 3.1. DEM simulation various loading conditions. The mechanical states are characterized locally by kinematic information including position, rotation and velocity of grains, as well as static information based on contact force between contacting particles. The open source software YADE (Smilauer et al. 2015) has been used.

Under consideration here is a quasi-2D soil sample in the form of an assembly of a single layer of 20,000 spherical particles contained within a surface domain of 1 m × 1.5 m. The particle sizes have an average diameter d 50 = 0.008 m and d max /d min = 2. The elastofrictional law introduced by Cundall and Strack 1979 has been adopted as the contact law. The cohesiveless contact parameters between two grains contain a normal and tangential linear spring of respective stiffness k n and k t , as well as a friction characterized by a friction angle ϕ = 35 • . k n /D s is given as 300 MPa, where

D s = 2R 1 R 2 /(R 1 + R 2 )
and R 1 , R 2 are the radii of particles in a given contact and k t /k n is 0.5. If the pressure of 100 kPa is considered, the average ratio of contact overlap and particle size < u n > /d 50 is around 0.16%, which is close to the value (10 -3 ) in literature (Kruyt and Rothenburg 2016;Theocharis et al. 2019). For sample preparation, two schemes can be used to obtain an isotropic compression of the specimen: the boundary moving scheme and the internal compacting scheme (Wautier, Bonelli, and Nicot 2019a). In this study, particles are enlarged, at first, while keeping the boundary walls fixed to generate a preliminary sample up to 90 kPa. Then, to achieve a more precise consolidation pressure, the sample undergoes an isotropic consolidation by imposing an equivalent incremental strain on the boundaries in the vertical and lateral directions up to a specific consolidation state. During this preparation process, contact friction angles of 2 • and 35 • are used to obtain dense and loose samples, respectively. Consolidation pressures of 10 kPa, 20 kPa, 40 kPa, 60 kPa and 100 kPa are considered. Three dense samples and five loose samples are thus prepared and named as indicated in Table 3.1. The void ratio is accounted for, based on solid and void surfaces on the 2D plane as shown in Fig. 3.1. During the shearing process, a contact friction angle ϕ = 35 • is adopted in all the samples. Note that the soil mechanics convention is adopted throughout the thesis with compression and contraction counted positive. As a result, dilatant/contracting volumetric strain are considered as negative/positive respectively.

For proportional strain tests, 100kPa-dense sample is used with an initial void ratio of 0.191. In the 2D simulation, dε 1 = λdε 2 and dε 2 = 10 -2 s -1 . As dε v = (λ + 1)dε 2 , λ > -1 corresponds to a contracting test, λ = -1 to a constant volume test, and λ < -1 to a dilatant test. We have adopted λ = -1.2, -1.3, -1.4 for the three dilatant proportional strain tests. Biaxial loading paths with corresponding pressures are applied to all the samples indicated in Table 3.1. The stress and strain states are then described in 2D as follows: deviatoric stress q = σ 2 -σ 1 , mean stress p = (σ 1 + σ 2 )/2, and volumetric strain ε v = ε 1 + ε 2 , where σ 1 and σ 2 are the principal stresses at horizontal and vertical directions and ε 1 and ε 2 are the principal strains. The direction 2 is always the vertical loading control direction ( ε2 = 10 -2 s -1 ), as shown in Fig. 3.1. When the pressure ranges from 20 kPa to 200 kPa, the Inertial Number I of the granular system ranges among 2.77 × 10 -5 -8.76 × 10 -5 , which corresponds to a quasi-static state (Da [START_REF] Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF].

In reality a 2D model of a granular material is quantitatively different from sand assemblies; for example void ratio values and coordination numbers in 2D assemblies are smaller than those in 3D [START_REF] Fu | Fabric evolution within shear bands of granular materials and its relation to critical state theory[END_REF], and stresses expressed in kPa relies on the use of an arbitrary out of plane dimension 1 . However, given that the complexity of the constitutive behaviour of granular assemblies stems mainly from the local properties and the disordered packing, both effects can be captured in 2D simulations (Nicot and Darve 2011b). Qualitative investigations of granular materials based on DEM simulations in 2D can then be considered as both effective and efficient. Proof of such effectiveness and efficiency can be found in various studies of the microscopic mechanism behind the mechanical responses of granular materials under different loading paths have been investigated [START_REF] Fu | Fabric evolution within shear bands of granular materials and its relation to critical state theory[END_REF][START_REF] Kruyt | On micromechanical characteristics of the critical state of two-dimensional granular materials[END_REF]Fu and Dafalias 2015;[START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]Liu et al. 2020a). Thus, the 2D assembly simulation is adequate enough for qualitative investigations of critical states under proportional strain and biaxial loadings performed in this study.

Mechanical response

3.2.1 Mechanical response along dilatant proportional strain loading paths

Stress-strain analysis

Figure 3.2 shows the evolution of the deviatoric stress q along the axial strain ε 2 under dilatant proportional strain loading path with λ = -1.2, -1.3, -1.4. In this figure, it can be observed that the deviatoric stress q grows quickly at first and then decreases after the peak. Gradually, q decreases to 0. When λ = -1.4, the sample undergoes the fastest dilatancy and the corresponding stress curve reaches the earliest peak at ε 2 = 0.013 and the earliest zero deviatoric stress at ε 2 = 0.061. When λ = -1.3 and λ = -1.2, it appears that the peaks are reached at ε 2 = 0.017 and ε 2 = 0.019 because of smoother and smoother dilatancy. Zero mean pressure p occurs at ε 2 = 0.085 and ε 2 = 0.142, respectively.

Kinematic pattern: from localization to diffuse

The existence of a stress peak (Fig. 3.2) corresponds to a generalized limit state (Nicot and Darve 2011a) which has been shown to be a proper failure state (Wan, Nicot, and Chapter 3. On the attraction power of critical state Darve 2017). Various failure modes, characterized by localized or diffuse patterns, can be encountered after the peak, along the descending branch. In biaxial tests, either diffuse or localized modes will appear (Nicot and Darve 2011a). In this study, however, these two modes appear successively along each proportional strain test. In this subsection, the incremental deviatoric strain distribution has been used to characterize the kinematic pattern, as introduced in literature (O'Sullivan, D. Bray, and [START_REF] O'sullivan | A new approach for calculating strain for particulate media[END_REF][START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]. For example, in the test with λ = -1.2, a typical diagonal shear band traversing the whole specimen appears at ε 2 = 0.008, as shown in Fig. 3.3 (I), corresponding to the point I in Fig. 3.2. After a process of dilatancy, the shear band begins to vanish in the strain state ε 2 = 0.110 shown in Fig. 3.3 (II). Fig. 3.3 (III) demonstrates the kinematic pattern when ε 2 = 0.142, where the mean pressure p is close to zero. An evolution from localization, when the strain largely concentrates in a partial domain of the material, to liquefaction, when the effective stress within the granular specimen is reduced to essentially zero, has been observed along the three dilatant proportional strain paths, as shown in Fig. 3.3. With the increase in volume along the loading path (λ < -1), the density of the sample gradually decreases and the kinematic pattern changes from localized to diffuse until a zero mean pressure p is reached.

Kinematic patterns characterize the nature of the failure mode, together with the spatial domain. It has been widely accepted that fabric-related measures should be taken either within the shear band or within the full sample when the shear deformation is diffuse, since measures within nonhomogeneous domains have no constitutive meaning [START_REF] Desrues | Void ratio evolution inside shear bands in triaxial sand specimens studied by computed tomography[END_REF][START_REF] Fu | Fabric evolution within shear bands of granular materials and its relation to critical state theory[END_REF]Fu and Dafalias 2015;[START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF][START_REF] Salvatore | Determination of the critical state of granular materials with triaxial tests[END_REF]. In this study, the index M ed [START_REF] Liu | Sustainability of internal structures during shear band forming in 2D granular materials[END_REF] has been adopted to define the shear band domain. M ed refers to an absolute difference of incremental deviatoric strain inside and outside the shear band with a trial width divided by their sum. Regarded as an optimization problem, the shear band width has been obtained according to the maximum M ed within a reasonable range of the trial shear band width. More details are available in Liu, Nicot, and Zhou 2018. In the following section, labels with * refer to measures within a homogeneous domain (shear band for dense specimen and whole sample for loose specimen). Taking the kinematic evolution as its base, the next section will examine the evolution of the stress and void ratio.

Mechanical response along biaxial loading paths

As recalled in the introduction, the critical state (CS) refers to a state where stresses, void ratios and fabrics tend to be steady at relatively large deformation when shear strain further increases (Theocharis et al. 2019). Biaxial tests used as references have been simu-lated to define the critical state line which refers to a collection of critical states obtained at different confining pressures. Typical macroscopic responses in biaxial tests are illustrated in Fig. 3.4, where deviatoric stress and volumetric strain evolve with respect to axial strain ε 2 . Figures on the left hand show the results from the dense samples. In these three figures, the deviatoric stresses show an increase before a peak is reached, and a decrease followed by a steady regime with small fluctuations. The stress peaks rise gradually with the increase in confining pressure from 40 kPa to 100 kPa. The volumetric strain shows a small contractancy, at first, and transfers to dilatancy before a steady state is reached. On the other hand, for the loose specimens, q and ε v continuously rise up to steady states, as shown on the right-hand panel of Fig. 3.4. Localized and diffuse patterns have developed within the dense and loose samples, respectively.

As highlighted in Subsection (3.2.1.2), it is important to focus on the domain within the shear band when considering void ratio and fabric indexes in order to characterize material scale properties by considering only homogeneous domains. The same point holds for stresses even if the stress heterogeneity is never reported in strain localization problems. Do stress patterns exhibit significant differences inside and outside the shear band? Is the stress path within the shear band the same as in the whole sample? Usually, stresses are obtained macroscopically from the external forces applied on the sample boundaries. To obtain the stresses within shear band, we have adopted a local definition based on Love-Weber stress at the grain scale. The mean stress tensor σp ij of each particle can be computed based on the contact forces applied, thanks to the Gauss theorem (Nicot et al. 2013b;[START_REF] Liu | Strong contacts, connectivity and fabric anisotropy in granular materials: A 3D perspective[END_REF]). More details are available in aforementioned literature. The stress tensor for a given domain Ω can be computed based on all N p particles in the domain as

σ * ij = 1 Ω Np p=1 σp ij V p .
In Figure 3.5, global stress components have been compared with local stresses computed within the shear band domain for the 100kPa-dense sample under a biaxial loading path. Stresses based on boundaries, including principal stresses σ 1 and σ 2 in the whole sample, are compared with σ * 1 and σ * 2 considering the domain of the shear band based on the grain scale stress definition. In addition, the relative orientation of the stress tensor characterized by (θ σ * 1 -θ σ 1 ) is presented. It is observed that the evolution of stresses within the shear band and the whole sample is far more similar than it is for the void ratio [START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]) in 2D. Thus, it is reasonable to adopt stresses based on sample boundaries instead of stresses within the shear band when investigating critical state. This approach has been adopted throughout the thesis. A precise characterization of the loading path within the shear band and its difference from the whole specimen is an issue that remains to be examined. Some preliminary results on stress rotation within shear bands can be found in (Liu et al. 2020a) and more investigations and discussions on this issue will be

p-q-e space analysis

The CSL is characterized by mean stress p, deviatoric stress q and void ratio e, usually plotted separately in planes (p, q) and (p, e). So, here, mechanical responses (p, q, e) from proportional strain tests and biaxial tests have been compared in p -q and p -e planes. As for the dilatant proportional strain tests, Fig. 3.6 shows that the path of the deviatoric stress and the mean stress (p, q) for each case, starts from the initial state of around (100 kPa, 0), evolves to the MSRL and, at the end, changes direction to turn back toward the zero-stress level along an asymptotic line after having formed a 'loop'. Result of a mixed biaxial test/proportional strain test is also presented in Fig. 3.6. The sample was led to CS by a biaxial loading under 100kPa lateral stress upon which the loading mode was switched to a dilatant proportional strain loading. It can be seen that the stress path evolves from the very beginning along the CSRL with decreasing p and q.

p-q plane

Figure 3.5 illustrates that the stresses in the sample are relatively homogeneous at least in a 2D condition. The vanishing of the principal stresses along the dilatant proportional tests occurs simultaneously with the evolution from localized to diffuse failure (Figure 3.3). The vanishing of the principal stresses offers two possibilities: (1) liquefaction occurs only within the shear band and the behavior outside the shear band is regarded to be close to an elastic solid under unloading; (2) the kinematic pattern evolves gradually to become diffuse throughout the whole sample. In addition, it is worth noting that all computed results display an asymptotic behavior in the p -q plane, approaching the CSRL independently of the imposed dilatancy ratio, as shown in Fig. 3.6. This kind of stress track, especially the loop, has scarcely been reported in the literature for dilatant proportional strain tests in either laboratory [START_REF] Ibraim | Strain path controlled shear tests on an analogue granular material[END_REF][START_REF] Sibille | Granular plasticity, a contribution from discrete mechanics[END_REF][START_REF] Daouadji | Experimental and numerical investigation of diffuse instability in granular materials using a microstructural model under various loading paths[END_REF] or DEM simulations [START_REF] Nicot | Micro-macro analysis of granular material behavior along proportional strain paths[END_REF], with exception of the work reported by (Chu, 2 ), and relative orientation of principal stress (θ * σ 1 -θ σ 1 ). The 100kPa-dense sample under biaxial loading path is adopted.

Lo, and Lee 1992; [START_REF] Guo | Stress dilatancy and fabric dependencies on sand behavior[END_REF][START_REF] Daouadji | Phase transformation states of loose and dense granular materials under proportional strain loading[END_REF].

To further refine the classification of the stress response along proportional strain paths, a variety of proportional strain paths have been additionally simulated using specimens with different densities. The results are shown in Fig. 3.7. Based on these results, the stress response path along proportional strain tests can be categorized into four types under different densities and volumetric strain rates, as illustrated in Fig. 3.8. Group 1 results from one typical contracting strain path along which p and q increase continuously, and Group 2 comes from a typical dilatant proportional strain path leading to zero mean pressure p [START_REF] Daouadji | Experimental and numerical investigation of diffuse instability in granular materials using a microstructural model under various loading paths[END_REF]Nicot et al. 2013a). Group 3 represents stress paths from a dense sample under proportional strain tests with a relatively small dilatant rate, whereas Group 4 demonstrates a stress path under a proportional strain test with a contracting rate in a dense sample, similar in trend to the stress path under an undrained loading.

It is worth paying attention to Group 3 because of the infrequently-reported 'stress loop'. A dense sample with imposed dilatancy was experimentally investigated by [START_REF] Ibraim | Strain path controlled shear tests on an analogue granular material[END_REF][START_REF] Chu | Strain-softening behavior of granular soil in strain-path testing[END_REF][START_REF] Daouadji | Phase transformation states of loose and dense granular materials under proportional strain loading[END_REF], and numerically by [START_REF] Guo | Stress dilatancy and fabric dependencies on sand behavior[END_REF] and Nicot, Sibille, and Hicher 2015. The p -q relation from [START_REF] Ibraim | Strain path controlled shear tests on an analogue granular material[END_REF] along such a loading path has been drawn in Fig. 3.8. The stress path was stopped early with an axial strain of around 5%. If the test had been conducted further, with the growth of volume and a constant number of particles, an increasing number of contacts would probably have opened that would have resulted in a stress drop at some point. Indeed, the complete increase and decrease in q along such a loading path reported in [START_REF] Daouadji | Phase transformation states of loose and dense granular materials under proportional strain loading[END_REF] produced a stress path that was more like a back-and-forth line than a loop. This may have been the case because it was a middle dense sample with a not very pronounced softening. Both the limited shearing strain and the use of the middle dense samples could explain why no such loops were observed, even for similar loading conditions [START_REF] Nicot | Micro-macro analysis of granular material behavior along proportional strain paths[END_REF]. Such loops were, however, reported in [START_REF] Guo | Stress dilatancy and fabric dependencies on sand behavior[END_REF][START_REF] Shi | Induced fabric anisotropy of granular materials in biaxial tests along imposed strain paths[END_REF]. The Group 3 type has also been obtained in 3D simulations not included in this chapter which focuses on 2D simulations, is presented in the perspectives in the last chapter. It would be worth to conduct a detailed 3D investigation in the future.

As for the type of Group 4, stress softening was observed both along isochoric and contracting proportional strain loading paths. The slight decrease in the deviatoric stress q is linked to strain localization within the sample, since shear bands were generated in these two simulations. As underlined in [START_REF] Guo | Stress dilatancy and fabric dependencies on sand behavior[END_REF], the slope of the stress path under a contractant proportional strain loading depends on the volumetric strain rate. But, according to the results in Fig. 3.6, the critical stress ratio in proportional strain tests is independent of λ when the axial strain is large enough, as shown in Fig. 3.6.

p-e plane

Void ratios e * within the shear band are tracked for dense samples, and global void ratios e * for loose samples. The void ratio is given in 2D by e * = (A sb -A s )/A s , where A sb and A s denote, respectively, the total and solid areas of the shear band. The evolution of e * is represented in Fig. 3.9, with respect to mean stress p.

As for proportional strain tests shown in Fig. 3.9 (a), the results converge to a master curve in the p -e * plane independent of the dilatant rate characterized by λ. When the starting point of the proportional strain loading is at critical state under a biaxial loading, the (p, e * ) curve will follow the master curve from the very beginning of the proportional strain loading.

As for biaxial tests with confining pressures of 40 kPa, 60 kPa and 100 kPa, the results illustrated in Fig. 3.9 (b) are consistent with experimental observations in [START_REF] Desrues | Void ratio evolution inside shear bands in triaxial sand specimens studied by computed tomography[END_REF]DEM simulations in Zhu et al. 2016. Namely, e * from dense samples, after dilatancy, evolves to meet e * of loose samples at the critical state. For a more precise representation of the CSL in the (p, e * ) plane, two biaxial test results with confining pressures at 10 kPa and 20 kPa have been included. These critical states with proportional strain tests in p -e * plane are compared. An interesting feature stands out: the fitting curve of p -e * from proportional strain tests crawls just along the critical state curve obtained from the biaxial tests. Even though the CS from nine biaxial tests can lead to the critical state line, the master curve obtained from fitting a curve to a cluster of p -e * data in dilatant proportional strain tests, compared to the nine CS obtained from biaxial tests, induces less uncertainty.

It is worth noting that the fitting curve is based on a power function rather than a logarithmic function. It has been proved experimentally [START_REF] Verdugo | The steady state of sandy soils[END_REF][START_REF] Li | State dependant dilatancy in critical-state constitutive modelling of sand[END_REF]) that, unlike in clay, the critical state line for granular materials cannot be estimated by a straight line in the e -logp plane. Thus, the commonly used CSL expression e c = A + B lnp c is generally not representative for granular materials, where the e c and p c are the critical void ratio and the critical pressure; A and B are two constants characterizing a given material. The relation of e c and p c for granular materials can be expressed as e c = D -E ( pc F ) G , where F is a reference pressure, usually the atmospheric pressure (101 kPa) for convenience, and D, E, and G are dimensionless constants that can be identified from experiments on a given material [START_REF] Li | Modeling of dilative shear failure[END_REF][START_REF] Yin | Stress-dilatancy behavior for sand under loading and unloading conditions[END_REF].

In summary, from Fig. 3.6 and Fig. 3.9, it can be observed that the master curves (p, q) and (p, e * ) from dilatant proportional strain tests agree well with the critical state lines from biaxial tests after a sufficiently large strain level has been imposed (or immediately if the initial state is already a critical state). Hence it is shown that the classical critical state surface defined in p -q -e space for a given granular material can be obtained by performing one single dilatant proportional strain test. To the best of our knowledge, no such results have been reported in the literature. For a further analysis of this relation, the following sections will investigate the micro/mesostructures of the samples under proportional strain and biaxial tests, as well as under mixed biaxial/proportional strain tests.

Fabric-related critical state locus

Even though the critical state is most often characterized by the relation between the mean stress p, the deviatoric stress q and the void ratio e when stress and strain rate directions are fixed (Theocharis et al. 2019), the evolution of the fabric is also an important feature to be added to the p -q -e space in order to stay within the framework of Anisotropic Critical State Theory (ACST) that requires the values of the fabric to be steady.

Fabric tensor analysis

A fabric tensor quantifies microstructural orientation-related characteristics of the material in a tensorial form (Fu and Dafalias 2015). In this subsection, the inter-particle contact normal directions are characterized by the second-order fabric tensor of

F c = 1 N c Nc k=1 n k ⊗n k ,
where N c is the total number of inter-particle contacts in the assembly; n k is the unit vector representing the normal direction of the kth contact. Only contacts in the shear band are considered. The norm of the normal at contact has to be normalized by the Lode angle in order to be unique in 3D. In the present 2D case there is no Lode angle; thus, no need for normalization. To make sure that the CS investigated within shear bands satisfies the third condition required by ACST as mentioned before, the orientation of the fabric tensor has been compared with the direction of the plastic flow which can be related to the stress tensor for monotonic radial loading (Li and Dafalias 2012). The deviatoric orientation is less than 2.5 degrees throughout the whole biaxial test on the 100kPa-dense sample. all inter-particle contacts investigated are exactly the same.

The evolution of α * c has been plotted in Fig. 3.10 with respect to mean stress p. Except for the proportional strain loading starting from the CS, all the specimens have similar fabric characteristics in terms of variable α * c close to 0 at the initial states, due to the isotropic initial states where the inter-particle contact normal direction distribution is uniform in all the possible directions. As shown in Fig. 3.10 (a), all α * c -p graphs from the proportional strain tests converge towards a master curve. A linear fit is given for the proportional strain tests and repeated in Fig. 3.10 (b) to compare with the results from the biaxial tests. The relation between proportional strain tests and biaxial tests characterized by (p, q, e * ) has also been observed in the p-α * plane, namely, the master curve obtained in proportional strain tests has gathered the critical states from biaxial tests. More importantly, all specimens under biaxial and proportional strain loading paths end up with the same anisotropy in the ultimate regime. The convergent feature at ultimate states along different loading paths is in line with the findings reported in the literature [START_REF] Fu | Fabric evolution within shear bands of granular materials and its relation to critical state theory[END_REF]Fu and Dafalias 2015).

Grain loop evolution

To obtain richer information on the critical state, we have analysed grain loops related statistics in this subsection. As illustrated in Chapter 2, grain loops, enclosed by contact branches, are formed by tessellating the material area [START_REF] Satake | A discrete-mechanical approach to granular materials[END_REF]. The side number of the loop influences considerably its deformability. In this section, loops will be categorized into four groups according to the side number, referred to as l i (i ∈ [3, 4, 5, 6p]), where '6p' refers to a side number equal to or greater than 6. The percentages of categories are given by r i = n i /n t , where n i is the population of the grain loop with a side number equal to i, and n t refers to the total population of loops. It is worth noting that the reason for categorizing r 6p is that r 7p (percentage of categories with a side number equal or greater than 7) evolves in a similar way as r 6 [START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF].The evolution of r * i (i ∈ [3, 4, 5, 6p]) within the shear band for dense specimens and the whole area for loose specimens is tracked along both the dilatant proportional strain loading path and the biaxial loading path.

The results as function of the mean stress p are presented in Fig. 3.11 -3.14. Small triangles and squares highlight the start and the end points of the curves and a fitting curve is given based on proportional strain tests in each figure. Two significant characteristics manifest themselves again: (1) the measures of r * i (i ∈ [3, 4, 5, 6p]) from proportional strain tests converge to form a master curve independent of the value of λ; (2) the master curve agrees very well with the critical states characterized by the grain loop indexes obtained For the shortest loop l 3 , r * 3 continuously decreases from the initial to the ultimate point along all different loading paths. As for the most deformable type of loop l 6p , it is in a quite different situation as shown in Fig. 3.12. The value of r * 6p for the dense samples increases during loading, but it decreases for the loose samples.

Different to the fabric tensor based on the statistics of inter-particle contact normals, investigated in Subsection 3.4.1, that has been widely accepted as a quantitative measure of fabric anisotropy in granular materials [START_REF] Li | Dissipation consistent fabric tensor definition from DEM to continuum for granular media[END_REF], the scalar r * i is rather a measure characterizing the deformability of granular materials, namely the void ratio. Unlike the fabric anisotropy normalized by the specific volume [START_REF] Li | Dissipation consistent fabric tensor definition from DEM to continuum for granular media[END_REF], the evolution of r * i depends on p as illustrated in Fig. 3.11 -3.14.

Mixed proportional strain and biaxial loading paths

A complex loading path combining proportional strain and biaxial paths has been simulated to confirm that proportional strain tests can indeed determine the position of critical state in the e -p -q space as shown in Fig. 3.6 and Fig. 3.9. It was conducted as follows: first, the proportional strain path with λ = -1.2 was imposed; six groups of state data at ε 2 = 0.0022, 0.0065, 0.0109, 0.0250 and 0.0294 were selected; then, a biaxial loading path was performed at each state while the lateral stress was maintained as it was. Such mixed loading paths allow us to assess the existence of critical states according to the definition given in the introduction (continuous shearing with no change in volume or mean pressure, while keeping aligned stress anisotropy with fabric anisotropy). The stress paths along the mixed paths are shown in Fig. 3.15, while the stress and strain responses are presented in Fig. 3.16.

In Fig. 3.15, the six transition states are labelled A, B, C, D, E and F and are marked by zoomed dots. The average critical (p, q) state is marked by squares. It can be observed that all six squares from continued biaxial tests lie along the stress track of proportional strain tests. As seen in Fig. 3.16 , when the loading path was adjusted to a biaxial test at ε 2 = 0.0022, a small compaction, followed by a large dilatancy, appears before a steady regime was reached; when the biaxial test was imposed at ε 2 = 0.0065 or 0.0109, the volume increased at first before it reacheed a steady state; when ε 2 = 0.0250 or 0.0294, the following biaxial loading path fluctuated around constant values for both deviatoric stress q (kPa)

q of proportional strain loading q of continued biaxial loading Figure 3.17: The relation between imposed dilatancy from proportional strain test and unconstrained dilatancy if switched to a biaxial condition according to results in Fig. 3.15 and Fig. 3.16. Note that soil mechanics convention is adopted with positive compression.

Figure 3.18: The relation between imposed dilatancy/contractancy in proportional strain tests and unconstrained dilatancy/contractancy expected in biaxial tests in the plane of axial strain vs. volumetric strain (left) and the possible stress paths in the plane of deviatoric stress vs. mean stress (right). The major solid lines represent the imposed dilatancy/contractancy, namely the proportional strain loading paths; the fins denote the incremental unconstrained dilatancy/contractancy if the loading is switched to biaxial conditions. and volumetric strain, which means that the switching points correspond to a critical state. The stress-strain evolution of biaxial tests further proves that at relatively large deformation, the proportional strain loading path will evolve along the critical state line defined in biaxial tests.

Based on the results in Figures 3.15 and 3.16, the definition of critical state can be expressed as follows: for a sample with fixed stress and strain rate directions, a p -q -e state is considered as a critical state if the application of a biaxial loading under the same lateral pressure, starting from this p -q -e state, leads to zero volume strain and zero deviatoric stress evolution. Note that for such a state, a rotation of the principal stress direction with p and q constant will nevertheless result in a change of e. This result observed for instance in Theocharis et al. 2019 shows that the critical state cannot be described only by the three internal state variables p -q -e; additional variables related to the microstructural anisotropy due to loading direction are necessary to define a state equation for the critical state line (Li and Dafalias 2012).

Figures 3.15 and 3.16 also provide information for comparing the unconstrained dilatancy obtained in biaxial tests with the imposed dilatancy in proportional strain tests2 . As illustrated in Fig. 3.17, the solid line represents the imposed dilatancy with respect to axial strain. The fixed dilatant rate is characterized by the slope of the line. The fins give the unconstrained dilatancy rate (namely a direction of incremental volumetric strain), that would be observed if the loading were switched to a biaxial one. At the beginning, the imposed dilatancy is larger than the unconstrained dilatancy (from point O to X), the mean stress p decreases slightly from p 0 to X. When the imposed dilatancy is smaller than the unconstrained dilatancy expressed as εunconstrained v < εimposed v < 0 (from point X to Y), the stress response p increases from X to Y. The unconstrained dilatancy decreases continuously thanks to the imposed dilatancy. The states at ε 2 = 0.0065 and 0.0109 in Fig. 3.16 belong to this domain. At points X and Y, where the unconstrained dilatancy equals the imposed dilatancy, the mean stress p reaches the peak. The state at ε 2 = 0.0207 in the proportional strain test is close to point Y. After point Y, the unconstrained dilatancy decreases continuously to become smaller than the imposed dilatancy, described as

0 > εunconstrained v > εimposed v till zero unconstrained dilatancy at point Z ( εunconstrained v = 0).
At point Z, the stress path reaches the CSL defined by biaxial tests. If a biaxial path is followed after point Z, εunconstrained v = 0 occurs as shown in the results from mixed tests switched at ε 2 = 0.0250 and 0.0294. The further imposed dilatancy after point Y leads the stress path to turn back to approach point O along the CSL. In conclusion, when the imposed dilatancy is smaller (larger) than the unconstrained dilatancy, the mean stress p increases (decreases).

The results in Fig. 3.7 suggest that the relation between imposed dilatancy/contractancy in proportional strain tests and dilatancy/contractancy expected in biaxial tests can be categorized as shown in Fig. 3.18. The strain path 1 represents a situation for which the imposed dilatancy is larger than the unconstrained dilatancy (contractancy can be regarded as the negative dilatancy) from the beginning. This strain path leads to a liquefaction presented as stress path 1. Dilatant proportional strain loading path 2 leads to a more complex response, as discussed before.

The proportional strain test with zero volume change (often referred to as undrained test) can be categorized as one of two types, depending on whether εunconstrained v < 0 at the beginning (path 3) or εunconstrained v > 0 (path 4). The unconstrained dilatancy along loading path 3 tends toward zero (horizontal direction), when approaching the critical state. Conversely, loading path 4, representing an isochoric test on a loose sample, results in liquefaction.

As for the contracting proportional strain path, two groups can be identified according to whether the stress path crosses the CSL as paths 5 and 6 do. At the beginning of curve 5, the unconstrained contractancy being larger than the imposed contractancy, the mean stress p decreases slightly. Shortly afterwards, the unconstrained contractancy becomes smaller than the imposed contractancy, which leads to the increase in p. The direction of the fins tending to become horizontal is synchronized with the p -q curve approaching the CSL when strain localization takes place within the sample. The representation of curve 5 is based on results from simulations of mixed biaxial/contracting proportional strain loading paths (these results are not presented in this chapter whose focus is upon dilatant proportional strain loading paths). Curve 6 represents the situation where the unconstrained contractancy is always smaller than the imposed contractancy, which leads to a continuous increase in p.

To summarize, the mean stress will increase in the dilatancy domain, when the imposed dilatancy rate is smaller than the unconstrained dilatancy rate. When the imposed dilatancy is larger than the unconstrained dilatancy, the mean stress decreases and tends to 0 for large volumetric strain. When the imposed dilatancy is first smaller and then larger than the unconstrained dilatancy, the stress path is likely to undergo a stress loop. In the contractancy domain, when the imposed contractancy is larger than the unconstrained contractancy, the mean stress increases. On the other hand, an imposed contractancy being smaller than the unconstrained contractancy will result in a decrease in p.

These numerical findings can also be interpreted in terms of dilatancy angles if we recall that the dilatancy angle is defined in 2D by sinψ = εv / εd = εv /(2 εaxial -εv ). As long as it can be assumed that granular materials follow a standard non-associated elasto-plastic behavior (which is true in 2D or axisymmetric conditions for instance [START_REF] Nicot | Basic features of plastic strains: from micromechanics to incrementally nonlinear models[END_REF], the dilatancy angle characterize the flow rule, and the unconstrained dilatancy can be viewed as a material property. This broadens the scope of the analysis summarized in Figure 3.18 insofar, as making it independent of the loading paths considered (at least for 2D or axisymmetric conditions). The present results are also consistent with previous experimental and theoretical studies [START_REF] Daouadji | Phase transformation states of loose and dense granular materials under proportional strain loading[END_REF][START_REF] Darve | Liquefaction phenomenon of granular materials and constitutive stability[END_REF]. In [START_REF] Daouadji | Phase transformation states of loose and dense granular materials under proportional strain loading[END_REF], the dilatancy rate defined by εv / εaxial can be related to the effective mean pressure. In the criterion of liquefaction given in Darve 1996, the dilatancy angle obtained from the flow rule and the mean pressure at a given initial void ratio have a non-linear negative correlation.

Conclusion and outlook

By simulating proportional strain tests, biaxial tests and complex loading paths combining biaxial and proportional strain paths using DEM, we have explored through a series of mechanical and fabric indexes the relation between proportional strain tests and critical states from biaxial tests.

In considering the results obtained from biaxial tests, we have observed that all structural variables (e * , α * c , r * i ) within shear bands for specimens experiencing a localized kinematic pattern will converge to the same values as for specimens experiencing a diffuse kinematic pattern at critical state along the same biaxial loading path. These results are consistent with the conclusion drawn by [START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF] that localization and diffuse patterns share the same fabric properties.

As for proportional strain tests, when the dilatant rate is relatively small and the axial strain is large enough, a stress path in the p -q plane will finally reach the critical state line defined by biaxial tests after experiencing a stress loop. The curves on the p -e * (void ratio), p -α * c (fabric anisotropy intensity defined from contact normal direction) and p -r * i (population of grain loops) planes also converge towards a master curve regardless of the magnitude of dilatancy characterized by λ = -1.2, -1.3, -1.4. In combining these results, it can be inferred that any dilatant volume change will always result in similar (p, q, e * ,

Conclusion and outlook

53

α * c , r * i ) states after the material memory has been erased.

More interestingly is the relation observed between the mechanical states reached along a proportional strain test and the critical states obtained from biaxial tests within the framework of ACST. Macroscopic responses (p, q, e * ) and fabric-related indices (α * c , r * i ) in homogeneous domains along biaxial tests evolve towards the evolution curve of those variables from proportional strain tests at critical state. Given these results, only one single test is therefore necessary to construct the classical critical state line in (pqe *α * cr * i ) space for any granular material. Therefore, the CS concept can be generalized to a wide class of loading paths which shows that CS acts as a general attractor irrespective of the loading path considered. This idea of proportional strain paths enriched considerably the data base for confirmation of CST and ACST.

The mixed proportional strain and biaxial loading paths also confirm this relation. A more general definition of CS can be given as follows: a p -q -e state, for a sample with constant stress and strain rate directions, after a given loading history, is considered at critical state if the application of a biaxial loading under the same lateral pressure starting from this p -q -e state leads to zero volume strain and zero deviatoric stress evolution.

In addition, the comparison between the imposed dilatancy/contractancy along proportional strain paths and unconstrained dilatancy/contractancy expected along biaxial paths provides valuable information to the interpretation of the various loading paths obtained in p -q plane for proportional strain loading paths. Along proportional strain tests, the volumetric strain, which tends to be steady, produces the stress path approaching the CSL.

Proportional strain loading paths are attracted by the critical state line. For samples initially at critical state (after a biaxial loading for instance), a proportional strain loading imposes the evolution of the sample state along the critical state line. Even for samples not initially at critical state, when subjected to proportional strain loading, the microstructure rearranges itself and reaches geometrical arrangements corresponding to critical state. This memory fading is probably a key ingredient affecting complex systems and may also be at the origin of emerging properties. It will be of great interest to further investigate the underpinning mechanisms of memory fading according to the kinematic pattern (localization vs. diffuse mode). In addition, from a microscopic point of view, critical state is a dynamic equilibrium, where only statistics are constant. Most constitutive models for granular materials describe the critical state from a phenomenological point of view ignoring the micromechanical dynamics. Allowing critical state to emerge in micromechanical constitutive models from a physical viewpoint will provide a challenging topic for the fu-ture. One limitation of this study is that all conclusions above are drawn based on 2D simulations. The extension to 3D, numerically and experimentally, will be attempted in the future. the geomechanics community, form the focus of Critical State Theory (CST) [START_REF] Casagrande | Characteristics of cohesionless soils affecting the stability of slopes and earth fills[END_REF][START_REF] Roscoe | On the yielding of soils[END_REF][START_REF] Been | Critical state of sands[END_REF] which stands as a cornerstone of constitutive theory for granular materials. Constitutive laws built on CST are however phenomenological and cannot account for the underlying conformational transitions which are responsible for the emergence of critical states. To this respect, only multiscale models may capture the microscale detailed balance responsible for critical state. But, many aspects of these transitions remain poorly understood, and no multiscale model exists in which i) the minimal set of conformational transitions responsible for critical state and ii) the rate of the reorganisation are embedded. Indeed, novel structures, transitions and regimes are still being uncovered in granular materials under shear, mostly thanks to the use of Discrete Element modeling (DEM), as well as advanced experimental equipment [START_REF] Rocks | Hidden topological structure of flow network functionality[END_REF].

One line of advance, which is pursued in this chapter, is to adopt a dynamical systems theory (DST) approach to improve fundamental understanding of the critical state by focusing on conformational transitions from the perspective of two important classes of mesostructures, namely force chains and grain loops [START_REF] Tordesillas | Transition dynamics and magic-number-like behavior of frictional granular clusters[END_REF][START_REF] Walker | Self-assembly in a near-frictionless granular material: conformational structures and transitions in uniaxial cyclic compression of hydrogel spheres[END_REF][START_REF] Pucilowski | Rattler wedging and force chain buckling: metastable attractor dynamics of local grain rearrangements underlie globally bistable shear banding regime[END_REF]). Here we similarly employ a DST approach to uncover insights on these structures' lifespan dynamics. Our approach can help enrich multiscale constitutive models, such as the H-model (Nicot and Darve 2011b) which can incorporate the physics of grain loops and force chains to the extent that the emergent critical states can be reproduced. More broadly, a proper understanding of the birth-death dynamics of these structures benefits a gamut of applications where granular materials are subjected to shear, including: the reduction of energy consumption and greenhouse gas emissions in processes like comminution (Ben-Nun, Einav, and Tordesillas 2010; [START_REF] Walker | Complex networks in confined comminution[END_REF], discovery and design of granular materials with enhanced properties, and geohazard forecasting [START_REF] Tordesillas | Spatiotemporal slope stability analytics for failure estimation (SSSAFE): linking radar data to the fundamental dynamics of granular failure[END_REF], to name a few examples.

Prior studies have shown that the emergent mesostructures of force chains and grain loops are important building blocks for self-organization in 2D and 3D systems [START_REF] Drescher | Photoelastic verification of a mechanical model for the flow of a granular material[END_REF][START_REF] Majmudar | Jamming transition in granular systems[END_REF][START_REF] Antony | Link between single-particle properties and macroscopic properties in particulate assemblies: role of structures within structures[END_REF][START_REF] Tordesillas | Force cycles and force chains[END_REF][START_REF] Walker | Self-assembly in a near-frictionless granular material: conformational structures and transitions in uniaxial cyclic compression of hydrogel spheres[END_REF][START_REF] Pucilowski | Rattler wedging and force chain buckling: metastable attractor dynamics of local grain rearrangements underlie globally bistable shear banding regime[END_REF]. Prior work has highlighted that some statistical descriptors of the microstructure (e.g. fabric tensor, proportions of chained particles, proportions of the different categories of loops) reach constant values at critical state [START_REF] Rothenburg | Critical state and evolution of coordination number in simulated granular materials[END_REF][START_REF] Fu | Fabric evolution within shear bands of granular materials and its relation to critical state theory[END_REF][START_REF] Kruyt | On micromechanical characteristics of the critical state of two-dimensional granular materials[END_REF][START_REF] Zhu | On a common critical state in localized and diffuse failure modes[END_REF][START_REF] Wang | Multiscale modelling of granular materials in boundary value problems accounting for mesoscale mechanisms[END_REF]Deng et al. 2021c;[START_REF] Tordesillas | Force cycles and force chains[END_REF]. Indeed birth-death dynamics and related concepts like "lifespan", is a promising way to think of the critical state regime. Granular shear simulations and experiments, including those focusing on stick-slip dynamics, have shown that this regime is governed by the competing mechanisms of growth/birth versus collapse/death of force chains and these events drive evolution towards, and during, stationary states of systems under shear [START_REF] Rechenmacher | Evolution of force chains in shear bands in sands[END_REF][START_REF] Tordesillas | Transition dynamics and magic-number-like behavior of frictional granular clusters[END_REF]Kuhn 2016b;Kuhn 2016a;[START_REF] Pucilowski | Rattler wedging and force chain buckling: metastable attractor dynamics of local grain rearrangements underlie globally bistable shear banding regime[END_REF]. Experimental evidence of detailed balance in sheared dense planar granular systems is recently reported in [START_REF] Sun | Experimental evidence of detailed balance in granular systems[END_REF]. On the other hand, a DST analysis of DEM simulations for a 3D system in the critical state regime uncovered bistable dynamics [START_REF] Pucilowski | Rattler wedging and force chain buckling: metastable attractor dynamics of local grain rearrangements underlie globally bistable shear banding regime[END_REF], which is governed by imperfect detailed balance and underlying conformational transitions having connections to the rattler dynamics reported in Wautier, Bonelli, and Nicot 2019b. Despite the significant attention paid to the topic, the open question(s) that arise from all these studies remains: Why do granular systems evolve towards a stationary limiting state under continued shear? What mechanisms attract the system towards the critical state?

To quantify whether jamming dominates unjamming in evolving systems, the biologyinspired concepts of lifespan and life expectancy are particularly useful. Starting from a reference microstructure conformation, the lifespan of a mesostructure corresponds to the "time" it has continuously existed since its formation. Its life expectancy corresponds to the remaining time it will survive before disappearing. By comparing the distribution of lifespans and life expectancies, one can assess whether the birth rate of a particular class of mesostructure dominates its death rate. These concepts have been applied to granular materials for instance by [START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]Wautier, Bonelli, and Nicot 2018a;Deng et al. 2021a. Indeed, these concepts enable the quantification of the renewal rate of the microstructure in both 2D and 3D conditions under an external forcing time parameter. Specifically, it measures how fast the system is losing the memory of its past state under an external forcing. Consequently, they are particularly well adapted to treat granular materials from the perspective of dynamical systems with memory. In all these contexts, the time may correspond to physical time or any relevant forcing parameter, depending on the system of interest (for biaxial tests considered in the present study, the axial strain is the time parameter).

The chapter is organized as follows. In section 4.1, the mechanical responses of granular specimens under biaxial tests with discrete element modeling (DEM) are analyzed. This is followed in Section 4.2 by a discussion on critical state from a stationarity perspective. To this respect, the concepts of lifespan and life expectancy of mesostructures are applied to force chains and grain loops. Birth and death rates of chained particles and grain loops are analyzed in Section 4.3. The relations between microstructure dynamics and loading conditions are presented in Section 4.4. The memory fading processes of the initial conformation and along the proportional strain tests in granular materials is analyzed in Section 4.5.

Numerical setup

A set of biaxial simulations are adopted. Confining pressures of 40kPa, 100kPa and 400kPa are considered. Note that stresses are obtained by considering an out of plane arbitrary dimension equal to 0.04 m. Dense and loose samples are prepared as indicated in Table 4.1. With p 0 equal 40 kPa, 100 kPa and 400 kPa, the ratio of the average contact overlap and particle size < u n > /d 50 are 6.6 × 10 -4 , 1.6 × 10 -3 and 5.6 × 10 -3 , respectively. For the loose samples (the friction angle is set to 35 • ), the values of < u n > /d 50 are 9.2 × 10 -4 , 2.1 × 10 -3 and 7.6 × 10 -3 . More details about sample preparation can be found in the previous chapter.

As the previous chapter, the stress and strain states are described in 2D conditions, with the deviatoric stress being equal to q = σ 2 -σ 1 , and the volumetric strain being defined as ε v = ε 1 + ε 2 , where σ 1 and σ 2 are the principal stresses and ε 1 and ε 2 are the principal strains. Axes 1 and 2 refer to lateral and vertical directions, respectively.

The so-called critical state

In the literature, it is common to report fluctuations of the deviatoric stress and the volumetric strain around steady-state values rather than a smoothly constant evolution under biaxial loading. The evolution of micro and meso indices like fabric tensor and ratio of loops population also exhibits fluctuations at critical state. In this section, macroscopic deviatoric stress and volumetric strain are put in relation with the evolution of force chains and grain loops populations, respectively.

Stress-strain analysis

Macroscopic responses of the tested specimens are illustrated in Fig. 4.1, in which deviatoric stress and volumetric strain evolve with respect to the axial strain. The figure shows 4.2. The so-called critical state 59 typical stress and strain responses of dense and loose samples under the biaxial loadings.

In the dense cases, the deviatoric stress q increases at first then decreases to a stationary value, while the specimen contracts first and then dilates before reaching a steady state. With an increase in the confining pressure there is an increase in the shear strength and a decrease in the dilation. On the other hand, the deviatoric stress and volumetric strain in the loose samples evolve monotonously in the course of loading. It can be observed that, at critical state, both deviatoric stress q and volumetric strain ε v fluctuate around a constant value. Based on these mechanical responses, both the stationarity and dynamical properties of CS will be investigated in a micro-mechanical framework.

The incremental deviatoric strain distribution can be analyzed to characterize kinematic patterns, as introduced in literature (O' Sullivan, D. Bray, and Li 2003;[START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]. The maps are computed based on the interpolation of particle incremental displacements over meshed specimen [START_REF] Liu | Sustainability of internal structures during shear band forming in 2D granular materials[END_REF][START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]. When dense specimens are considered, shear bands may develop, traversing the whole specimen as shown in Fig. 4.2. It is the spatial domain within shear bands for dense cases and the whole sample for loose cases that are considered in this study [START_REF] Fu | Fabric evolution within shear bands of granular materials and its relation to critical state theory[END_REF][START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]. The spatial extension of the shear band was determined based on the method proposed in [START_REF] Liu | Sustainability of internal structures during shear band forming in 2D granular materials[END_REF].

Force chains and deviatoric stress

Force chains have attracted considerable interest owing to their significant role in force transmissions in granular systems [START_REF] Drescher | Photoelastic verification of a mechanical model for the flow of a granular material[END_REF][START_REF] Radjai | Force distributions in dense two-dimensional granular systems[END_REF]. Three conditions proposed by [START_REF] Peters | Characterization of force chains in granular material[END_REF], as recalled in Chapter 2, are adopted to detect a force chain .

In this chapter, the particles belonging to force chains are called chained particles. The proportions of chained particles is given by r cha = N cha /N tot , where N cha and N tot are the number of chained particles and the total number of particles within the investigated domain, respectively. The ratio within the shear band r * cha and the whole sample r cha are compared in Fig. 4.3.

It can be seen that the evolution with strain of both r * cha and r cha tends to reach a constant value (with fluctuations) when the deviatoric stress and the volumetric strain become constant at CS. r * cha is larger than r cha after the generation of the shear band, with .1: Evolution of deviatoric stress q and volumetric strain ε v with respect to the axial strain ε 2 along biaxial tests. Note that soil mechanics convention is adopted with positive compression and positive contraction. greater fluctuations. This may be results in a stronger rearrangement of chained particles within the shear band. It is worth noting that the ratio of chained particles fluctuating around a mean constant value highlights that important force chain reorganizations take place at the critical state.

Grain loops and volumetric strain

In a 2D granular material, meso-loops are composed of a set of contacting grains forming a closed loop. It is usually categorized according to the edge number (from L3 with three particles to L6p for loops containing six or more particles), a factor which has a significant influence on the behavior of a loop, especially its deformability. The percentage of each category reaches a steady state that corresponds to the critical state regime at the specimen scale [START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]. In this section, we link the grain loop populations (L3, L4, L5, L6p) to the volumetric evolution along the biaxial test by introducing the volumetric strain of the different grain loops.

The incremental strain tensor of each type of loop i can be expressed as follows:

d εLi = 1 |Ω| k |Ω k Li |d εk Li i ∈ (3, 4, 5, 6p) (4.1)
where Ω is the volume of the global domain considered and Ω k Li is the volume of a given loop k belonging to the category L i . The incremental strain d εk Li of each loop 'k' can be computed assuming a uniform deformation within that loop according to the strain definition proposed by [START_REF] Liu | Sustainability of internal structures during shear band forming in 2D granular materials[END_REF] The contribution of local incremental strains to the global incremental strain can be expressed by Equation 4.2:

< d ε >= 1 |Ω| Ω d εdS = d εL3 + d εL4 + d εL5 + d εL6p (4.2)
Figure 4.4 presents the volumetric strain of each category of loops and their integration during the biaxial loading as well as the evolution of the global volumetric strain.

Figure 4.4 shows a good agreement between the macroscopically computed volumetric strain and the mesoscopically averaged one. Despite reaching a zero volume change on average, it can be seen that each loop category reaches constant dilative or contractive rates at critical state. The dilatancy from L 3 , L 4 and L 5 , and the contractancy from L6p should be related with the inherent property of each loop: the grains belonging to L i (i ≤ 5) have limited space to move inward to induce contractancy, whereas L6p loops have larger void space to move inward. As a result, critical state can only be achieved from a collective process by loops transforming from one category to another.

Hidden dynamics at critical state

As highlighted in the previous sections, the jamming and unjamming behavior at CS includes dynamic mechanisms at the mesoscale [START_REF] Rechenmacher | Evolution of force chains in shear bands in sands[END_REF][START_REF] Tordesillas | Transition dynamics and magic-number-like behavior of frictional granular clusters[END_REF]Kuhn 2016b;[START_REF] Sun | Experimental evidence of detailed balance in granular systems[END_REF]. To this respect, the dynamical rate are investigated in this section through computing the lifespan and life expectancy of chained particles (i.e. particles belonging to force chains) and grain loops. The incremental evolutions of the sets of chained particles and grain loops are checked for small axial strain increments of δε 2 = 0.1%.

Generating and vanishing process of chained particles

All chained particles are tracked from birth, i.e. the first time they become part of a given force chain, to death, i.e. the first time they no longer belong to any force chain. This enables to construct probability density functions (PDF) for the lifespan and the life expectancy of chained particles at any time of the simulation. The lifespan and life expectancy PDFs are displayed in Fig. 4.5 for two arbitrary microstructure conformations (M and N) in the critical state regime. Because the sample is known to experience shear band localization (see Fig. 4.2), the analysis is restricted to particles belonging to the shear band in Fig. 4.5.

At the reference conformation M (resp. N) in Fig. 4.5 (a), the PDFs for the lifespan and life expectancy demonstrate the generating and vanishing processes of chained particles within the shear band domain. Few chained particles existing at M or N comes from the initial conformation, which proves that the initial arrangement of particles has been erased when the sample reaches CS. The symmetric shape of the lifespan and life expectancy PDFs shows that the rearrangement of chained particles reaches an equilibrium at CS, which is expected for a stationary regime. More interestingly the shape of the PDFs for the reference conformations M and N is almost the same. Force chains are always rearranging at critical state and the rate of renewal of the set of chained particles is constant. These features can be observed in 100kPa -loose with lower renewal rates from Fig. 4.5 (b).

The PDFs have been averaged to get rid of statistical noise coming from the finite size of the chained particle population at different strain levels belonging to the critical regime (ε 2 = 0.0610, 0.0708, 0.0806, 0.0904, 0.1002 and 0.1100 for the 100kPa-dense sample and ε 2 = 0.0605, 0.0708, 0.0802, 0.0906, 0.1005 and 0.1106 for the 100kPa-loose sample). Since force chains have a limited life duration at critical state, the six microstructure conformations can be considered as statistically independent (not correlated), which allows us to build a single PDF from the six data sets. The corresponding lifespan and life expectancy PDFs are presented in Figure 4.6 with respect to incremental evolution of the axial strain (∆ε 2 ). Thus, ∆ε 2 at current state is 0. ∆ε 2 > 0 and ∆ε 2 < 0 correspond to the lifespan and life expectancy of chained particles, respectively.

The renewal rates in dense and loose samples appear to be different at first glance, which would contradict the existence of a unique critical state, whatever the initial sample density. To find the origin of the apparent faster renewal rate within the shear band in the dense sample, it is meaningful to compare the incremental strain magnitude dε m = dε 2 1 + dε 2 2 in the shear band domain for the dense sample, to the incremental strain magnitude in the whole sample domain for the loose sample. Similar to the incremental strain of each type of loops computed in Subsection 3.3, the incremental strain within the shear band can be computed by averaging the incremental strains of loops within the shear band over the area of the shear band. When averaged over the whole sample domain dε m is the same in dense and loose samples at critical state, as shown in Fig. 4.7. However, because of the strain localization in the dense sample, a larger dε * m develops within the shear band. The ratio between < dε * m > and < dε m > at CS inside the shear band is around 2.8. In other words, an axial strain increment of 1% from the boundaries of the sample imposes a strain The averaged PDFs of chained particles with respect to the incremental strain evolution (∆ε 2 ), over six microstructure conformations at different strain levels belonging to the critical regime (ε 2 = 0.0610, 0.0708, 0.0806, 0.0904, 0.1002 and 0.1100 for the 100kPa-dense sample and ε 2 = 0.0605, 0.0708, 0.0802, 0.0906, 0.1005 and 0.1106 for the 100kPa-loose sample). increment of 2.8% in the shear band domain because of strain localization. This means that the actual forcing parameter (the "time") inside the shear band is 2.8 time faster than for the whole sample.

As a result, the axial strain applied at sample scale is not the right forcing parameter in the dense case as the strain is 2.8 times more intense in the shear band domain. In order to account for this observation, we need to stretch the curves of 100kPa-dense* by the ratio 2.8 to make the lifespan and life expectancy PDF comparable between the dense and loose sample cases. Then, it can be observed in Fig. 4.8 that the PDF in dense and loose samples are coinciding. Similar results can be observed for different confining pressures, as shown in Appendix A. It can be concluded that the strain magnitude (observed in homogeneous domain of the samples) is the major factor affecting the rate of rearrangement of chained particles. We can also conclude that the uniqueness of the critical state (namely whatever the initial density, statistical descriptors converge towards the same value in biaxial tests of similar confining pressure) that holds for statistical descriptors, also holds for the underlying dynamical processes.

Generating and vanishing process of grain loops

To further understand the critical state dynamics, the generating and vanishing of grain loops are analyzed for loops composed of 3, 4, 5 and 6 or more particles. Lifespan and life expectancy PDFs for grain loops are computed. Similar observations as reported in the previous subsections are visible in Fig. 4.9 -4.11: (1) PDFs show symmetrical shapes at CS which suggests birth and death rates of grain loops are equal at CS; (2) PDFs within the shear band from 100kPa-dense coincide with those obtained in the loose case after being stretched by the appropriate ratio of 2.8. These results extends the conclusion drawn by [START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF], that the whole loose sample and the inner part of the shear band of the dense sample share the same statistical distribution of mesostructures once the critical state regime is reached.

Based on , a novel way for understanding CS from the perspective of a mesoscale dynamics can be inferred: the critical state results from the generating and vanishing of meso-clusters until a dynamic equilibrium is reached under continuous shearing. Even though dynamic rearrangements permanently occur, at any time, the material admits the same fabric characterized by the fabric tensor, the distribution of loop orders etc.. The particle rearrangements at critical state can be quantitatively characterized by the lifespan and life expectancy of force chains and grain loops. These observations still hold when other confining pressures (40 kPa and 400 kPa) are considered as shown in Appendix A.

Microstructure reorganization dynamics under different confining pressures

In this section, the influence of the confining pressure on the dynamics in the three loose samples is analyzed. As shown in Figs. 4.12 and 4.13, the renewal rate decreases with the increase in the confining pressure. The increase in confining pressure helps the mesostructures live longer (with respect to the forcing parameter of the system, i.e. the incremental axial strain). This is consistent with the fact that larger confining pressure leads to larger elastic energy storage and larger overlap at contacts, which delays contact opening and contact sliding. The averaged PDFs of grain loops with respect to the incremental strain evolution (∆ε 2 ), over six microstructure conformations at different strain levels belonging to the critical regime (ε 2 = 0.0610, 0.0708, 0.0806, 0.0904, 0.1002 and 0.1100 for the 100kPa-dense sample and ε 2 = 0.0605, 0.0708, 0.0802, 0.0906, 0.1005 and 0.1106 for the 100kPa-loose sample). .13: Lifespan and life expectancy PDFs for grain loops with respect to the incremental strain evolution (∆ε 2 ) under different confining pressures with samples 40kPa-loose, 100kPa-loose and 400kPa-loose considered.

Microstructure reorganization dynamics under different confining pressures 75

In order to quantify the rate of microstructure reorganization with respect to the axial strain, it is interesting to propose a fit for the PDFs presented so far. As a first simple guess, an exponential fit with one characteristic strain may been proposed as P DF (∆ε 2 ) = a εc e -|∆ε 2 |/ε c . ε c accounts for the typical size of the axial strain increment ∆ε 2 needed to renew the whole mesostructures. The smaller ε c , the faster the microstructure renewal under the external forcing ∆ε 2 . However, such a fit cannot capture well the shape of the whole PDFs. Depending on ε c , the fit is able to account either for the beginning or for the tail of the PDFs. This observation motivated the use of a double exponential fit in the following form:

P DF (∆ε 2 ) = P c 1 ε c1 e - |∆ε 2 | ε c1 + (1 -P c ) 1 ε c2 e - |∆ε 2 | ε c2 (4.3)
where ε c1 and ε c2 (ε c1 < ε c2 ) correspond to two characteristic strains accounting for the dynamic evolution of the system. ε c1 and ε c2 account for the typical amplitude of the axial strain increments ∆ε 2 needed to renew the mesostructures on the short term and long term, respectively. The characteristic probability P c ∈ [0, 1] can show the relative contribution of the short and long term mechanisms. When P c < 0.5, it is the long term mechanism that dominates; the short term mechanism will be more important when P c > 0.5; the two mechanisms are evenly matched when P c = 0.5. Figs. 4.14 and 4.15 show that double exponential fit is able to account for the shape of all the lifespan and life expectancy PDFs. The values of ε c1 and ε c2 for different confining pressures are presented in Table 4.2 and shown in Fig. 4.16 with respect to the critical mean stress.

As shown in Table 4.2 and Fig. 4.16, large loops (L5 and L6p) live shorter (faster reorganization) than small loops (L3 and L4) and force chains. It is consistent with the assumption that loops connected to force chains open prior to force chain bending (Wautier, Bonelli, and Nicot 2018b).

When it comes to the effect of confining pressures, the renewal rate of meso-structures which can be characterized by 1 ε c1 and 1 ε c2 decreases with the increase in the confining pressure. This effect differs according to the type of meso-structures. L3 and L4 seem to be more sensitive to the confining pressures, especially in the long term mechanism, which could be due to the fact that these small, less compressible loops can store larger (stabilizing) elastic energy before breaking. On the other hand, large loops (L5 and L6p) are more compressible, then less sensitive to confining pressure, as the increase in p does not have as much stabilization effect as for smaller loops. Figure 4.16: Two characteristic strains (ε c1 and ε c2 ), the characteristic probability (P c ) and the ratio ε c1 ε c1 with respect to the mean stress at the critical state. The critical mean stress under different confining pressures (40 kPa, 100 kPa and 400 kPa) are 59 kPa, 143 kPa and 559 kPa, respectively.

Microstructure reorganization dynamics under different confining pressures 79

As for chained particles, even though force chains can share particles with large loops, their stability (and thus lifespan) is strongly dependent on the existence of small order loops around them. Thus, the effect of confining pressure on force chains lifespan stem from the different responses of loops under the confining pressure according to their category (small vs. large ones), which may explain why the slope of the ε c1 and ε c2 for chained particles is smaller than the one for small loops but larger than the one for large loops.

In addition, the long term mechanism increasingly has dominate effect, since P c , always being less than 0.5, decreases with the growth of critical pressures.

The form of the proposed fit may appear arbitrary at first glance; but a double exponential PDF is the signature of the existence of two mechanisms responsible for microstructure reorganization acting over two different "time" scales. A first mechanism provokes microstructure reorganizations over small axial strain scales of typical magnitude ε c1 , while a second mechanism induces microstructure reorganizations over much larger axial strain scales of typical magnitude ε c2 .

The existence of two reorganization mechanisms is indeed consistent with the visualization of the incremental deviatoric strain map in Fig. 4.2 for the loose sample. The map is far from being homogeneous over the whole sample domain, and zones with large increments of deviatoric strain are visible (in the form of red zones in Fig. 4.2.). Contrary to shear bands, these zones are not persisting when the axial strain increases. Indeed, these zones may be seen as shear transformation zones (STZ) [START_REF] Schall | Structural rearrangements that govern flow in colloidal glasses[END_REF] and characterized by the recent concept of shear chain [START_REF] Darve | Slip lines versus shear bands: two competing localization modes[END_REF].

As a result, we conjecture that ε c1 relates to the axial strain needed for a shear transition zone to develop, while ε c2 > ε c1 corresponds to the axial strain needed to develop shear transition zones everywhere in the sample domain. This conjecture can be initially verified by the color maps of lifespan and life expectancy of loops in the sample 100kPa-loose, as shown in Figs. 4.17 and 4.18, respectively.

The short and long threshold are set to 4δε 2 and 30δε 2 , respectively. In Fig. 4.17 (a), red loops are newly-generated at the current conformation. The loops that have lived longer than 4δε 2 are filled in grey. Fig. 4.17 (a) shows that the microstructure is experiencing a series of localized reorganizations (visible as clusters of loops of the same colors) affecting a limited proportion of the sample domain. In the long term, these series of reorganizations

Microstructure dynamics under evolving conditions: memory effects in granular materials 81

will affect the whole sample domain as shown in Fig. 4.17 (b), where loops that have been existing not less than 30δε 2 are very limited. In Fig. 4.18 (a), grey loops are those that will live longer than four steps. It can also be observed in Fig. 4.18 (b) that beyond the long threshold, only a small part of the loops will still be alive.

Microstructure dynamics under evolving conditions: memory effects in granular materials

The critical state has been paid great attention at both macroscopic and microscopic scales with the emphasis of the stationary section. However, how the CS emerges from the initial state under shearing from a micromechanical perspective has been highlighted much less often. Tordesillas 2007 pictured the whole fabric evolution process through introducing the jamming and unjamming concepts in granular materials under shearing. Shear band generating and expanding processes are characterized by microstructure dynamics in 2D simulations. It has been widely accepted that the critical state is independent on the initial void ratio. Matsushima and Blumenfeld 2017 found that the microstructural characteristics are governed mainly by the packing procedure, and the effects of intergranular friction and initial states are details that can be scaled away. But the fading process of the initial memory of the sample is unclear. It has been proved that the sample under the proportional strain test can evolve along the critical state line normally defined by a series of biaxial tests in Chapter 3. But how the memory of the structure has been erased needs further exploration.

It is the external loading and internal rearrangements that work together to form the stationary state. As commented in Matsushima and Blumenfeld 2017, the self-organization leaves a fingerprint that manifests in the emergence of intriguing universality-like behavior. The structural characteristics are not completely random and bear the hallmark of the self-organization process that brought the medium to mechanical equilibrium. Thus, the purpose of this chapter is to explore the competition between microstructure dynamics and loading condition, through memory effects characterized by the life duration of grain clusters. We focus on grain loops and force chains.

Particularly, the life duration of the mesostructures introduced in the previous chapter is recalled in this section. Three tests are considered: biaxial tests in a dense and a loose sample, and a dilatant proportional strain test in a dense sample. In Fig. 4.19 (a), it can be observed that a part of the chained particles at the initial state is transferred to be non-chained particles at an early state (before stress stabilizes) during the loading, while the other remains continuously in force chains until the end of the test. The permanent particles will be proved to be located outside of the shear band domain later in this section.

Fading process of the initial memory

In general, a chained particle being converted to be non-chained can be caused by (1) the force chain being transferred to a normal chain without collapse, buckling or bending because of redistribution of contact forces, or (2) collapse, buckling or bending of the force chain. It is important to distinguish the two modes. The former one does not relate to particle rearrangements, while the latter does.

In Fig. 4.19 (a), the decreasing process of the PDF can be divided into three stages according to the curve slope, namely the rate of transition from chained particles to nonchained particles: a rapid drop at stage I, a gentle decrease at stage II and nearly no change at stage III. The transition point between the first two slopes is close to the characteristic point, which corresponds to the equality between incremental elastic strain energy and incremental plastic strain energy [START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF], and the appearance of shear band (Liu et al. 2020a). The first drop should be mostly related to the stress anisotropy resulting from the redistribution of contact forces. A number of vertical force chains are built in the assembly. The second decrease corresponds to the softening process resulting from strain localization, during which the buckling and bending of force chains occurring within the shear band, erasing the memory at the initial state. The stationary stage corresponds to the CS characterized by constant deviatoric stress and volumetric strain. It demonstrates that once the sample reaches the CS, chained particles from initial state nearly stop disappearing (the remaining ones belong indeed to the sample domain out of the shear band). As discussed in the previous chapter, few chained particles existing at the critical state come from the initial state, which demonstrates that the initial arrangement of particles has been erased when the sample reaches CS within the shear band. Thus, the chained particles living throughout the whole loading processes almost all belong to the outside shear band domain. The decrease in the PDF, namely the vanishing of chained particles at the initial state, is related to the building up of the stress anisotropy. The assembly builds a vertical force network. This network will not live long as the lateral dilatation of the sample imposed by the biaxial loading increasingly open the lateral contacts supporting the force chains. At the characteristic point, plasticity mechanisms become dominant, force chain bend and collapse, which results in the gentle decreasing part of the PDF.

By contrast, in Fig. 4.19 (b) few chained particles survive from the initial state remaining to the end of test in a loose sample. This is consistent with the diffuse failure mode that grain reorganization occurs everywhere. Comparing the beginning stage of the two PDF curves, the first drop in the dense sample is stronger and quicker than that in the loose sample. This is because the strain localization drives an intense and fast reorganization within the band (the strain rate is 3.2 times larger as analyzed in subsection 4.3.1).

Grain loops

To further understand the critical state dynamics, the life expectancy of grain loops at the initial states is analyzed in this subsection. The disappearance of a grain loop composing of a cluster of grains is confirmed by checking whether some grains come in or out of it in the next simulation time step. Thus, the life duration of grain loops can be a direct evidence of structural self-organization. The PDF for life expectancy of grain loops at initial state is presented in Fig. 4.20.

Similar interesting observations that have been reported in the previous Subsection can be found in Fig. 4.20: (1) in the dense sample, around 80% (accumulation of PDF) of the loops existing at the initial state vanish gradually before the end of the test, and the rest remains until the end, located outside the shear band; while the loose sample totally reorganize. (2) The vanishing process is much smoother and gentler in the loose sample than in the dense sample. In Fig. 4.20 (a), the PDF curve for Loop3 is nonlinear before approaching a steady value around 0. The strongest vanishing point of L3 locates close to the characteristic point of the volumetric strain. Change from contractancy to dilatancy corresponding to the breakage of L3. The non-linearity is becoming gentler with the increase in the order of loops. As for L6p, the largest vanishing occurs at the initial point. In the dense sample, L6p is the most sensitive type of loops at the contractive stage. In the loose sample, PDFs for all types of loops are in similar trends. A Slight difference among them can found in the initial values. L6p has be largest initial value compared to the others. In both dense and loose samples, l6p is the weakest category at the beginning stage where the deviatoric stress level is low.

To demonstrate the memory fading process intuitively, the life expectancy characterized by incremental axial strains has been shown in the color map in Fig. 4.21. Figure 4.21 (a) and (c) presents the loops destroyed within 10 increments δε 2 of axial strain (δε 2 = 0.1%) for the dense and loose samples, respectively. A band can be observed in Fig. 4.21 (b). Yellow loops corresponds to loops alive until the end of the test. They are located outside the band, coexisting with loops having short life expectancy. Particle rearrangements occur early and intensely within the shear band, and over a longer time scale outside the shear band. The life expectancy of mesostructures is a promising index that can be used to identify the shear band. The vanishing of loops in the loose sample is relatively homogeneous in space, as shown in Fig. 4.21 (d).

Memory fading process along proportional strain test

It has been proved that one single dilatant proportional strain test can present a critical state line normally defined by a series of biaxial tests. Lifespan dynamics of cluster conformations would help to refine the mesostructure reorganization along the decrease in mean stress during a proportional dilatant strain loading, thus, to better understand the critical state line as an attractor.

In Fig. 4.22, life expectancy PDFs (probability density functions) of different types of grain loops are presented under the dilatant proportional strain loading considering the whole sample, with deviatoric stress and volumetric strain curves. The decrease in each PDF is rapid initially and reach at a stationary state close to zero. There is little grain loops remaining until the end of the test.

The lifespan and life expectancy PDFs at the deviatoric stress descending stage are explored within the shear band. The two axial strains, 0.0610 and 0.1094, are considered as reference states, as shown in Fig. 4.23. The generating and vanishing processes for grains loops appear to be symmetrical, which can be further observed in Fig. 4.25 presenting the absolute value of PDFs with respect to the axial strain evolution (∆ε 2 ). In addition, lifespan and life expectancy PDFs are explored with the whole sample considered, as shown in PDFs. The gap could be the origin of the leakage as a imperfect balance. The leakage is believed to mainly happen around the blurred boundary of the shear band area since the lifespan and life expectancy PDFs within the shear band domain are symmetric. The life expectancy curve is above the lifespan curve which may due to the decrease in the mean stress resulting in slower structural reorganization.

A series of life expectancy color maps are shown in Fig. 4.27, with the initial state (a), the peak of deviatoric stress (b), states belonging to the convergent stage (c-f) being considered as marked in Fig. 4.24.

In the first four maps, the shear band can be observed. Within the shear band, the color become darker and darker from (a) to (d) corresponding to faster and faster rate of grain rearrangement, which will be further discussed. Outside the shear band, the domain is similarly darker and darker. However, the shape of the band has little change from (a) to (d). The band is no more evident at the state (e) and (f) where the vanishing of grain loops inside and outside the shear band is almost at the same rate. These maps demonstrate that the process from localization to liquefaction along the dilatant proportional strain test does not result from the expanding of the shear band to the whole sample, but is caused by the unloading process with decrease of stresses and the disappearing of the shear band.

Conclusion

The critical state normally emphasized as a stationary state in geomechanical fields was revisited and characterized based on a series of analysis on the rearrangement of particles from a dynamic perspective at the mesoscale. We have shown that the critical state results from the balance competing between generating and vanishing of meso-clusters. We have established that force chains and grain loops have a regular, short life duration (0.4-3% compared to the whole 15% loading) at critical state, corresponding to a rapid fading of the material memory along its loading history. The critical state not just "forgot" the initial microstructure conformation but its recent past in the stationary regime in a sample under continuous shearing. We have shown that the uniqueness of the critical state that holds for statistical descriptors also holds for the underlying dynamical processes.

By fitting the data, we have shown that force chains and loops of different sizes have different lifespan and life expectancy which increases with the mean stress level. Larger loops have a shorter life due to their deformability and force chain life duration is similar to the life duration of small order loops which is consistent with the known relationship between force chains and their small order supporting loops. Moreover, the microstructure reorganization has been proved to rely on two mechanisms that act over very different axial strain increments. A local mechanism related to sheared mesostructures provoke a rapid (i.e.for small axial strain increments) reorganization of the microstructure while the characteristic strain related to the nucleation of new sheared mesostructures account for the microstructure reorganization over larger strain increments. The recent concept of shear chain [START_REF] Darve | Slip lines versus shear bands: two competing localization modes[END_REF] will probably help precise the elastic and plastic mechanisms responsible for the observed dynamics of critical state. Complementary investigations deserve being conducted to analyze the velocity of strain waves in the sample, since the dynamic load transfer characteristics can be also interpreted within the context of wave propagation theory [START_REF] Sadd | DEM simulation of wave propagation in granular materials[END_REF]. Such a future work may help seeing whether the unjammed microstructure keeps a partial memory of its previous jammed state.

The new findings presented in this chapter might help to find clues to explicitly incorporate CS features in micro-mechanical constitutive models, as well as to limit the number of different loading paths needed to construct the data base in data driven models because of limited hysteresis effect. In relation with one particular micromechanical model, it would be possible to consider the rearrangement of the hexagonal loop meso-structures on the H-model (Nicot and Darve 2011b;Wautier et al. 2021) to achieve a better description of the critical state, based on the renewal rate of mesostructures as quantified in this study.

If the present study deals with critical state dynamics in granular material, the concept, tools and methodology are rather generic and may be applied to a great variety of other evolving complex systems to reveal their hidden dynamics.

Chapter 5

Critical state and the H-model Based on the extensive investigations on the critical state conducted in the previous chapters, it was shown that mesostructural transformations under external loading is a central ingredient that should be included in multiscale constitutive modeling. In this chapter, the framework of the H-model is reviewed, and inspected at mesoscale on a biaxial test in Section 5.1. In Section 5.2, a method is proposed to introduce dynamical reorganization between generating and vanishing of meso-clusters through a deactivation/reactivation procedure. This process is inspired by results obtained in Chapter 4. Based on the updated framework of H-model, the global stress-strain in a biaxial test is simulated, and the corresponding meso-structural analysis is presented. Additionally, predictions under two different densities and two different confining pressures are compared separately to assess the performance of the updated H-model with respect to critical state prediction.

Strain localization hypothesis

The global incremental strain δε ij is localized to each hexagon according to the equations as follows:

δl 1 = -l 1 δε ij n i n j δl 2 = -l 2 δε ij t i t j (5.1)
These formula corresponds to the relative length variation in direction n and t. Note that no distortion is taken into account in the model for the sake of preserving the symmetry of the H-cell.

Directional behavior

Each H-cell representing one direction can be as illustrated in Fig. 5.2. The state of a H-cell can be defined by three geometrical parameters: the opening angle α, the inter-granular distances d 1 and d 2 . These parameters relate to the cell dimensions l 1 and l 2 as

l 1 = d 2 + 2d 1 cosα l 2 = 2d 1 sinα, (5.2)
The mechanical behavior of each hexagon depends on both the contact law between adjoining grains and the mechanical equilibrium with external forces. An elastic-frictional contact law is adopted. Such a contact law relies on three parameters: a normal stiffness k n , a tangential stiffness k t , and an inter-granular friction angle ϕ g . The normal and tangential contact forces can be expressed in an incremental form as

δN i = k n δu i n δT i = min[|T i + k t u i t |, tanϕ g (N i + δN i )]ξ -T i , (5.3)
where ξ is the sign of the quantity T i + k t u i t , u i n and u i t are the normal and tangential relative displacements of the grains at contact i. Given the symmetries in the H cell, we only need to consider the contact forces called N 1 , N 2 and T 1 as shown in Fig. 5.2, The corresponding incremental displacements are obtained by differentiating the expression of the branch vector joining the grains in contact:

δu 1 n = -δd 1 δu 2 n = -δd 2 δu 1 t = -δd 1 δα, (5.4)
Mechanical balance of grain 2 reads, along directions n:

N 2 = N 1 cosα + T 1 sinα + G 2 = N 1 cosα + T 1 (sinα + i G 2 ) (5.5)
where i G 2 = 0 when G 2 = 0 in the original version of H model (Nicot and Darve 2011b), i G 2 = 1 when G 2 = T 1 in the updated version (Wautier et al. 2021). Equation 5.5 differentiates into

cosαδN 1 -(i G 2 + sinα)δT 1 -δN 2 + (N 1 sinα -T 1 cosα)δα = 0 (5.6)
Combining these equations at the H-cell scale, the incremental evolution δl 1 and δl 2 can be related to δd 1 , δd 2 and δα in matrix form as

A   δd 1 δd 2 δα   =   δl 1 δl 2 λ   , A =   2cosα 1 -2d 1 sinα 2sinα 0 2d 1 cosα A 31 -1 A 33   (5.7)
where

A 31 = cosα if no sliding occurs cosα + ξtanϕ g (i G 2 + sinα) if sliding occurs (5.8) A 33 = [N 1 sinα -T 1 cosα -k t d 1 (i G 2 + sinα)]/k n if no sliding occurs (N 1 sinα -T 1 cosα)/k n if sliding occurs (5.9) λ = 0 if no sliding occurs ξ(1 + sinα)(tanϕ g N 1 -T 1 ) if sliding occurs (5.10)
Based on the updated variables δd 1 , δd 2 and δα, the contact evolution (δN 1 , δN 2 and δT 1 ) can be computed according to Equation 5.7.

Stress averaging

The macroscopic stress tensor σ at a material point is calculated by a directional averaging of all contact forces acting in a collection of unit H-cells (belonging to a virtual REV corresponding to this material point). This is done in two steps by first defining a meso stress σ meso for each H-cell and then averaging all the meso-stresses of the H-cell collection.

Here and below, the meso-stress tensor is defined at the scale of the unit H-cell using the Love-Weber formula:

σ meso = 1 V meso c∈Vmeso f c l c (5.11)
where f c is the contact force at contact c, l c the branch vector joining the centers of the two grains in contact and V meso the volume of the unit H-cell. Note that Bagi correction term accounting for boundary contact forces to the H-cell is not included in the computation here [START_REF] Bonelli | On the definition of an average strain tensor for twodimensional granular material assemblies[END_REF]Nicot et al. 2013b;[START_REF] Yan | Definition and symmetry of averaged stress tensor in granular media and its 3D DEM inspection under static and dynamic conditions[END_REF]. This is part of the model assumption.

Expressed in the local frame (n, t), equation 5.11 gives with soil mechanics conventions (positive compression).

     V meso σ meso nn = 4N 1 d 1 cos 2 α + 4T 1 d 1 cosαsinα + 2N 2 d 2 V meso σ meso tt = 4N 1 d 1 sin 2 α -4T 1 d 1 cosαsinα
V meso σ meso nt = V meso σ meso tn = 0

(5.12)

Finally, the macroscopic stress tensor for the REV can be determined as an integration of all meso-stresses of individual unit H-cells, knowing their statistical distribution via the probability density function ω(θ). Thus, for the entire collection of unit H-cells,

σ = 1 V π θ=0 ω(θ)V meso (θ)σ meso (θ)dθ (5.13) such that V = π θ=0 ω(θ)V meso (θ)dθ and π θ=0 ω(θ)dθ = 1 (5.14)
If σ meso is expressed as a matrix σ n,t meso in the local frame n, t, the stress matrix σ in the global frame e 1 , e 2 is readily obtained by conversion from local to global configuration, i.e. σ e 1 ,e 2 = 1 V π θ=0 ω(θ)V meso (θ)P -1 σ n,t meso (θ)P dθ (5.15)

where P is the transformation matrix, P = cosθ sinθ -sinθ cosθ (5.16)

The probability density function ω(θ) describes the statistical distribution of unit H-cells and accounts for possible anisotropy. It is assumed that the initial distribution of unit H-cells can be approximated for instance as

ω(θ) = 1 π [1 + α ω cos2(θ -β ω )] (5.17)
where α ω is a parameter ranging from 0 to 1 describing anisotropy of the hexagon cell distribution, and β ω its major principal direction. In practice, the distribution is discretized in n θ directions.

The global volume in Equ.5.15 reads as V = π θ=0 ω(θ)V meso (θ). Note that the choice of V meso is not unique as long as the strain localization equation 5.1 is updated accordingly to ensure that δVmeso Vmeso = T r(δσ). There are three main options for the volume of H-cell as shown in Fig. 5.4. From the smallest to the largest, they are labelled as V 1 meso , V 2 meso and Figure 5.4: Three options of the meso-volume: hexagonal domain (blue), rectangular domain going through the grain centers (green), bounding box domain (red)

V 3 meso . With these different meso-domains being considered, there are some changes in the porosity and the strain localization equation 5.1. Note that the strain localization equation 5.1 is consistent with the second definition of V meso (V meso = V 2 meso ).

• To consistent with Equation 5.1, it is the rectangular domain (green in Fig. 5.4) should be considered. The volume reads:

V 2 meso = l 1 l 2 = 2d 1 sinα(d 2 + 2d 1 cosα) (5.18)

• When the bounding box domain (red in Fig. 5.4) is considered, the volume reads:

V 3 meso = (l 1 + d 0 )(l 2 + d 0 ) (5.19)
In this case, Equ.5.1 should be updated with l 1 + d 0 and l 2 + d 0 instead of l 1 and l 2 . The geometric compatibility (Equation 5.2) matrix A (Equ.5.7) hold.

• When the hexagonal domain (blue in Fig. 5.4) is considered, the meso-volume V h meso reads,

V 1 meso = 2d 1 sinα(d 2 + d 1 cosα) (5.20)
A virtual rectangle with the same area of the hexagon can be considered to make the same volume adopted in both strain localization and stress averaging processes. A rectangle of l * 1 × l * 2 can be proposed. With this hypothesis, Equ.5.2 should be replaced by:

l * 1 = d 2 + d 1 cosα l * 2 = 2d 1 sinα
(5.21)

The matrix A needs to be updated as A * with two elements replaced: A * 11 = cosα and A * 13 = -d 1 sinα.

In the original version of the H-model, V meso is considered as the area of the hexagon, even though this volume is not consistent with the volume l 1 l 2 for strain localization in Equation 5.1. The option is inspired by the fact that granular materials in 2D can be tessellated by contact branches.

The general scheme of the standard 2D H-model can be summarized as follows:

1. Based on given (δε 1 , δε 2 ), update the H-cell dimensions l 1 , l 2 of each H-cell according to a kinematic localization hypothesis.

2. Update d 1 , d 2 , α through Equation 5.7.

3. Update N 1 , N 2 , T 2 based on the incremental evolution of the H-cell geometry.

4. A meso-stress is defined based on the application of Love-Weber formula for each H-cell.

5. The macroscopic stress is eventually obtained by statistical averaging of all the mesostresses.

In conclusion, the 2D H-model relies on three contact parameters (k n , k t and ϕ g ), on one geometrical parameter (the initial opening angle α 0 ), and on α ω and β ω for tuning the microstructure anisotropy.

Pathological cases

Pathological types of H-cell may be obtained owing to contact loss or gain during external loading. There are two modes of abnormal H-cell due to contact open as shown in Fig. 5.5. The meso-stress on the broken H-cell is 0. The H-cell can be over-compressed when α ≤ 30 • or α ≥ 90 • , as shown in Fig. 5.6. In this case, additional contacts are generated. The overcompressed H-cells only orientate along 90 • and for a large level of axial strain. in Python (See in Appendix B). Under this type of loading, a constant lateral stress condition is imposed together with a constant strain rate in the axial direction: δε 1 = const and δσ 2 = 0. An iterative calculation process is run to find incremental strains (δε 1 , δε 2 ) that gives the target confining pressure. The parameters [START_REF] Veylon | Modélisation numérique du mécanisme de liquéfaction des sols: application aux ouvrages hydrauliques[END_REF]) are shown in Table 5.1. Three confining pressures 100 kPa, 200 kPa and 400 kPa are considered. The number of directions θ i is set to 360.

Biaxial test at the material point scale

Stress and strain curves are presented in Fig. 5.7. It can be observed that, with the increase in p 0 , there are larger maximum stress ratios and stronger dilatancy. The trend is consistent to the results reported in Veylon 2017. If we pursue the loading up to a large level of axial strain (ε 11 was 40%), the behavior is obviously not physically relevant. As shown in Fig. 5.8, the deviatoric stress decreases to a very low value and the volumetric strain turns to contractancy at a large level of axial strain. This is due to the fact that a large number hexagons are broken and do not participate to stress transmission. Details are presented in the following section.

Mesoscale inspection of the H-model during biaxial loading

The mesostructure of granular materials is always changing both in time and space. In the H-model, directional statistics are introduced to get the average behavior at macroscale. In each direction θ i , H-cells have a specific geometry. Before introducing any additional ingredient to the H-model, it is important to clarify the limiting aspects of the model. The biaxial loading with p 0 = 200 kPa is considered in this subsection. The evolutions of opening angle α, the contact branches d 1 and d 2 and the meso-stress V meso σ meso ij are discussed as follows, for different H-cell directions θ i .

The contact branches d 1 and d 2

The lengths d 1 and d 2 can be used to detect contact loss as shown in Fig. 5.5. The directional distribution of the ratio d 1 /d 0 between the length and its initial value is shown in Fig. 5.9 (a). Five specific strain states are considered as marked in Figure 5.8. The evolution of d 1 /d 0 along the axial strain ε 11 is presented in Fig. 5.9 (b). Eight directions ranging from 0 • to 90 • are considered. V meso meso 22 1e7 0.0°1 5.0°3 0.0°3 5.0°4 0.0°4 5.0°6 0.0°9 0.0°F igure 5.14: Evolution of V meso σ meso 22 : (a) directional analysis at different axial strains ε 11 = 0.0, 0.005, 0.032, 0.128 and 0.347 labelled as O, A, B, C and D, respectively, (b) macro strain evolution along specific directions.

From Fig. 5.9 (a), we can observe that, from O to D (see Fig. 5.8), an increasing range of directions has d 1 /d 0 larger than 1, which means that more and more hexagons are broken because of contact loss at c 1 (see Fig. 5.2). At state D, there are only a small fraction of directions around 40 • and 140 • with the contact c 1 being unbroken. In Fig. 5.9 (b), the orientations close to 0 • experience compression for c 1 at the beginning of the loading; The decrease in d 1 /d 0 is smaller and smaller along the orientation θ ranging from 0.0 • to 90 • in the beginning of the loading (ε < 5%); the hexagon oriented to 90 • is the earliest one to be broken. When θ is around 40 • , d 1 has a complex non-monotonic behavior; contacts can re-exist after being open.

As for d 2 , the curves in Fig. 5.10 have some similarities and some differences to those in Fig. 5.9. Similarly, with the increase in axial strain, an increasing number H-cells have d 2 /d 0 larger than 1, mainly open contacts; at the early state of the loading, the length of d 2 decreases (corresponding to compressed contacts) in the direction 0 • , and increases (corresponding to extended contacts) in the direction 90 • . One important difference is that the value of d 2 /d 0 oriented along 0 • is always smaller than 1, namely the corresponding contacts always exist during the loading.

More and more H-cells are broken around 0 • and 90 • because d 1 or d 2 extend which results in a decrease in stress level along the loading.

The opening angle α

The value of α is closely related to the deformabitity of the mesostructure and to fabric anisotropy. The H-cell has the potential to dilate and contract with the evolution of α as shown in Fig. 5.11.

The evolution of α is shown in Fig. 5.12. It can be observed there are two directions (around 50 • and 130 • ) along which α has little change over the loading. Between the two directions, α decreases continuously corresponding to the evolution of 1 ○ in Fig. 5.11; while beyond this domain it increases as the trend 2 ○. The pathological case occurs resulting from over-compression with α < 30 • in the direction 90 • .

The deformed processes imply that the force chain composed by any three nearby grains always bends along the loading direction. This is different to what has been observed in DEM simulations where force chains tend to lie along the direction of loading after particle rearrangements. During the rearrangements, grain loops not only are deformed, but also are decomposed to form new loops. The generation and vanishing activities of loops is certainly a missing puzzle piece in the H-model to account for large deformation.

5.1.3.3

The meso-stress V meso σ meso ij The distribution of meso-stresses V meso σ meso 11 and V meso σ meso 22 in the global frame (e 1 , e 2 ) is presented in Fig. 5.13 -5.14. Only unbroken hexagons are considered to contribute to the global stresses, thus, the results are highly dependent on the length of d 1 and d 2 . It can be observed that the hexagons nearby the vertical direction govern the contribution to global stresses at the early regime. At a large level of axial strain, only a small range of directions around 40 • are still active, which is also observed in other revised versions of H-model, for example the model including a shear mechanism [START_REF] Zhu | Including a meso-structure in multi-scale modeling of granular soils[END_REF] and the 3D H-model [START_REF] Xiong | Multiscale modeling of granular materials in application to geotechnical engineering problems[END_REF].

Comments

The H-model is capable of considering the deformation not only from the relative displacement at the contact scale between grains, but also from the particle rearrangement through deforming elementary hexagonal patterns of adjoining particles. The memory of history is stored in the structure of the hexagons. Thus, the model is capable of making hardening emerge from the deformed H-cells, and the following possible softening through the breakage of H-cells, which is consistent with the micromechanical mechanisms of hardening and softening observed from DEM simulations [START_REF] Zhu | Meso-structure organization in two-dimensional granular materials along biaxial loading path[END_REF]. Theoretically, the model also has the potential to capture dilatancy and contractancy through the changes in the opening angle. Making such macroscopic characteristics emerge without introducing any empirical law at the macroscopic scale is the most significant feature of the H-model, as a breakthrough with respect to phenomenological constitutive models.

However, the model cannot describe the critical state in a satisfying way because of the continuous breakage of H-cells without creation of new H-cells. As investigated in the previous chapter, the life duration of a loop is relatively short, and the generation and vanishing process take place all along the loading though. Microstructure reorganization in granular materials not only results in loop deformation but also induces loop transformations. Only taking into account H-cell deformation is necessary but insufficient to capture the complex reorganization process in granular materials. Even though the hexagons are allowed to be broken and reformed later on in the model, there is little newly-generated hexagon after the initial state.

Emergence of critical state in the H-model

Method

It has been shown, in the previous chapters, that the dynamical competition between generating and vanishing of meso-clusters always exists even at the critical state. Thus, the H-model might be improved through introducing such a dynamical process. In order to achieve the goal, some questions should be answered first.

• How to define the death and birth of mesostructures in the H-model? This processes can be implemented through deactivating and reactivating H-cells. When a H-cell is deactivated, its contacts are set to an open configuration through releasing the overlaps at contacts c 1 and c 2 . Thus, the H-cell is reset with d 1 = d 2 = d 0 . The deactivated H-cell does not contribute to the global stresses at this loading step. During the reset, the volume of the H-cell is kept constant by adjusting the opening angle accordingly. In real granular materials, the contact network rearranges instantaneously resulting in different meso-structures with constant volume.

In granular materials, a sample under loading evolves and adapts through grain rearrangements. The fabric of the assembly evolves to be coincident with the loading direction [START_REF] Li | Micro-macro quantification of the internal structure of granular materials[END_REF]. Based on this fact, the reactivated H-cells could be adjusted from the old configuration. In the H-model, the fabric is mainly controlled by the opening angles α at mesoscale and the distribution function ω(θ) at the macroscale. To work at mesoscale, a scheme for resetting α is proposed for reactivated hexagons with the porosity being kept, as shown in Fig. 5.15. In the resetting process, the shape of the hexagon is stretched along the compressed direction.

• Which hexagons should be deactivated/reactivated? As a first step, we propose to conduct the deactivation/reactivation scheme on the H-cells that are chosen randomly at each loading step. The percentage p h of hexagons renewed over an axial strain increment δε 11 can be set according to the renewal rate ξ characterized by the axial strain. Thus, p h = δε 1 ξ . As obtained in the previous chapter, the renewal rate of L6p at critical state is ranging from 0.3% to 2.79% as shown in Section 4.4.

Figure 5.15: Schematic model to reset α with the porosity being kept. The reset 1 is triggered when the opening angle smaller than 60 • tends to decrease further, and the reset 2 when the opening angle larger than 60 • tends to increase further.

Random reset has been used in complex chaotic system in physics. It has been proven by Chen 1996 in chaotic systems, that if their velocities are randomly reset at a regular time interval ∆t according to a prescribed probability distribution, then it might happen that all the particles will eventually converge to some unique yet erratic orbit in the physical space, depending on the system style. Similar convergent effect should be expected on the H-model after randomly resetting the hexagons.

Preliminary results

A biaxial test is predicted using the updated H-model. The flowchart of the computation incorporating the deactivation/reactivation process is shown in Fig. 5.16. One confining pressure p 0 = 100 kPa is considered. The renewal rate ξ is set to 1.5%. Parameters are shown in Table 5.1. The simulated deviatoric stress (q = σ 11 -σ 22 ) and volumetric strain (ε v = ε 11 + ε 22 ) are presented in Fig. 5.17. It can be observed that the stress-strain curves tend to reach a stationary state with fluctuations. Note that it is the volume V 2 meso considered in both strain localization and stress averaging processes, as noted in green in Fig. 5.4.

The directional distribution of geometrical variables d 1 /d 0 , d 2 /d 0 and α of the H-cell is analyzed at four different strain levels ε 11 = 0.000, 0.005, 0.150 and 0.350 labelled as It is believed that this is the reason why a stationary state emerge can be obtained at the macroscopic scale in the prediction from the deactivation and reactivation scheme.

An additional simulation with a smaller initial opening angle (40 • ) has been conducted as shown in Fig. 5.21. It can be observed in Fig. 5.21 that the decrease in the opening angle results in a higher q peak, a larger softening, and stronger dilatancy. This is consistent with the fact that a denser sample has a larger initial stiffness and stronger dilatancy. When the initial angle decreases from 45 • to 40 • , the porosity decreases as shown in Fig. 5.11.

Conclusions

The H-model is capable of considering the deformation not only from contact displacements but also from particle rearrangements through deformed H-cells. It is able to make the macromechanical characteristic of granular materials emerge, e.g.dilatancy/contractancy and hardening/softening, without introducing any empirical law at macroscopic scale.

Following the spirit of directional distribution in the H-model, a method is proposed to make the critical state emerge through a deactivation/reactivation procedure of hexagons based on the dynamical view of the critical state in Chapter 4. It has been shown that the deviatoric stress and volumetric strain along a biaxial test reach a steady state. Directional analysis has been conducted. It is demonstrated that the deactivation/reactivation scheme can help to restrict the continuously increasing number of opening contacts. The updated H-model seems to have the ability to generate similar structural configurations leading to a stationary state at the macroscopic scale.

Further efforts need to be devoted to polish and precise the deactivation/reactivation scheme through analyzing the sensitivity of the parameters. Especially for the mesostructural reorganization rate ξ, the value is inconstant before reaching the critical state according to the results in the previous chapter. Detailed comparisons like anisotropy between DEM simulations and the modeling in this chapter can help to further evaluate the model. In addition, prediction of experimental data using the H-model is worth putting forward in the future.

Chapter 6

General conclusion and perspectives

Conclusion

This PhD work was devoted to micromechanical investigations and multiscale modeling of granular materials with an emphasise on the critical state concept.

Dilatant proportional strain loading and relationships with critical state

A relation has been established between the mechanical states reached along a proportional strain test and the critical states obtained from biaxial tests within the framework of ACST. Macroscopic responses and fabric-related measures in homogeneous domains along biaxial tests converge towards the evolution curve of those variables from proportional strain tests at critical state. Given these results, only one single test is therefore necessary to construct the classical critical state line in (pqe *α * cr * i ) space for any granular material. Therefore, the CS concept can readily be obtained for a wide class of loading paths. Indeed, CS acts as a general attractor irrespective of the loading path considered. The mixed proportional strain and biaxial loading paths further confirm this relation. A more general definition of CS can be given as follows: a p -q -e state, for a sample under constant stress and strain rate directions, after a given loading history, is considered as a critical state if the application of a biaxial loading under the same lateral pressure starting from this p -q -e state leads to zero volume strain and zero deviatoric stress evolution. In addition, the comparison between the imposed dilatancy/contractancy along proportional strain paths and unconstrained dilatancy/contractancy expected along biaxial paths provides valuable information to interpret the various loading paths obtained in p -q plane for proportional strain loading paths.

Dynamical view of critical state

We have shown that the critical state results from imperfect detailed balance competing between generating and vanishing of meso-clusters. We have established that force chains and grain loops have a regular, short life duration (0.4-3% compared to the whole 15% loading in a biaxial test) at critical state, corresponding to a rapid fading of the material memory along its loading history. Critical state "forgot" not just the initial microstructure conformation but its recent past in the stationary regime in a sample under continuous shearing. We have shown that the uniqueness of the critical state that holds for statistical descriptors also holds for the underlying dynamical processes. By fitting the data, we have shown that force chains and loops of different sizes have different lifespan and life expectancy which increase with the mean stress level. Larger loops have a shorter life due to their deformability and force chain life duration is similar to the life duration of small order loops. This is consistent with the known relationship between force chains and their small order supporting loops. Moreover, the microstructure reorganization has been proved to rely on two mechanisms that act over very different axial strain increments. A local mechanism related to sheared mesostructures provokes a rapid (i.e. for small axial strain increments) reorganization of the microstructure while the characteristic strain related to the nucleation of new sheared mesostructured accounts for the microstructure reorganization over larger axial strain increments.

Enriched H-model with emerging critical state

Following the spirit of directional distribution in the H-model, a method is proposed to approach the microstructural mechanisms taking place when approaching the critical state, by introducing a deactivation/reactivation procedure on H-cells, based on DEM results about the dynamical view of critical state. It has been shown that the deviatoric stress and volumetric strain along a biaxial test reach a steady state. It is demonstrated that the deactivation/reactivation scheme can help to limit contact loss and occurrence of other pathological cases in the deformation of H-cells. The updated H-model seems to have the ability to generate similar structural configurations leading to a stationary regime at the macroscopic scale in biaxial tests which corresponds to critical state.

Perspectives

Micromechanical investigations

• Bursts of kinetic energy and global reorganization Studies proved that outbursts of kinetic energies can be seen as signatures of mechanical instabilities and inertial transitions in granular materials [START_REF] Gaume | Quasistatic to inertial transition in granular materials and the role of fluctuations[END_REF]Wautier, Bonelli, and Nicot 2018b;[START_REF] Clerc | Meso-scale signatures of inertial transitions in granular materials[END_REF]. While evolving at critical state, these bursts were shown to have a finite size [START_REF] Clerc | Meso-scale signatures of inertial transitions in granular materials[END_REF], and inertial transitions are limited in space and time. Following the results on the two mechanisms forcing microstructure reorganization in this thesis, it will be interesting to give a close look to outbursts of kinetic energy, and to characterize the reorganization rate within the domain of the bursts using the lifespan PDF. So that we can further refine how successive local inertial events can result in a global microstructure reorganization. Investigations on the evolution of grain loops and force chains at the critical state and during bursts of kinetic energy can help to improve the understanding of the co-evolution of the structures and stresses.

• Further clarify the mechanism of critical state through studying transition of conformations

In this thesis, using PDFs of lifespan and life expectancy of mesostructures, we characterized the dynamical reorganization rates in granular materials during loading. The vanished mesostructures are transferred to other configuration. At local scale, we can track the move of each particle and explore the transitions of conformations and the possibility of reforming disappeared mesostructures. This could help to understand the leakage or escapes of the imperfect detailed balance (as mentioned in Chapter 4) as well as memory fading process during different loading conditions in granular materials, thus further clarify the mechanism of the critical state as an attractor.

• Further verify the conclusions obtained in the thesis in 3D situation

It will be interesting to further verify the conclusions in this thesis in 3D. Some preliminary results in 3D about dilatant proportional strain loading and their relationship to critical state are shown in Fig. 6.1 and 6.2. Three triaxial tests were simulated to obtain standard critical states. Three dilatant proportional tests within different dilatant rates characterized by κ are considered with dε v = (1 + 2κ)dε 2 . A "stress loop" can be obtained when κ = -0.7 in Fig. 6.1.

A complex loading mixing a proportional dilatant strain test and a triaxial test with the lateral stress kept unchanged has been simulated in 3D and the corresponding stressstrain responses are presented in Fig. 6.2. The transition points are marked by dots. It can be observed that ε v and q along the continued triaxial test tend to be stable. This result suggests the possibility of extending the generalized definition of critical state to 3D conditions. More efforts are worth being devoted into micromechanical investigations in the 3D simulations on force chains, and grain loops benefiting from the concept of grain loops that was extended to 3D conditions recently (Chueire et al. n.d.).

The concept of lifespan and life expectancy presented in the thesis should be valid to any types of meso-structure well defined in both 2D and 3D conditions. However, the amplitude of the reorganization rate of mesostructures in 3D are expected to be different.

• Refine the results using angular particles

In this thesis, spherical particles was used. Simulations with angular particles merit attention in future studies. The real shape of particles can be simulated using Level-Set (LS-DEM) shape description through avatar paradigm and X-ray computed tomographic (XRCT) images [START_REF] Kawamoto | All you need is shape: predicting shear banding in sand with LS-DEM[END_REF][START_REF] Duriez | A Level Set-Discrete Element Method in YADE for numerical, micro-scale, geomechanics with refined grain shapes[END_REF]. In a system with angular particles, the reorganization rate would be different since the angularity has effect on the rolling at contacts. We expect that the mesostructures will live longer for angular particles.

H-model

• Further improve the deactivation/reactivation scheme.

The deactivation/reactivation scheme can be further improved from different perspectives. The volume of H-cell was kept constant during the reset operation through adjusting the opening angle α, namely V adjust meso = V before meso . There are other possible options to reset the H-cell, for example, by controlling V adjust meso -V before meso V before meso = δε 11 +δε try 22 . A full implementation of the life expectancy PDF in the selection of the H-cells to reset may also improve the procedure. The relation of the material porosity, the loading rate and the confining pressure etc. to the structural reorganization rate ξ can be further refined. Another potential direction to improve the scheme is to find the limit condition of a loop in a sample from DEM simulations first, and then control the deactivation and reactivation of each H-cell according to the configuration itself.

• Introduce the deactivation/reactivation scheme into the 3D H-model. The H-model has been extended to 3D with mesostructures of ten grains consisting of two imbricated hexagons in perpendicular directions [START_REF] Xiong | Multiscale modeling of granular materials in application to geotechnical engineering problems[END_REF]. It would be interesting to introduce the deactivation/reactivation scheme into the 3D H-model.

Multiscale modeling of landslides

The H-model has been implemented in different computational frameworks to simulate a variety of geotechnical problems where multi-physical couplings can occur [START_REF] Xiong | A multiscale work-analysis approach for geotechnical structures[END_REF][START_REF] Xiong | A novel multi-scale large deformation approach for modelling of granular collapse[END_REF]Wautier et al. 2021). The H-model has been coupled with Finite Element Method (FEM), based on which, a strip footing problem in plane strain conditions has been analyzed as a non-homogeneous boundary value problem [START_REF] Xiong | A multiscale work-analysis approach for geotechnical structures[END_REF]. The H-model has been implemented into a Finite Difference Method (FDM) solver to perform a stress analysis for an earth dam as the upstream water level gradually increases (Wautier et al. 2021). The model has been used with smoothed particle hydrodynamics (SPH) method to simulate collapse of granular material [START_REF] Xiong | A novel multi-scale large deformation approach for modelling of granular collapse[END_REF]). One collapse is simulated based on the H-model and SPH method, as presented in Figure 6.3. The multiscale approach using the H-model is helpful in improving the understanding of non-homogeneous boundary value problems and should be an effective computational tool for the analysis of geotechnical boundary value problems.

There will be great interest to apply the improved version of H-model developed in this thesis coupling with a numerical platform, for example, SPH method to landslide modeling. The Bas-Verger landslide in the Pays d'Auge plateau (Normandy, France) could be an appropriate case to start with (Fig. 6.4). The Bas-Verger landslide is a shallow planar slide. The landslide surface morphology is attainable from previous studies based on field mapping and aerial image interpretation [START_REF] Fressard | Multi-method characterisation of an active landslide: Case study in the Pays d'Auge plateau (Normandy, France)[END_REF]. It is promising to carry out simulations along cross-sections to analyze the potential triggering processes of the landslide. Comparisons between results from multiscale modeling and from a more classical approach used for landslide hazard analysis can be helpful to further demonstrate the interest of the new approach as well as the pros and cons of each method.

  . . . . . . . . . . . . . . . . . . . . . . . 2.4 Elasto-frictional contact law: components in the contact (left), behavior in the tangential direction (right). . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Force chain definition (Peters et al. 2005). . . . . . . . . . . . . . . . . . . 2.6 Grain loops categorized by edge length in 2D granular assembly. . . . . . . 3.1 Quasi-2D DEM specimens for (a) proportional strain tests where loading parameters are ε 1 and ε 2 and (b) biaxial tests where loading parameters are ε 2 and σ 2 . The direction 2 is always the vertical loading control direction at the REV scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Evolution of the deviatoric stress q along dilatant proportional strain tests with different magnitudes of dilatancy characterized by λ = -1.2, -1.3, -1.4. Point I, point II and point III refer to onset of a well marked shear band at ε 2 = 0.010, slightly blurred shear band at ε 2 = 0.109, and the pressure p close to zero at ε 2 = 0.142, respectively. The corresponding kinematic patterns of these three points are shown in Fig.3.3 . . . . . . . . . . . . . . 3.3 Incremental deviatoric strain maps (dε d ) estimated for axial strain increments of 0.11%. Three axial strains are considered along the proportional strain path λ = -1.2: (I) onset of a well marked shear band at ε 2 = 0.010, (II) slightly blurred shear band at ε 2 = 0.109, (III) just before liquefaction (p is close to 0) at ε 2 = 0.142. The corresponding points in q -ε 2 plane are shown in Fig.3.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Deviatoric stress and volumetric strain of biaxial tests in dense (left panel) and loose (right panel) samples. Three biaxial loading paths with σ 0 = 40 kPa, 60 kPa and 100 kPa are considered. Note that soil mechanics convention is adopted with positive compression. . . . . . . . . . . . . . . .

3. 8

 8 Four categories of stress path along proportional strain tests according to DEM simulation results shown in Fig.3.7. . . . . . . . . . . . . . . . . . . . 3.9 p-e * evolution in proportional strain (a) and biaxial (b) tests. The start and end points are marked by triangles and squares, respectively. A fit based on a power function is shown for proportional strain tests and repeated to compare with the results from biaxial tests, as all dilatant proportional strain tests converge towards a master curve. The equation of the fitted curve is e * = 0.2571 -0.0227( p 100 ) 0.6274 . . . . . . . . . . . . . . . . . . . . . 3.10 p -α * c evolution for proportional strain (a) and biaxial (b) tests. The start and end points are marked by triangles and squares, respectively. A linear fit is shown for the dilatant proportional strain tests and repeated on the biaxial tests, as all the proportional strain tests converge towards a master curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 p -r * 3 for proportional strain (a) and biaxial (b) test. The start and end points are marked by triangles and squares, respectively. An exponential fit is shown for the dilatant proportional strain tests and repeated on the biaxial tests, as all the proportional strain tests converge towards a master curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 p -r * 6p for proportional strain (a) and biaxial (b) tests. The start and end points are marked by triangles and squares, respectively. An exponential fit is shown for the dilatant proportional strain tests and repeated on the biaxial tests, as all the proportional strain tests converge towards a master curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List of Figures xi 3.13 p -r * 4 for proportional strain (a) and biaxial (b) test. The start and end points are marked by triangles and squares, respectively. An exponential fit is shown for the dilatant proportional strain tests and repeated on the biaxial tests, as the all proportional strain tests converge towards a master curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 p -r * 5 for proportional strain (a) and biaxial (b) test. The start and end points are marked by triangles and squares, respectively. An exponential fit is shown for the dilatant proportional strain tests and repeated on the biaxial tests, as all the proportional strain tests converge towards a master curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 p -q evolution along combining proportional strain and biaxial loading paths. The proportional strain loading with λ = -1.2 are conducted at first until a certain axial strain level labeled as A, B, C, D, E and F, after that a biaxial loading is performed by keeping the corresponding lateral stress unchanged. The transition states and the average critical states are marked by dots and squares, respectively. The corresponding stress-strain responses are shown in Fig.3.16 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  deviatoric strain (IDS) maps at the critical states, estimated for axial strain increments of 0.1%. Two samples are considered: (a) 100kPadense sample (ε 2 = 6.1%) and (b) 100kPa-loose sample (ε 2 = 8.02%). . . . 61 4.3 Proportions of chained particles.Two samples are considered: (a) 100kPadense sample and (b) 100kPa-loose sample. Note that r * cha represents the value when the spatial domain of shear band are considered, and r cha when the whole sample is considered. The deviatoric stresses are recalled. . . . . 62 4.4 Volumetric strain of different types of loops and the global assembly during the biaxial tests in 100kPa-dense sample (a) and 100kPa-loose sample (b). 64

  Color maps of the lifespan of grains loops for the reference axial strain ε 2 = 8.02%. Two thresholds are considered corresponding to the two characteristic strains: (a) 0.4% and (b) 3%. . . . . . . . . . . . . . . . . . . . . 4.18 Color maps of the life expectancy of grains loops for the reference axial strain ε 2 = 8.02%. Two thresholds are considered corresponding to the two characteristic strains: (a) 0.4% and (b) 3%. . . . . . . . . . . . . . . . . . . 4.19 Life expectancy PDFs (probability density functions) of chained particles at the initial state of the whole sample in the dense (a) and loose (b) cases. The deviatoric stress curve (dot-dashed line) and the volumetric strain (dashed line) curves are recalled. . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv List of Figures 4.20 Life expectancy PDFs (probability density functions) of grain loops at the initial state of the whole sample in the dense (a) and loose (b) cases. The deviatoric stress curve (dot-dashed line) and the volumetric strain (dashed line) curves are recalled in the subfigure of Loop 4. . . . . . . . . . . . . . 85 4.21 Color map of life expectancy of grain loops at the initial state under the biaxial test in the dense (a,b) and loose (c,d) samples. Two thresholds are used in the color scale: 10 (a,c) and 120 (b,c) times the reference incremental axial strain δε 1 = 0.1% which corresponds respectively to 1 % and 10 % of axial strain from the initial state in the biaxial loading . . . . . . . . . . . 86 4.22 Life expectancy PDFs (probability density functions) of grain loops under the dilatant proportional strain loading considering the whole sample. . . . 88 4.23 Lifespan and life expectancy PDFs (probability density functions) of grain loops at the convergent state under the dilatant proportional strain loading considering a domain within the shear band. . . . . . . . . . . . . . . . . . 89 4.24 The six states chosen for color maps along the proportional strain tests . . 90 4.25 The absolute value of PDFs for Lifespan and life expectancy of grain loops under the dilatant proportional strain loading within the domain of the shear band. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.26 The absolute value of PDFs for Lifespan and life expectancy of grain loops under the dilatant proportional strain loading with the whole sample considered. A fit is given to the PDFs using the double exponential in Equation 4.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.27 Color map of life expectancy of grain loops at different conformations: (a) the initial state ε 2 = 0, (b) the peak of deviatoric stress ε 2 = 0.0185, (c-f) belonging to the convergent stage under the dilatant proportional strain test in the dense sample ε 2 = 0.0610, 0.0810, 0.0910, 0.1094. . . . . . . . . . . . 92 5.1 The directional distributions and the coordinates at the REV scale (e 1 , e 2 ) and the mesoscale (n, t). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.2 The symmetrical description of geometry (a) and forces (b) at the elementary hexagonal pattern of adjoining particles (H-cell). There are two types of contact c 1 and c c , the corresponding contact branches are d 1 and d 2 . . . . . 96 5.3 General homogenization scheme relating both incremental stress and strain tensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 List of Figures xv 5.4 Three options of the meso-volume: hexagonal domain (blue), rectangular domain going through the grain centers (green), bounding box domain (red) 5.5 Two modes of pathological H-cell due to contact loss. . . . . . . . . . . . . 5.6 Two modes of over-compressed H-cell with the corresponding opening angle. 5.7 The stress ratio and volumetric strain curves based on the parameters from Veylon 2017. Three confining pressures are considered. . . . . . . . . . . . 5.8 The deviatoric stress and volumetric strain along the biaxial loading with p 0 = 100 kPa until a relative large axial strain. . . . . . . . . . . . . . . . . 5.9 Evolution of the ratio between the length d 1 and its initial value d 0 : (a) directional analysis at different axial strains ε 11 = 0.0, 0.005, 0.032, 0.128 and 0.347 labelled as O, A, B, C and D, respectively, (b) macro strain evolution along specific directions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Evolution of the ratio between the length d 2 and its initial value d 0 : (a) directional analysis at different axial strains ε 11 = 0.0, 0.005, 0.032, 0.128 and 0.347 labelled as O, A, B, C and D, respectively, (b) macro strain evolution along specific directions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.11 Porosity ϕ 1 meso of the hexagon with respect to the opening angle α; sketch of the two possible deforming processes. . . . . . . . . . . . . . . . . . . . . 5.12 Evolution of the opening angle α: (a) directional analysis at different axial strains ε 11 = 0.0, 0.005, 0.032, 0.128 and 0.347 labelled as O, A, B, C and D, respectively, (b) macro strain evolution along specific directions. . . . . . . 5.13 Evolution of V meso σ meso 11 : (a) directional analysis at different axial strains ε 11 = 0.0, 0.005, 0.032, 0.128 and 0.347 labelled as O, A, B, C and D, respectively, (b) macro strain evolution along specific directions. . . . . . . . . . 5.14 Evolution of V meso σ meso 22 : (a) directional analysis at different axial strains ε 11 = 0.0, 0.005, 0.032, 0.128 and 0.347 labelled as O, A, B, C and D, respectively, (b) macro strain evolution along specific directions. . . . . . . . . . 5.15 Schematic model to reset α with the porosity being kept. The reset 1 is triggered when the opening angle smaller than 60 • tends to decrease further, and the reset 2 when the opening angle larger than 60 • tends to increase further. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.16 Flowchart of the computation of biaxial loading incorporating the deactivation/reactivation scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Deviatoric stress and volumetric strain with respect to the axial strain using the updated H-model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.18 Polar diagrams of the ratio between the length d 1 and its initial value d 0 : at different axial strains ε 11 = 0.000, 0.005, 0.150 and 0.350 labelled as O, A, B and C, respectively (see Fig.5.17). . . . . . . . . . . . . . . . . . . . . . 118 5.19 Polar diagrams of the ratio between the length d 2 and its initial value d 0 : at different axial strains ε 11 = 0.000, 0.005, 0.150 and 0.350 labelled as O, A, B and C, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.20 Polar diagrams of the opening angle α: at different axial strains ε 11 = 0.000, 0.005, 0.150 and 0.350 labelled as O, A, B and C, respectively (see Fig. 5.17).120 5.21 Deviatoric stress and volumetric strain with different initial opening angles 40 • and 45 • being considered, predicted by the updated H-model. . . . . . 121 6.1 Stress paths in p -q plane along triaxial tests and proportional strain tests with different dilatant rates characterized by κ (dε v = (1 + 2κ)dε 2 ). . . . . 126 6.2 Stress-strain responses along mixed proportional strain and triaxial loading paths in 3D. At the transition state, the lateral stress is kept unchanged. The loading of 100kPa-dense (κ = -0.7) as shown in Fig.6.1 is adopted here. εv = (1 + 2κ) ε1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 6.3 simulation of the collapse of granular column using SPH particles and the H-model: (a) initial configuration and (b) color map of the deformation in direction 1 after collapse. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 6.4 Aerial view of the Bas-Verger landslid in 2006 with associated localization of the field measurements conducted on the landslide (Fressard et al. 2016). 131 A.1 PDFs for chained particles under two confining pressures: (a) 40kPa, PDFs from the 40Kpa-dense are stretched by strain magnitude ratio 1.9; (b) 400kPa, the 400kPa-dense stretched by 3.2. . . . . . . . . . . . . . . . . . . 134 A.2 PDFs for grain loops under two confining pressures: (a) 40kPa, PDFs from the 40Kpa-dense are stretched by strain magnitude ratio 1.9; (b) 400kPa, the 400kPa-dense stretched by 3.2. . . . . . . . . . . . . . . . . . . . . . . 135 List of Tables 3.1 Initial void ratios e 0 of prepared samples. The sample '100kPa-dense' is used in proportional strain tests) . . . . . . . . . . . . . . . . . . . . . . . 4.1 Initial void ratios e 0 of prepared samples. . . . . . . . . . . . . . . . . . . . 4.2 The fitting parameters for different confining pressures. . . . . . . . . . . . 5.1 Parameters used in biaxial tests Veylon 2017 . . . . . . . . . . . . . . . . . and objectives . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Structure of presentation . . . . . . . . . . . . . . . . . . . . . . . . 3

  of presentationThis thesis summarizes three years of research that led to two journal publications (JMPS and PRE) and three conference proceedings (Alert materials workshop 2019, 2021, and Powder and Grains 2021) at the meeting point of geomechanics and statistical physics.

Figure 2

 2 Figure 2.1: Typical deviatoric stress and volumetric strain responses in a dense (left) and a loose (right) granular assembly under a biaxial loading path.

Figure 2

 2 Figure2.2: images of a system of photoelastic disks for an isotropically compressed state (left) and a sheared state (right)[START_REF] Luding | Contact force measurements and stress-induced anisotropy in granular materials[END_REF]

Figure 2 . 4 :

 24 Figure 2.4: Elasto-frictional contact law: components in the contact (left), behavior in the tangential direction (right).

Figure 2 . 5 :

 25 Figure 2.5: Force chain definition (Peters et al. 2005).

Figure 2 . 6 :

 26 Figure 2.6: Grain loops categorized by edge length in 2D granular assembly.

  simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Mechanical response . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2.1 Mechanical response along dilatant proportional strain loading paths 27 3.2.2 Mechanical response along biaxial loading paths . . . . . . . . . . . . 28 3.3 p-q-e space analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.3.1 p-q plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.3.2 p-e plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4 Fabric-related critical state locus . . . . . . . . . . . . . . . . . . . 39 3.4.1 Fabric tensor analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.4.2 Grain loop evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.5 Mixed proportional strain and biaxial loading paths . . . . . . . 46 3.6 Conclusion and outlook . . . . . . . . . . . . . . . . . . . . . . . . . 52

  Figure 3.1: Quasi-2D DEM specimens for (a) proportional strain tests where loading parameters are ε 1 and ε 2 and (b) biaxial tests where loading parameters are ε 2 and σ 2 . The direction 2 is always the vertical loading control direction at the REV scale.

  Figure 3.2: Evolution of the deviatoric stress q along dilatant proportional strain tests with different magnitudes of dilatancy characterized by λ = -1.2, -1.3, -1.4. Point I, point II and point III refer to onset of a well marked shear band at ε 2 = 0.010, slightly blurred shear band at ε 2 = 0.109, and the pressure p close to zero at ε 2 = 0.142, respectively. The corresponding kinematic patterns of these three points are shown in Fig.3.3

Figure 3 .

 3 Figure 3.6 illustrates the stress paths of the biaxial tests and the proportional strain tests. The three straight lines result from biaxial tests with dense samples. The initial states are marked with triangles and ultimate points are highlighted by squares labelled as A, B and C. The critical stress ratio line (CSRL) is given as a dashed line based on the three

Figure 3 . 4 :Figure 3

 343 Figure 3.4: Deviatoric stress and volumetric strain of biaxial tests in dense (left panel) and loose (right panel) samples. Three biaxial loading paths with σ 0 = 40 kPa, 60 kPa and 100 kPa are considered. Note that soil mechanics convention is adopted with positive compression.

Figure 3

 3 Figure3.6: Stress paths of biaxial tests and proportional strain tests in p -q plane. The start and end points are marked by triangles and squares, respectively. The critical stress ratio curve and the maximum stress ratio curve are drawn according to the critical states and the maximum value from biaxial tests, respectively.

Figure 3 . 8 :Figure 3

 383 Figure 3.8: Four categories of stress path along proportional strain tests according to DEM simulation results shown in Fig.3.7.

Figure 3

 3 Figure 3.10: p-α * c evolution for proportional strain (a) and biaxial (b) tests. The start and end points are marked by triangles and squares, respectively. A linear fit is shown for the dilatant proportional strain tests and repeated on the biaxial tests, as all the proportional strain tests converge towards a master curve.

Figure 3 Figure 3 Figure 3 Figure 3

 3333 Figure 3.11: p -r * 3 for proportional strain (a) and biaxial (b) test. The start and end points are marked by triangles and squares, respectively. An exponential fit is shown for the dilatant proportional strain tests and repeated on the biaxial tests, as all the proportional strain tests converge towards a master curve.

  Figure3.15: p -q evolution along combining proportional strain and biaxial loading paths. The proportional strain loading with λ = -1.2 are conducted at first until a certain axial strain level labeled as A, B, C, D, E and F, after that a biaxial loading is performed by keeping the corresponding lateral stress unchanged. The transition states and the average critical states are marked by dots and squares, respectively. The corresponding stressstrain responses are shown in Fig.3.16

Figure 3

 3 Figure 3.16: Stress-strain responses along combined proportional strain and biaxial loading paths. The switch states and ultimate states are marked by dots and squares, respectively. The corresponding stress paths are shown in Fig.3.15

Chapter 4 .

 4 Dynamical view of the critical state 5

Figure 4

 4 Figure 4.1: Evolution of deviatoric stress q and volumetric strain ε v with respect to the axial strain ε 2 along biaxial tests. Note that soil mechanics convention is adopted with positive compression and positive contraction.

Figure 4 .

 4 Figure 4.2: Incremental deviatoric strain (IDS) maps at the critical states, estimated for axial strain increments of 0.1%. Two samples are considered: (a) 100kPa-dense sample (ε 2 = 6.1%) and (b) 100kPa-loose sample (ε 2 = 8.02%).

Figure 4 . 3 :

 43 Figure 4.3: Proportions of chained particles.Two samples are considered: (a) 100kPa-dense sample and (b) 100kPa-loose sample. Note that r * cha represents the value when the spatial domain of shear band are considered, and r cha when the whole sample is considered. The deviatoric stresses are recalled.

Chapter 4 .

 4 Dynamical view of the critical state

Figure 4 . 4 :

 44 Figure 4.4: Volumetric strain of different types of loops and the global assembly during the biaxial tests in 100kPa-dense sample (a) and 100kPa-loose sample (b).

Figure 4 . 5 :

 45 Figure 4.5: Lifespan and life expectancy PDFs (probability density functions) of chained particles in two samples: (a) 100kPa-dense sample where only the spatial domain within the shear band is considered and (b) 100kPa-loose sample with the whole domain considered.For each sample, two reference conformations belonging to the stationary regime are adopted, labeled as M and N. The deviatoric stress curve (dot-dashed line) and the volumetric strain (dashed line) curves are recalled.

Figure 4 . 7 :

 47 Figure 4.7: The strain magnitudes within the shear band domain (100kPa-dense*), whole dense sample (100kPa-dense) and the whole loose sample (100kPa-loose). The strain magnitude ratio between 100kPa-dense* and 100kPa-dense is around 2.8.

Figure 4 . 8 :

 48 Figure 4.8: PDFs for chained particles from 100kPa-loose sample and from 100kPa-dense sample stretched by the strain magnitude ratio 2.8.

Figure 4 . 9 :

 49 Figure 4.9: Lifespan and life expectancy PDFs (probability density functions) of grain loops in two samples with respect tho the axial strain ε 2 : (a) 100kPa-dense sample where only the spatial domain within the shear band is considered and (b) 100kPa-loose sample with the whole domain considered. For each sample, two reference conformations belonging to the stationary regime are adopted. The deviatoric stress curve (dot-dashed line) and the volumetric strain (dashed line) curves are recalled in the subfigure ofLoop 4.

Figure 4 .

 4 Figure 4.10: The averaged PDFs of grain loops with respect to the incremental strain evolution (∆ε 2 ), over six microstructure conformations at different strain levels belonging to the critical regime (ε 2 = 0.0610, 0.0708, 0.0806, 0.0904, 0.1002 and 0.1100 for the 100kPa-dense sample and ε 2 = 0.0605, 0.0708, 0.0802, 0.0906, 0.1005 and 0.1106 for the 100kPa-loose sample).

Figure 4 .

 4 Figure 4.11: PDFs for grain loops from 100kPa-loose sample and from 100kPa-dense sample stretched by the strain magnitude ratio 2.8.

Figure 4 .

 4 Figure 4.12: Lifespan and life expectancy PDFs for chained particles under different confining pressures with samples 40kPa-loose, 100kPa-loose and 400kPa-loose considered.

Figure 4

 4 Figure 4.13: Lifespan and life expectancy PDFs for grain loops with respect to the incremental strain evolution (∆ε 2 ) under different confining pressures with samples 40kPa-loose, 100kPa-loose and 400kPa-loose considered.

Figure 4 .

 4 Figure 4.14: Double-exponential fits given for the chained particle PDFs with respect to the absolute value of incremental evolution |∆ε 2 | under different confining pressures.

Figure 4 .

 4 Figure 4.15: Double-exponential fits given for the loop PDFs with respect to the absolute value of incremental evolution |∆ε 2 | under different confining pressures.

Figure 4 .

 4 Figure 4.17: Color maps of the lifespan of grains loops for the reference axial strain ε 2 = 8.02%. Two thresholds are considered corresponding to the two characteristic strains: (a) 0.4% and (b) 3%.

Figure 4 .

 4 Figure 4.18: Color maps of the life expectancy of grains loops for the reference axial strain ε 2 = 8.02%. Two thresholds are considered corresponding to the two characteristic strains: (a) 0.4% and (b) 3%.

  particlesIn Fig.4.19, lifespan and life expectancy PDFs for chained particles are displayed for the initial state of the whole sample in both dense and loose cases.

Figure 4 .

 4 Figure 4.19: Life expectancy PDFs (probability density functions) of chained particles at the initial state of the whole sample in the dense (a) and loose (b) cases. The deviatoric stress curve (dot-dashed line) and the volumetric strain (dashed line) curves are recalled.

Figure 4 .Figure 4 .

 44 Figure 4.20: Life expectancy PDFs (probability density functions) of grain loops at the initial state of the whole sample in the dense (a) and loose (b) cases. The deviatoric stress curve (dot-dashed line) and the volumetric strain (dashed line) curves are recalled in the subfigure of Loop 4.

Figure 4 .

 4 26. Corresponding fits are given based on the double exponential equation 4.3. It can be seen that there is a narrow gap between the lifespan and life expectancy 0

Figure 4 .

 4 Figure 4.22: Life expectancy PDFs (probability density functions) of grain loops under the dilatant proportional strain loading considering the whole sample.

Figure 4 .

 4 Figure 4.23: Lifespan and life expectancy PDFs (probability density functions) of grain loops at the convergent state under the dilatant proportional strain loading considering a domain within the shear band.

Figure 4 .

 4 Figure 4.24: The six states chosen for color maps along the proportional strain tests

Figure 4 .Figure 4 .Figure 4 .

 444 Figure 4.25: The absolute value of PDFs for Lifespan and life expectancy of grain loops under the dilatant proportional strain loading within the domain of the shear band.

  Figure 5.1: The directional distributions and the coordinates at the REV scale (e 1 , e 2 ) and the mesoscale (n, t).

Figure 5 . 3 :

 53 Figure 5.3: General homogenization scheme relating both incremental stress and strain tensors

  Figure 5.5: Two modes of pathological H-cell due to contact loss.

Figure 5 . 6 :

 56 Figure 5.6: Two modes of over-compressed H-cell with the corresponding opening angle.

Figure 5 . 8 :

 58 Figure 5.8: The deviatoric stress and volumetric strain along the biaxial loading with p 0 = 100 kPa until a relative large axial strain.

  of the ratio between the length d 1 and its initial value d 0 : (a) directional analysis at different axial strains ε 11 = 0.0, 0.005, 0.032, 0.128 and 0.347 labelled as O, A, B, C and D, respectively, (b) macro strain evolution along specific directions. 10: Evolution of the ratio between the length d 2 and its initial value d 0 : (a) directional analysis at different axial strains ε 11 = 0.0, 0.005, 0.032, 0.128 and 0.347 labelled as O, A, B, C and D, respectively, (b) macro strain evolution along specific directions.

Figure 5 .Figure 5

 55 Figure 5.11: Porosity ϕ 1 meso of the hexagon with respect to the opening angle α; sketch of the two possible deforming processes.

Figure 5 .

 5 Figure 5.16: Flowchart of the computation of biaxial loading incorporating the deactivation/reactivation scheme.

Figure 5 .Figure 5 Figure 5 Figure 5

 5555 Figure5.17: Deviatoric stress and volumetric strain with respect to the axial strain using the updated H-model.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 6.1.1 Dilatant proportional strain loading and relationships with critical state123 6.1.2 Dynamical view of critical state . . . . . . . . . . . . . . . . . . . . . 124 6.1.3 Enriched H-model with emerging critical state . . . . . . . . . . . . . 125 6.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.2.1 Micromechanical investigations . . . . . . . . . . . . . . . . . . . . . 125 6.2.2 H-model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 6.2.3 Multiscale modeling of landslides . . . . . . . . . . . . . . . . . . . . 129

Figure 6 . 1 :

 61 Figure 6.1: Stress paths in p -q plane along triaxial tests and proportional strain tests with different dilatant rates characterized by κ (dε v = (1 + 2κ)dε 2 ).

  Figure6.2: Stress-strain responses along mixed proportional strain and triaxial loading paths in 3D. At the transition state, the lateral stress is kept unchanged. The loading of 100kPa-dense (κ = -0.7) as shown in Fig.6.1 is adopted here. εv = (1 + 2κ) ε1

  Figure 6.3: simulation of the collapse of granular column using SPH particles and the Hmodel: (a) initial configuration and (b) color map of the deformation in direction 1 after collapse.

Figure 6 . 4 :Figure A. 2 :

 642 Figure 6.4: Aerial view of the Bas-Verger landslid in 2006 with associated localization of the field measurements conducted on the landslide (Fressard et al. 2016).

  4.11 PDFs for grain loops from 100kPa-loose sample and from 100kPa-dense sample stretched by the strain magnitude ratio 2.8. . . . . . . . . . . . . . . . 4.12 Lifespan and life expectancy PDFs for chained particles under different confining pressures with samples 40kPa-loose, 100kPa-loose and 400kPa-loose considered. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.13 Lifespan and life expectancy PDFs for grain loops with respect to the incremental strain evolution (∆ε 2 ) under different confining pressures with samples 40kPa-loose, 100kPa-loose and 400kPa-loose considered. . . . . . . 4.14 Double-exponential fits given for the chained particle PDFs with respect to the absolute value of incremental evolution |∆ε 2 | under different confining pressures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Double-exponential fits given for the loop PDFs with respect to the absolute value of incremental evolution |∆ε 2 | under different confining pressures. . .

  Shearing behavior of granular materials . . . . . . . . . . . . . . . 5 2.1.1 Dilatancy and contractancy . . . . . . . . . . . . . . . . . . . . . . .

	Chapter 2
	A brief literature review
	Contents
	2.1

6 2.1.2 Softening and hardening . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.3 Critical State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Discrete element view of granular materials . . . . . . . . . . . . 11 2.2.1 Discrete Element Method . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 Micromechanical analysis in granular materials . . . . . . . . . . 15 2.3.1 Particle/contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.2 Force chains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.3 Grain loops . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.4 Multiscale modeling of granular materials . . . . . . . . . . . . . . 20

Table 3 .

 3 1: Initial void ratios e 0 of prepared samples. The sample '100kPa-dense' is used in proportional strain tests)

	Sample 40kPa-dense 60kPa-dense 100kPa-dense 10kPa-loose
	e 0	0.197	0.195	0.191	0.266
	Sample 40kPa-loose 60kPa-loose	100kPa-loose	20kPa-loose
	e 0	0.267	0.265	0.261	0.267

Table 4 .

 4 1: Initial void ratios e 0 of prepared samples.

	Sample 40kPa-loose 100kPa-loose 400kPa-loose
	e 0	0.267	0.261	0.234
	Sample 40kPa-dense 100kPa-dense 400kPa-loose
	e 0	0.197	0.191	0.164

Table 4 .

 4 2: The fitting parameters for different confining pressures.

	Sample	Parameters Chained particles	L3	L4	L5	L6p
	40kPa-loose	ε c1	0.25%	0.30% 0.33% 0.30% 0.25%
		ε c2	2.29%	1.61% 1.92% 1.81% 1.59%
		P c	0.165	0.229 0.228 0.288 0.290
		ε c1 /ε c2	9.02	5.44	5.78	5.97	6.41
	100kPa-loose	ε c1	0.35%	0.43% 0.41% 0.34% 0.34%
		ε c2	3.30%	2.57% 2.79% 2.40% 2.25%
		P c	0.146	0.205 0.162 0.209 0.271
		ε c1 /ε c2	9.35	5.92	6.89	7.03	6.71
	400kPa-loose	ε c1	0.61%	0.77% 0.64% 0.59% 0.50%
		ε c2	5.33%	5.56% 5.60% 4.35% 3.73%
		P c	0.129	0.124 0.113 0.169 0.209
		ε c1 /ε c2	8.76	7.25	8.74	7.41	7.52

  Contents 5.1 Review and analysis of the H-model . . . . . . . . . . . . . . . . . 97 5.1.1 The H-model in brief . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.1.2 Biaxial test at the material point scale . . . . . . . . . . . . . . . . . 102 5.1.3 Mesoscale inspection of the H-model during biaxial loading . . . . . 105 5.2 Emergence of critical state in the H-model . . . . . . . . . . . . . 114 5.2.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.2.2 Preliminary results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

Table 5 .

 5 1: Parameters used in biaxial tests Veylon 2017k n (N/m) k t /k n ϕ g α 0 α ω β ω

	p 0 (kPa)

v Figure 5.7: The stress ratio and volumetric strain curves based on the parameters from Veylon 2017. Three confining pressures are considered.

In the present study, the out of plane dimension is then taken equal to the thickness of the sample box (0.04 m), which allows kPa rather than kN/m to be used for the stress unit on the boundary walls. The thickness of the box has a proportional influence on the stress values, which should not impact the following results, at least qualitatively.

Note that the terms unconstrained and imposed are used to recall that the volumetric strain is a response variable in a biaxial loading while it is a control variable in proportional strain tests.

Acknowledgements

First and foremost, I would like to express my sincere gratitude to Dr. François Nicot and Dr. Antoine Wautier for their clear guidance and continuous encouragement during the three years. I really appreciate the wonderful research environment they created in our group through organizing in-group and out-group meetings regularly. I have learnt so much from them, from technical methods to critical thinking skills. I have been inspired

to keep an open mind and stay curious. I know the experience working with them will

Complementary results on critical state dynamics

In addition to the results presented in the main part of the chapter, figures A.1 and A.2 provides complementary results on a) the detailed balance achieved at critical state for varying confining pressures and b) the uniqueness of the PDF for lifespan and life expectancy in dense and loose samples provided the strain rate is rescaled in the dense case to account for strain localization.