Enn

Phuong Thao

Ls Vinh

Vinh Thân

Du

Keywords: DEM, wet granular materials, one-dimensional compression, isotropic compression, X-ray computed tomography, grain-scale analysis, plastic behavior, microstructure v Modélisation des élements discrets, matériaux granulaires mouillés, compression oedométrique, compression isotrope, microtomographie aux rayons X, analyse à l'échelle du grain, comportement plastique, microstructure PVA, polyvinyl alcohol, ESEM, Environmental Scanning Electron Microscopy, DIC, Digital Image Correlation, V-DIC, Volumetric Digital Image Correlation, *

Acknowledgment

Ce travail de thèse n'aurait être possible sans appui de beaucoup de gens. Je tiens à remercier dans cette page un certain nombre de peronnes qui m'ont aidé à mener à bien ce travail.

Je tiens tout d'abord à remercier Monsieur Gaël Combe, qui m'a fait l'honneur de présider le jury de ma soutenance de thèse ainsi qu'a examiné mon travail. Ses remarques et questions sont très intéressantes et constructives.

Je voudrais remercier vivement Messieurs

Jean-Yves Delenne et Bertrand François qui ont accepté la longue et lourde tâche d'être rapporteurs de ma thèse. Leurs remarques fructueuses et présieuses m'ont aidé à améliorer le niveau scientique de ce travail. Je souhaite exprimer ma profonde gratitude à mes trois directeurs de thèse Anh Minh Tang, Jean-Noël Roux et Jean-Michel Pereira pour leur conance, leur suivi et leur gentilesse. Merci à Jean-Michel pour tes conseils précieux et ton aide. Merci sincèrement à Jean-Noël pour ton modèle de la simulation, tes conseils précieux, tes encouragements ainsi que ta patience dans les moments que je n'ai pas pu bien comprendre ton modèle. Si ce travail est une réussite, c'est une très grande partie grâce à toi. Un grand merci à Minh pour avoir accepté ma demande à faire ma thèse à l'École des Ponts ParisTech. C'est un grand honneur pour moi. Tu m'as beaucoup aidé, tu m'as guidé et accompagné tout le long de l'élaboration de ce travail. Tes soutiens et encouragements m'ont permis mener à bien ce travail. Je voudrais également adresser mes remerciements aux Messieurs Michel Bornert et Patrick Aimedieu de l'équipe multi-échelles, qui m'ont aidé de faire des essais de microtomographie aux rayons X. Leurs conseils sont très utiles

2.1

Cumulative grain size distribution of glass beads . . . . . . . . . . . . . .

2.2

Specimen preparation for cell φ50 . . . . . . . . . . . . . . . . . . . . . .

2.3

Specimen preparation for cell φ20 . . . . . . . . . . . . . . . . . . . . . .

2.4

Schematic view of the system for preparing the specimen . . . . . . . . .

2.5

Initial void ratio e 0 and initial porosity n versus height of free fall h f . .

2.6

Schematic view of oedometric compression system . . . . . . . . . . . . .

3.1

Grain size distribution of the polydisperse system . . . . . . . . . . . . .

3.2

Static normal force, F N = F e N + F cap , as a function of the interparticle distance h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.3

Coulomb cone limiting the value of the tangential force . . . . . . . . . .

3.4

Motion of two spheres . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.5

Three spheres connected with liquid bridges . . . . . . . . . . . . . . . .

3.6

Capillary force law F cap (h), for two dierent meniscus values, according to the Maugis approximation and the Soulié formula . . . . . . . . . . .

3.7

Monodisperse system on FCC lattice . . . . . . . . . . . . . . . . . . . .

3.8

Typical specimens of (a) the monodisperse (Φ 0 = 0.30) and (b) the polydisperse (Φ 0 = 0.24) systems. . . . . . . . . . . . . . . . . . . . . . . . .

Overview

Granular materials appear in various natural and industrial processes. They are composed of individual grains that vary in shape, size, composition, surface texture, etc.

These granular characteristics can remarkably aect the packing and the contact distribution properties of granular assemblies. Therefore, granular materials may exhibit quite dierent behaviors depending on the applied external forces, the size and geometry of the grains, the density and composition of particles, the type of interactions at the grain scale and many other material properties. Because of these complexities, it is generally dicult to relate in a straightforward manner the macroscopic properties of granular materials to the microscopic properties of their constituents. In the case of wet granular materials, the existence of liquid menisci between particles plays a key role in the overall behavior of the assembly. Capillary cohesion bestows to these materials specic mechanical features that do not exist with dry grains, such as the ability to form stable structures with very low density, and a strong sensitivity to stress intensity as well as to stress direction. So far, many experimental and numerical studies of bonded granular materials in general and wet granular materials in particular have been carried out to tackle the diculties exposed here above.

Published studies cover a wide range of applications. For instance, some investigate the mechanical properties of cohesive soils (clays and silts) [START_REF] Wood | Soil Behaviour and Critical State Soil Mechanics[END_REF]2,[START_REF] Atkinson | An Introduction to Mechanics of Soils and Foundations: Through Critical State Soil Mechanics[END_REF][START_REF] Tripodi | Elastoplastic nite element model development and validation for low pressure uniaxial compaction of dry cohesive powders[END_REF][START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF], metallic powder processing [START_REF] Poquillon | Cold compaction of iron powders -relations between powder morphology and mechanical properties: Part i: Powder preparation and compaction[END_REF] or modeling and treatment of ceramic powders [START_REF] Cooper | Compaction behavior of several ceramic powders[END_REF][START_REF] Reed | Principles of Ceremic Processing[END_REF][START_REF] Falgon | Identication procedure of a hardening law for powder compaction[END_REF]. Assemblies of wet beads were observed [START_REF] Pierrat | Tensile strength of wet granula materials[END_REF][START_REF] Gröger | Modelling and measuring of cohesion in wet granular materials[END_REF][START_REF] Fournier | Mechanical properties of wet granular materials[END_REF], and descriptions of such materials at microscopic scale were made possible by using numerical and experimental tests [START_REF] Gröger | Modelling and measuring of cohesion in wet granular materials[END_REF][START_REF] Fournier | Mechanical properties of wet granular materials[END_REF]. A 2D numerical investigation into the microstructure and mechanical properties of cohesive powders is reported in [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF][START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF]. In addition, some numerical studies of cohesive materials have shown that the stability of loose structures formed by particles packed under gravity rely on adhesive forces (see for instance Dong et al. [START_REF] Dong | Role of interparticle forces in the formation of random loose packing[END_REF]). Other numerical works studied loose packings stabilized by cohesion and their collapse when subjected to increasing loads during oedometric compression [START_REF] Kadau | Contact dynamics simulations of compacting cohesive granular systems[END_REF][START_REF] Kadau | Pore stabilization in cohesive granular systems[END_REF][START_REF] Wolf | Compaction of cohesive powders[END_REF][START_REF] Bartels | The eect of contact torques on porosity of cohesive powders[END_REF]. Some studies focused on the failure of bounded particle assemblies in static [START_REF] Delenne | Mechanical behaviour and rupture of cohesive granular media[END_REF][START_REF] Delenne | Mechanical behaviour and failure of cohesive granular materials[END_REF] or dynamic conditions [START_REF] Thornton | How do agglomerates break?[END_REF]. Other works on wet bead packs in which cohesion stems from liquid bridges joining neighbouring particles investigated the structure of poured samples [START_REF] Yang | Numerical study of the packing of wet coarse uniform spheres[END_REF], or their shear strength [START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF].

2 Overview Shear ow of cohesive granular materials has also been simulated [START_REF] Brewster | Plug ow and the breakdown of bagnold scaling in cohesive granular ows[END_REF][START_REF] Aarons | Shear ow of assemblies of cohesive and non-cohesive granular materials[END_REF][START_REF] Rognon | Rheophysics of cohesive granular materials[END_REF][START_REF] Rognon | Dense ows of cohesive granular materials[END_REF][START_REF] Khamseh | Flow of wet granular materials: a numerical study[END_REF].

Besides, the collapsible behavior of loessic soils [START_REF] Muñoz Castelblanco | Some aspects of the compression and collapse behaviour of an unsaturated natural loess[END_REF][START_REF] Muñoz Castelblanco | The water retention properties of a natural unsaturated loess from northern france[END_REF][START_REF] Jiang | Summary of collapsible behaviour of articially structured loess in oedometer and triaxial wetting tests[END_REF], as well as wet sand with capillary bridges [START_REF] Mason | Critical angle of wet sandpiles[END_REF][START_REF] Mitarai | Wet granular materials[END_REF][START_REF] Bruchon | Full 3d investigation and characterisation of capillary collapse of a loose unsaturated sand using x-ray ct[END_REF] has been investigated experimentally. The contribution of articial solid bridges to the micro-macroscopic behavior of soils, especially in the case of cemented sands, has been studied using both numerical simulations and experimental observations [START_REF] Wang | A particulate-scale investigation of cemented sand behavior[END_REF][START_REF] Wang | Characterization of cemented sand by experimental and numerical investigations[END_REF][START_REF] Jiang | Modeling shear behavior and strain localization in cemented sands by two-dimensional distinct element method analyses[END_REF][START_REF] Jiang | A simple distinct element modeling of the mechanical behavior of methane hydrate-bearing sediments in deep seabed[END_REF][START_REF] De Bono | Discrete element modelling of one-dimensional compression of cemented sand[END_REF][START_REF] De Bono | Dem of triaxial tests on crushable sand[END_REF][START_REF] Mcdowell | On the micro mechanics of one-dimensional normal compression[END_REF]. These works evidenced the link between the macroscopic mechanical behavior and changes in the microstructure. However, most studies on wet granular materials have focused on dense states, and few works have focused on the grain-scale behavior.

Currently, the grain-scale behavior of wet granular materials is usually investigated using numerical simulations, such as the discrete element method [START_REF] Cundall | A disrete numerical model for granular assemblies[END_REF], in 2D and/or 3D [START_REF] Delenne | From liquid to solid bonding in cohesive granular media[END_REF][START_REF] Delenne | Liquid clustering and capillary pressure in granular media[END_REF][START_REF] Kadau | Pore stabilization in cohesive granular systems[END_REF][START_REF] Khamseh | Flow of wet granular materials: a numerical study[END_REF][START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF][START_REF] Richefeu | Force transmission in dry and wet granular media[END_REF][START_REF] Rognon | Rheophysics of cohesive granular materials[END_REF][START_REF] Than | Basic mechanical properties of wet granular materials: A DEM study[END_REF]. These studies characterized precisely the microstructure in terms of coordination number of contacts, coordination number of distant interactions, coordination number of compressive and tensile bonds, radial distribution function, force chains, distribution of forces between particles, etc. The evolution of these microstructural descriptors with the externally applied pressure was also investigated. However, few works have quantitatively compared simulations (e.g. with the discrete element method)

with experimental results (e.g. with microstructural observations using the X-ray computed tomography) in order to validate the numerical method [START_REF] Fu | Investigation of particle packing in model pharmaceutical powders using x-ray microtomography and discrete element method[END_REF][START_REF] Golchert | Application of X-ray microtomography to numerical simulations of agglomerate breakage by distinct element method[END_REF][START_REF] Moreno-Atanasio | Combining X-ray microtomography with computer simulation for analysis of granular and porous materials[END_REF][START_REF] Wang | A particulate-scale investigation of cemented sand behavior[END_REF][START_REF] Wang | Characterization of cemented sand by experimental and numerical investigations[END_REF].

Following the two-dimensional model for cohesive powders which is based on the discrete element method and developed by Gilabert et al. [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF][START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF], we propose a threedimensional model to investigate the mechanical behavior of a wet granular material at very loose states of density. On the same material, we also perform one-dimensional compression tests combined with microstructural observations using X-ray computed tomography. Thereby, we provide a comprehensive view of the mechanical behavior, as well as a further insight into microstructure changes of wet granular materials at very loose states under growing applied external forces. In the present study, we only focus on the properties of wet granular materials in the pendular state.

This thesis is organized in four chapters, as follows:

In Chapter 1, the mechanical behavior of granular materials is briey described. Several typical experimental and numerical observations of wet granular materials are then presented in details. The combinations of experiments and simulations which have been studied in literature are also presented in order to provide a comprehensive view of wet granular materials.

In Chapter 2, we present two experimental tests on wet spherical beads at very loose states. The behavior is characterized at the macroscopic scale as well as at the grain (and pore) scale, by combining one-dimensional compression tests and X-ray computed tomography.

In Chapter 3, following the 2D model from Gilabert et al. [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF][START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF], we propose a 3D model to investigate the mechanical behavior of wet spherical glass beads at very loose Overview 3 states under isotropic compression. Three systems, including monodisperse systems without and with rolling and pivoting resistances and a polydisperse system without rolling and pivoting resistances, are simulated. The plastic response of the material, as well as the evolution of microstructure and force transmission, along compression and decompression paths are then characterized. The inuence of various micromechanical parameters have also been investigated.

In Chapter 4, we rst propose a 3D observation of wet spherical glass beads in an oedometrically compressed model. We then validate the compaction behavior and microstructural changes by comparing the results of experiments and simulations.

Finally, some conclusions and perspectives are drawn.

CHAPTER 1

Literature review

In this chapter, the mechanical behavior of granular materials is briey described, with a particular focus on wet granular materials. Several non-destructive techniques to observe the micro-macroscopic behavior of such materials are then presented. Numerical investigations, which provide a detailed access to the microstructure of these materials, are then described. A comparison of experimental and numerical results are nally presented. Granular materials are ubiquitous in various elds, and natural or industrial processes including soft matter physics, soil mechanics, powder technology, agronomic transformations, and geological processes. They are composed of individual grains that vary in shape, size, and surface texture. Because of the wide variety of physicochemical and morphological grain properties, the fundamental behavior of granular materials should be investigated under dierent working environments. The analysis of physical systems involving granular materials requires a clear understanding of their behavior not only at the single particle level, but should also consider multi-physics problems involving multiscale phenomena from molecular to macroscopic scales. A brief description of granular materials in general is now given.

Satake and Tobita [START_REF] Iwashita | Mechanics of Granular Materials: An Introduction[END_REF] dened granular materials as grains in contact and the surrounding voids. The micromechanical behavior of granular materials is therefore inherently discontinuous and heterogeneous. The macroscopic (overall or averaged) behavior of granular materials is determined not only by how discrete grains are arranged in space, but also by what kinds of interactions are operating among them. In order to understand the mechanical behavior of granular materials from a microscopic point of view, the spatial distribution and orientation of grains and their contact conditions should be rst specied.

Lanier [START_REF] Lanier | Mécanique des milieux granulaires[END_REF] described granular materials as an intermediate class of materials between uids and solids, that not only ow like uids (snow avalanches, emptying silos, etc.) but also resist to deviatoric stresses (like solids). This type of materials is in the category of complex structured materials, anisotropic, and strongly heterogeneous. Their mechanical behavior depends on the interactions occurring at the particle level. Noteworthy is the fact that this behavior is partially reversible only within a very restricted domain and that it is non-linear.

After Radjaï and Lanier (in Cambou et al. [START_REF] Cambou | Micromechanics of granular materials[END_REF]), granular materials consist of densely packed solid particles and a pore-lling material which can be a uid or a solid matrix.

The particles interact via elastic repulsion, friction, adhesion and other surface forces.

By nature, the length scales involved in these contact interactions are much smaller than the particle size. External loading leads to particle deformations as well as cooperative particle rearrangements. For instance, the particle deformations are of particular importance in powder metallurgy but the particles may be considered as quasi-rigid beyond the elastic deformation.

From a macroscopic scale point of view, the phenomena that make granular materials interesting are pattern formation, mixing or segregation, clustering (granular gases), avalanches, rotating ows, granular convection, and jamming/ unjamming. For example, Chapter 1. Literature review if a granular material is heaped on an inclined plane, then the large scale state of the system depends on the angle of the plane. For large angles of this plane, the granular material ows like a non-Newtonian liquid. For small angles of this plane, the granular material will behaves like a solid and remain stationary. The critical value of the plane angle delimiting the two regimes (or phases) depends on the preparation history, and on the transition between the phases. The fundamental behavior of granular materials has been widely studied in the literature [START_REF] Antony | Granular Materials: Fundamentals and Applications[END_REF]. The following features have been investigated:

(i) quasi-static deformation characteristics of granular materials at both small and large strains, with discrete particles considered in the analysis;(ii) eects of mechanical periodic excitation on a granular medium; (iii) advances in the constitutive modeling of granular materials using a continuum approach; (iv) interactions of particles at elevated temperatures, including adhesion forces and sintering; (v) critical state starting from distinct initial conditions using experiments and computer simulations. In summary, there is a lot of interesting and complex physics to be understood by studying just the large scale properties of granular assemblies. As the composition of the assembly becomes more complicated, its behavior becomes even richer.

From a meso-microscopic scale point of view, Radjaï and Lanier (in Cambou et al. [START_REF] Cambou | Micromechanics of granular materials[END_REF]) noted that the geometrical changes of granular texture were at the origin of the complex rheology of granular materials. These changes were highly nonlinear, involving creation and loss of contacts, rotation frustration and frictional sliding. They depended on the dissipative nature of contact interactions and steric exclusions among particles. In quasi-static deformation, various features of the plastic behavior such as shear strength and dilatancy could be traced back to the evolution of the granular texture. Therefore, these authors proposed to analyze and distinguish the variables as follows: granular texture, kinematics, and force transmission in granular materials.

The granular texture is disordered with many dierent variants depending on the composition (particle shapes and sizes), interactions and assembling procedure. The granular texture evolves with loading. The inhomogeneous distribution of contact forces reects granular disorder in static equilibrium. At the lowest order, the relevant scalar parameters concern the connectivity of this network. At higher orders, the anisotropy of the texture is described by fabric tensors. The granular texture evolves mainly due to contact losses and gains. The fraction of lost and gained contacts depends on the contact orientation. During the deformation process, the gained contacts tend to move in the direction of the major principal stress, conversely, the lost contacts are generally observed in the direction of the minor principal stress [START_REF] Oda | The mechanism of fabric changes during compressional deformation of sand[END_REF][START_REF] Antony | Granular Materials: Fundamentals and Applications[END_REF].

The kinematics are directly linked to the deformation of granular media. In other words, the kinematics are the displacement of particles. Each grain has six degrees of freedom (three components of translation and three components of rotation). As an example, the analysis of the displacement eld is shown in Fig. 1.1 [START_REF] Andò | Grain-scale experimental investigation of localised deformation in sand: A discrete particle tracking approach[END_REF]. A triaxial compression test was carried out and combined with the X-ray microtomography to get The contact forces are considered like the static variable of granular assemblies. These forces appear when the granular materials are subjected to an external load. The contact network and the corresponding contact forces result from the reaction of the granular material to the applied load and the transmission of internal forces. The inhomogeneity of these forces in granular assemblies were rst observed optically in assemblies of photoelastic particles, which have the property to develop birefringence upon the application of stresses [START_REF] Dantu | Contribution à l'étude mécanique et géométrique des milieux pulvérulents[END_REF], as shown in Fig. 1.2. This precursory experimental study allowed to observe the load bearing contacts and the establishment of force chains within a granular material. Since then, other studies have tried to use this technique to measure forces in such granular assemblies [START_REF] Jaeger | Granular solids, liquids, and gases[END_REF]. However, it is dicult to measure accurately tangential forces experimentally.

Later, numerical simulations enabled to analyze and provide detailed evidence of force chains, the classication of force networks in strong and weak networks, and the exponential distribution of strong forces [START_REF] Radjaï | Friction-induced self-organization of a one-dimensional array of particles[END_REF][START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]. Moreover, the force probability density Chapter 1. Literature review function from simulation showed that the weak forces (below the average force) in a sheared granular system have a nearly uniform or decreasing power-law shape which is in agreement with experiments [START_REF] Mueth | Force distribution in a granular medium[END_REF][START_REF] Radjaï | Force distributions in dense two-dimensional granular systems[END_REF].

Bonded granular materials

In general, two types of bonding are distinguished in the literature: solid bonding and liquid bonding. While liquid bonding arises from the presence of capillary meniscii at the grain contacts, solid bonds between grains may originate from dierent sources, such as the process of sedimentation in natural soil and rock deposits or the precipitation of a solid, having either a natural or an articial origin. During sedimentation, cementation is often observed during the early diagenesis process [START_REF] Pettijohn | Sand and sandstone[END_REF]. In the case of solid precipitation, on the one hand, natural cementation takes various forms such as calcite, silica, iron oxides, or even clays [START_REF] Wood | Soil Behaviour and Critical State Soil Mechanics[END_REF]2,[START_REF] Atkinson | An Introduction to Mechanics of Soils and Foundations: Through Critical State Soil Mechanics[END_REF][START_REF] Tripodi | Elastoplastic nite element model development and validation for low pressure uniaxial compaction of dry cohesive powders[END_REF][START_REF] Santamarina | Soils and Waves: Particulate Materials Behavior, Characterization and Process Monitoring[END_REF][START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF]. On the other hand, articial bonded materials are often encountered in applications involving materials improved by mixing with Portland cement or lime.

Despite the complicatedness of the formation of the bonds, the eects of general features of solid bonding on the properties of granular materials were widely observed based on experimental ndings: (i ) the strength (dynamic and static) and small-strain stiness are enhanced [START_REF] Dupas | Static and dynamic properties of sand-cement[END_REF][START_REF] Acar | Low strain dynamic properties of articially cemented sand[END_REF][START_REF] Clough | Cemented sands under static loading[END_REF][START_REF] Saxena | Dynamic moduli and damping ratios for cemented sands at low strains[END_REF][START_REF] Huang | Properties of articially cemented carbonate sand[END_REF]; (ii ) the stress-strain and volumetric response become relatively brittle and more dilative, respectively [START_REF] Clough | Cemented sands under static loading[END_REF][START_REF] Abdulla | Behavior of Cemented sands -I. Testing[END_REF][START_REF] Schnaid | Characterization of cemented sand in triaxial compression[END_REF]; and (iii ) quasipreconsolidation pressure or yield stress can be observed in loading path [START_REF] Leroueil | The general and congruent eects of structure in natural soils and weak rocks[END_REF][START_REF] Airey | Triaxial testing of naturally cemented carbonate soil[END_REF]. Furthermore, the mechanical responses of cemented soils are found to depend on the amount and nature of the cementing agent [START_REF] Ismail | Eect of cement type on shear behavior of cemented calcareous soil[END_REF]. At the macroscopic scale, cementation induces a strength enhancement and the occurrence of volumetric dilation [START_REF] Wang | A particulate-scale investigation of cemented sand behavior[END_REF][START_REF] Wang | Characterization of cemented sand by experimental and numerical investigations[END_REF]. In addition, at the grain scale, the contact behavior of bonded materials was also observed by several studies using the idealized granules bonded [START_REF] Jiang | Contact behavior of idealized granules bonded in two dierent interparticle distances: An experimental investigation[END_REF][START_REF] Jiang | An experimental investigation on the mechanical behavior between cemented granules[END_REF], and the articially structured bonds [START_REF] Jiang | Summary of collapsible behaviour of articially structured loess in oedometer and triaxial wetting tests[END_REF].

Furthermore, liquid bonds were also observed in several experimental investigations.

The eect of capillary bonding on the mechanical behavior of granular materials is of primary importance in powder technology [START_REF] Iveson | Fundamental studies of granule consolidation part 1: Eects of binder content and binder viscosity[END_REF][START_REF] Iveson | Nucleation, growth and breakage phenomena in agitated wet granulation processes: a review[END_REF][START_REF] Nokhodchi | An overview of the eect of moisture on compaction and compression[END_REF] and transformations of geomaterials [START_REF] Cardell | Saltinduced decay in calcareous stone monuments and buildings in a marine environment in sw france[END_REF][START_REF] Lubelli | The role of sea salts in the occurrence of dierent damage mechanisms and decay patterns on brick masonry[END_REF][START_REF] Lourens | Salt crystallization as damage mechanism in porous building materials -a nuclear magnetic resonance study[END_REF]. The eects of inter-grain cohesion were also observed in order to present a general observation of the mechanical properties of wet granular materials [START_REF] Fournier | Mechanical properties of wet granular materials[END_REF][START_REF] Mitarai | Wet granular materials[END_REF].

Finally, the dynamics of wet granular matter were clearly reported in the study by Herminghaus [START_REF] Herminghaus | Dynamics of wet granular matter[END_REF].

Besides, bonded granular materials have been extensively investigated using numerical simulations in various contexts, such as: metallic powder processing [START_REF] Poquillon | Cold compaction of iron powders -relations between powder morphology and mechanical properties: Part i: Powder preparation and compaction[END_REF][START_REF] Poquillon | Cold compaction of iron powdersrelations between powder morphology and mechanical properties: Part ii. bending tests: results and analysis[END_REF], modeling and treatment of ceramic powders [START_REF] Cooper | Compaction behavior of several ceramic powders[END_REF][START_REF] Reed | Principles of Ceremic Processing[END_REF][START_REF] Falgon | Identication procedure of a hardening law for powder compaction[END_REF]. Assemblies of wet beads were observed [START_REF] Pierrat | Tensile strength of wet granula materials[END_REF][START_REF] Gröger | Modelling and measuring of cohesion in wet granular materials[END_REF][START_REF] Fournier | Mechanical properties of wet granular materials[END_REF], in which some microscopic observations are possible by combining numerical and experimental tests [START_REF] Gröger | Modelling and measuring of cohesion in wet granular materials[END_REF][START_REF] Fournier | Mechanical properties of wet granular materials[END_REF]. A 2D numerical investigation into the microstruc-ture and mechanical properties of cohesive powders was carried out [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF][START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF]. Furthermore, the micro-macroscopic behavior of granular materials containing articial solid bridges, especially cemented sands, was extensively investigated [START_REF] Wang | A particulate-scale investigation of cemented sand behavior[END_REF][START_REF] Wang | Characterization of cemented sand by experimental and numerical investigations[END_REF][START_REF] De Bono | Discrete element modelling of one-dimensional compression of cemented sand[END_REF][START_REF] De Bono | Dem of triaxial tests on crushable sand[END_REF][START_REF] De Bono | Investigating the micro mechanics of cemented sand using dem[END_REF][START_REF] Mcdowell | On the micro mechanics of one-dimensional normal compression[END_REF][START_REF] Jiang | An investigation on loose cemented granular materials via dem analyses[END_REF][START_REF] Jiang | Study of mechanical behavior and strain localization of methane hydrate bearing sediments with dierent saturations by a new dem model[END_REF].

In addition, the eects of liquid bonds on the mechanical behavior of granular matter have been extensively studied in the past, and several models of capillary cohesion have been proposed [START_REF] Mikami | Numerical simulation of cohesive powder behavior in a uidized bed[END_REF][START_REF] Willett | Capillary bridges between two spherical bodies[END_REF][START_REF] Soulié | Inuence of liquid bridges on the mechanical behaviour of polydisperse granular materials[END_REF]. The contribution of capillary bridges to the micromacroscopic mechanical behavior has been investigated by means of numerical simulations [START_REF] Kim | Modeling of tensile strength on moist granular earth material at low water content[END_REF][START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF][START_REF] Richefeu | Force transmission in dry and wet granular media[END_REF][START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF][START_REF] Scholtès | Micromechanics of granular materials with capillary eects[END_REF][START_REF] Guo | Capillary interaction-induced rolling resistance between elliptical particles and its inuence on grain column length at pendular state[END_REF][START_REF] Khamseh | Flow of wet granular materials: a numerical study[END_REF][START_REF] Than | Basic mechanical properties of wet granular materials: A DEM study[END_REF].

The most relevant known features of bonded granular materials at the macroscale are the yield criterion, the stress-strain response, the critical state, the inuence of cohesion, etc. Here, we present the inuence of intergranular cohesion on the overall behavior.

Macroscopic inuence of cohesion

A non-null cohesion in the Mohr-Coulomb failure criterion has an important consequence since it introduces a stress scale in the behavior of the material. While noncohesive materials are essentially sensitive to ratios of stress components, cohesive ones respond dierently according to the absolute magnitude of stresses. This is quite remarkable when dealing with compression tests in which all stress ratios remain constant.

During an oedometric compression test, the vertical stress σ 1 increases, while the radial strain is maintained null. Hence, the stresses are also observed to correspond to coecient of lateral stress ratio σ 3 /σ 1 . In addition, it is well recognized that the density of non-cohesive materials is slightly evolve with growing stresses, while cohesive ones undergo some irreversible compression under growing average pressure P. In soil mechanics [2], the irreversible compression curves are usually described using a linear relationship between the void ratio e (i.e. e = -1 + 1/Φ with Φ = 1 -n, Φ is the solid fraction and n is the porosity) and the logarithm of P:

e = e 0 -C c log P P 0 , (1.1)
in which, e 0 and P 0 are the initial void ratio and initial pressure, respectively. This relation applies to cohesive materials such as clays and silts [START_REF] Wood | Soil Behaviour and Critical State Soil Mechanics[END_REF][START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF] and also to cohesive powders with rather wide intervals of void ratio (e varying normally by a factor of 2).

It was also applied to critical states, characterized by a specic value of the stress ratio σ 3 /σ 1 , and the compaction curve is often modeled with the same slope C c (called the compression index). Relation (1.1) only applies for the monotonically growing pressures.

The irreversibility of the plastic compression phenomenon implies that it does not apply to pressure values lower than the highest pressure reached in the past, namely the preconsolidation pressure. The reports on bonded materials are similar but the change of density is considerably lower.

Wet granular materials

The mechanical properties of granular materials are remarkably aected when a small amount of water is added. For instance, a wet sand can be sculptured into quite stable structures, for instance a sand castle, while this is impossible when it is dry. However, this phenomenon signicantly depends on the size of particles, as shown in Fig. 1.

(van der

Waals attraction is computed for an inter-particle separation of 30 Å, the skeletal force is shown for σ = 10 kPa and σ = 1 MPa) [START_REF] Santamarina | Soil Behavior at the Microscale: Particle Forces[END_REF]. These compressive forces mobilize the electrical repulsion forces and bring particles together until compression and repulsion are balanced. Changing the pore uid can alter the inter-particle distance at equilibrium; the upper part of the gure shows the strain caused by changing the pore uid from fresh-water to seawater concentrations. of the gure shows the strain (axis on right) caused by changing the pore uid ionic concentration from fresh-water to seawater conditions. Note slopes: skeletal 2:1, weight 3:1, capillary and van der Waals 1:1 (from Santamarina [START_REF] Santamarina | Soil Behavior at the Microscale: Particle Forces[END_REF]).

Here, several observations can be made as follows (from Santamarina [START_REF] Santamarina | Soil Behavior at the Microscale: Particle Forces[END_REF]):

• Particle weight looses relevance with respect to capillary forces for particles smaller than d ≈ 3 mm (point 1), and with respect to van der Waals attraction for particles smaller than d ≈ 30 µm (point 2).

• Capillary forces can exceed the contribution of σ = 10 kPa connement for particles smaller than d ≈ 20 µm (point 3) and the contribution of σ = 1 MPa for d < 0.2 µm (point 4).

• Particles are considered coarse when skeletal forces due to boundary loads prevail. This is the case for particles larger than d ≈ 20 µm (point 3).

1.1. Mechanical behavior of granular materials 13

• Particles are fine when contact-level capillary and electrical forces gain relevance. This is the case when particles are smaller than d ≈ 1 -10 µm.

Therefore, we can see that the smaller the particle size, the stronger the capillary forces. Even the moisture in the atmosphere may create tiny liquid bridges at the contact point between particles. The presence of these liquid bridges induces capillary forces that attract the particles to each other. In reality, because of the asperities on the particles surfaces, molecular interactions scale down to small adhesive forces, and capillary forces can be large and dominant.

Unlike in the nanoscale range of adhesive forces, they can apply from an interparticle distance of the order of the particle size. Such unique properties of capillary interactions lead to remarkable changes in mechanical properties of granular materials in the appearance of the liquid in their interstitial spaces. Therefore, the presence of liquid menisci plays a key role in the overall mechanical properties of wet granular materials. Based on the liquid content, wet granular materials can be classied in dierent regimes [START_REF] Newitt | A contribution to the theory and practice of granulation[END_REF][START_REF] Iveson | Nucleation, growth and breakage phenomena in agitated wet granulation processes: a review[END_REF][START_REF] Mitarai | Wet granular materials[END_REF], as given in Table . 1.1.

Table 1.1: Schematic diagrams for classication of wet granular materials with various amount of liquid. In the third column, the black circles represent the grains and the grey regions represent the interstitial liquid (from Mitarai and Nori [START_REF] Mitarai | Wet granular materials[END_REF]).

• Pendular state: When a small amount of liquid is added to the granular material, the liquid initially collect near the contact point of particles. In this regime, the Chapter 1. Literature review liquid content is limited to the domain where menisci can form without coalescence; in this domain, the capillary forces are limited to pair-wise interactions.

• Funicular state: As the liquid content is increased further, the neighbor liquid bridges start to coalesce. The liquid lls some pores and multiple grains can be in contact with a given volume of liquid.

• Capillary state: At higher values of liquid content most of the pores lls with liquid and large contiguous wet clusters forms.

• Slurry state: In this regime all the pore space are fully saturated and the particles are completely immersed in the liquid.

In the rst three regimes, the role of capillary force is very important while it is negligible in the slurry state.

Many numerical and experimental studies reported that the presence of a small amount of interstitial liquid strongly aects the yield stress of the materials. The effects of the presence of water on the apparent Coulomb cohesion was widely studied in the literature. Richefeu et al. [START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF] studied the Coulomb cohesion in pendular state for dierent values of liquid content by using both numerical and experimental observations.

They observed that the Coulomb cohesion increases nonlinearly with liquid content and saturates to a maximum value, c m 600 Pa at w m 0.03. They found a constant value of friction angle ϕ (=tan -1 µ), regardless the level of liquid content.

In soil mechanics, Terzaghi [START_REF] Terzaghi | Principles of soil mechanics: a summary of experimental studies of clay and sand[END_REF] proposed the concept of eective stress to describe the mechanical properties of water-saturated soils. Terzaghi's principle states that all measurable eects of a change of stress of the soil, that is compression, distortion and change of shearing resistance, are exclusively due to changes in eective stress [START_REF] Nuth | Eective stress concept in unsaturated soils: Claracation and validation of a unied framework[END_REF].

Eective stress σ ij is dened as

σ ij = σ ij -u w δ ij , (1.2)
with u w is the pore water pressure and the Kronecker's delta δ ij : δ ii = 1, δ i =j = 0.

In the early works of Bishop [START_REF] Bishop | The principle of eective stress[END_REF], this principle was extended to unsaturated soils, as:

σ ij = (σ ij -u a δ ij ) + χ(u a -u w )δ ij , (1.3)
where χ is called the eective stress parameter or Bishop's parameter, vary from 0 for dry soils to 1 for saturated soils, and u a is the pore air pressure. The terms σ ij -u a δ ij = σ net-ij and u a -u w = s, dene the net stress and the matric suction, respectively. Independently of the stress state, the capillary forces are assumed to react in an isotropic way, generating the same pressure in all directions. However, Scholtès et al. reported on their recent studies [START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF][START_REF] Scholtès | Micromechanics of granular materials with capillary eects[END_REF] that this assumption is not valid. These authors performed numerical simulations of the triaxial compression of unsaturated granular materials. They assumed 1.2. Experimental investigations 15 that the components of the stress tensor includes two terms: the contributions of contact forces σ c ij (eective stress), and an isotropic stress due to capillary interactions σ cap ij , assuming an isotropic distribution of liquid bridges in the material. The additive eect of capillary forces in stress tensor is veried for an initial isotropic conguration, although deformation creates a slight anisotropy in liquid bridge distribution. Therefore, the contribution of capillary forces to the total stress is not an isotropic pressure, which is not consistent with the Bishop form of the eective stress, or any expression ignoring the deviatoric eect of the forces in liquid bridges.

Moreover, at the microscopic level, the inuence of liquid bridges on the tensile strength of unsaturated granular materials was observed in the pendular state by Kim

& Hwang [START_REF] Kim | Modeling of tensile strength on moist granular earth material at low water content[END_REF] and Gröger et al. [START_REF] Gröger | Modelling and measuring of cohesion in wet granular materials[END_REF]. They numerically studied the tensile strength under isotropic pressure. Kim and Hwang [START_REF] Kim | Modeling of tensile strength on moist granular earth material at low water content[END_REF] measured the actual magnitude of tensile strength induced by water in moist granular soils with low water contents (w < 4%).

They found that the magnitudes of the measured tensile strength are signicantly different from zero. The tensile strength generally increases with increasing water content and relative density. They also proposed a model to estimate approximately the tensile strength of an unsaturated granular material. Gröger et al. [START_REF] Gröger | Modelling and measuring of cohesion in wet granular materials[END_REF] reported that there was a large inuence of the surface roughness on the tensile stresses in wet particle systems.

By means of simulated tensile tests, they have shown that even when the tensile strength is reached, elastic body contacts supporting the external tensile load still exist. They also simulated shear tests to compare with tensile test simulations and concluded that for a small tensile load, the yield locus appeared to be the straight extension of the graphs for positive loads [START_REF] Gröger | Modelling and measuring of cohesion in wet granular materials[END_REF].

In addition, Gilabert et al. [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF] reported the inuences of assembling procedure of model cohesive powders under low pressure. They also studied the plastic compression, structural changes of this materials under isotropic loads with the same value of tensile strength between grains [START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF]. Under low pressure, the tensile strength plays a key role to stabilize the loose structures. Moreover, in recent study of Delenne et al. [START_REF] Delenne | Liquid clustering and capillary pressure in granular media[END_REF], they numerically investigated the process of growth and coalescence of liquid clusters in a granular material as the amount of liquid increases under the isotropic compaction. The most important nding of their work is that a peak grain pressure induced by capillary forces in a granular packing occurs inside the funicular state. It presents the transition from a primary coalescence process, where the volume of the largest cluster remains small, to a secondary coalescence process governed by the increase of liquid cluster volumes carrying a larger capillary stress [START_REF] Delenne | Liquid clustering and capillary pressure in granular media[END_REF]. A Fluorescence Microscopy is an optical microscope that uses uorescence and phosphorescence instead of, or in addition to, reection and absorption to study properties of organic or inorganic substances. The FM refers to any microscope that uses uorescence to generate an image, whether it is a simpler set up like an epiuorescence microscope, or a more complicated design such as a confocal microscope, which uses optical sectioning to get better resolution of the uorescent image. However, in the FM, uorophores lose their ability to uoresce as they are illuminated in a process called photobleaching.

Experimental investigations

Photobleaching occurs as the uorescent molecules accumulate chemical damage from the electrons excited during uorescence. Photobleaching can severely limit the time over which a sample can be observed by uorescent microscopy. Besides, the specimen must be uorescent to be suitable for uorescence microscopy and FM is also limited to 2D observations.

One of the newer and most promising qualitative methods for studying, and, where possible, quantifying the arrangements of aggregations/particles and voids in unsaturated soils is the ESEM. ESEM is a special type of scanning electron microscope that works under controlled environmental conditions and requires no conductive coating on the specimen. This makes it possible to examine wet samples and to preserve their natural characteristics for further testing, which is an obvious advantage of ESEM compared to the conventional SEM. This method enables to observe the surface of moist materials without pre-treatment, at dimensions from millimeters to micrometers. Nevertheless, when saturated samples are observed, the presence of a water lm covering the solid grains makes image analysis dicult to carry out, so that little information about the sample can be collected.

Getting more and more success over recent years, XRCT is a nondestructive technique for visualizing interior features within solid objects, and for obtaining digital information on their three-dimensional geometries and properties. A XRCT image is typically called a slice, as it corresponds to what the object being scanned would look like if it were sliced open along a plane. An even better analogy is a slice from a loaf of bread, because just as a slice of bread has a thickness, a XRCT slice corresponds to a certain thickness of the object being scanned. So, while a typical digital image is composed of pixels (picture elements), a XRCT slice image is composed of voxels (volume elements). Taking the analogy one step further, just as a loaf of bread can be reconstituted by stacking all of its slices, a complete volumetric representation of an object is obtained by acquiring a contiguous set of XRCT slices.

Experimental investigations 21

The gray levels in a XRCT slice image correspond to X-ray attenuation, which reects the proportion of X-rays scattered or absorbed as they pass through each voxel. X-ray attenuation is primarily a function of X-ray energy and the density and composition of the material being imaged. The great majority of XRCT systems use X-ray tubes, although tomography can also be done using a synchrotron or gamma-ray emitter as a monochromatic X-ray source.

X-ray source

Important tube characteristics are the target material and peak X-ray energy, which determine the X-ray spectrum that is generated; current, which determines X-ray intensity;

and the focal spot size, which impacts the spatial resolution of the tomography.

Most XRCT detectors utilize scintillators. Important parameters are scintillator material, size and geometry, and the means by which scintillation events are detected and counted. In general, smaller detectors provide better image resolution, but reduced count rates because of their reduced area compared to larger ones. To compensate, longer acquisition times are used to reduce the signal to noise ratio. Common scintillation materials are cesium iodide, gadolinium oxysulde, and sodium metatungstate.

The diagram in Fig. 1.4 illustrates the cone-beam conguration for laboratory XRCT scanners. In cone-beam scanning, the linear array is replaced by a planar detector, and the beam is no longer collimated. Data for an entire object, or a considerable thickness of it, can be acquired in a single rotation. The data are reconstructed into images using a cone-beam algorithm. Besides, planar-beam and parallel-beam congurations are also Chapter 1. Literature review the most common congurations for CT scanners. Further details of the equipment and applications can be found in Banhart [START_REF] Banhart | Advanced Tomographic Methods in Materials Research and Engineering[END_REF], and Hsied [START_REF] Hsied | Computed Tomography -Principles, Design, Artifacts, and Recent Advances[END_REF].

However, as the resolution is limited to about 1000-2000 voxels along the object crosssectional diameter, high resolution requires small objects. Finite resolution causes some blurring of material boundaries. The calibration of gray levels to attenuation coecients is complicated by the presence of polychromatic X-rays. For the large (decimeter) scale, geological specimens cannot be penetrated by low-energy X-rays, reducing resolving capability. Not all features have suciently large attenuation contrasts for useful imaging (e.g. carbonate fossils in carbonate matrix, quartz versus plagioclase. . . ). Image artefacts (beam hardening) can complicate data acquisition and interpretation. Large data volumes (in the order of a few gigabytes) can require considerable computer resources for visualization and analysis.

Typical observations

One of the most popular and basic method to observe the plasticity of granular materials is the standard one-dimensional compression test. Here, an example of the onedimensional compression test, performed in 1990 by Golightly [START_REF] Mcdowell | On the Yielding and Plastic Compression of Sand[END_REF], to analyze the yielding and plastic hardening of carbonate and silica sands, is shown in Fig. 1.5. At low stresses (region 1), the behavior is quasi-elastic, beyond yield (region 2), an approximately linear normal compression line emerges.

We present here typical studies using the XRCT to investigate the mechanical properties of wet granular material, from macroscopic to microscopic level.

Several authors reported on the liquid capillary bridges between grains and their eects on the mechanical behavior of wet granular materials [START_REF] Mason | Critical angle of wet sandpiles[END_REF][START_REF] Kohonen | On capillary bridges in wet granular materials[END_REF][START_REF] Fournier | Mechanical properties of wet granular materials[END_REF]. The specimens were prepared in the pendular state (with very small water content, < 1 %) by using the random loose/ dense packing method. An example of the liquid bridges between grains was clearly observed in Fig. 1.6 by using the Fluorescence Microscope by Fournier et al. [START_REF] Fournier | Mechanical properties of wet granular materials[END_REF]. The average number of bridges per bead was then obtained by zooming with the microscope through the sample and counting the bridges. Besides, they also described experiments on the dynamic deformation of a wet granular material, which reveal a type of shear-thinning behavior that they attributed to the shear-rate dependence of liquid bridge volumes.

More recently, ESEM has been used to observe wet granular materials [START_REF] Stokes | In situ mechanical testing of dry and hydrated breadcrumb in the environmental scanning electron microscope (esem)[END_REF][START_REF] Maeda | Evaporation and instabilities of microscopic capillary bridges[END_REF][START_REF] Musso | The role of structure in the chemically induced deformations of febex bentonite[END_REF][START_REF] Romero | Microstructure investigation in unsaturated soils: A review with special attention to contribution of mercury intrusion porosimetry and environmental scanning electron microscopy[END_REF][START_REF] Lourenço | Formation and evolution of water menisci in unsaturated granular media[END_REF]. This method allows observing water inside the materials. An example [START_REF] Lourenço | Formation and evolution of water menisci in unsaturated granular media[END_REF] of the use of this method to show clearly the water menisci at contacts of sand grains, is illustrated in Fig. 1.7. Fig. 1.7(a) shows predominantly concave menisci (corresponding to tensile capillary pressure) and Fig. 1.7(b) shows predominantly convex menisci (corresponding to compressive capillary water pressure). Lourenço et al. [START_REF] Lourenço | Formation and evolution of water menisci in unsaturated granular media[END_REF] observed the curvature of water menisci in both idealized and natural granular materials with vari- Maeda et al. [START_REF] Maeda | Evaporation and instabilities of microscopic capillary bridges[END_REF] also studied the formation and disappearance of liquid bridges between two surfaces of grains. They reported the slow evaporation and rupture of liquid bridges relatively to low vapor pressure. They also studied the mechanical (as opposed to the thermodynamic) stability of liquid bridges. These results shown clearly the importance of adhesion between surfaces in the presence of liquid bridges. For example, the rupture distance of such bridges is an important parameter in describing the dynamic mechanical properties of moist granular materials. This problem was reported in Ref. [START_REF] Lian | A Theoretical Study of the Liquid Brid Forces between Two Rigid Spherical Bodies[END_REF] and then conrmed by Seville et al. [START_REF] Seville | Interparticle forces in uidisation: A review[END_REF].

Early study performed by Crestana et al. [START_REF] Crestana | Static and dynamic three-dimensional studies of water in soil using computed tomographic scanning[END_REF], in 1985, used the XRCT to measure the water content of soil. They also showed that the XRCT can be applied to measure and follow dynamically the motion of water in soil in three dimensions. An example using the third-generation of the XRCT for spatial and real-time (dynamic) measurement is shown in Fig. sand [START_REF] Bruchon | Analyse par microtomographie aux rayons X de l'eondrement capillaire dans les matériaux granulaires[END_REF]. Bruchon et al. [START_REF] Bruchon | Full 3d investigation and characterisation of capillary collapse of a loose unsaturated sand using x-ray ct[END_REF] characterized in 3D the capillary collapse phenomenon during water inltration into a loose unsaturated sand. They used the volumetric digital image correlation tool to assess the volume of water-lled pores and deformation of the granular skeleton. Wildenschild et al. [START_REF] Wildenschild | Using x-ray computed tomography in hydrology: Systems, resolutions, and limitations[END_REF] also used this technique to characterize phase distribution and pore geometry for partially saturated media. They obtained qualitative and quantitative results with dierent scanning systems and sample sizes. Furthermore, Wildenschild and Sheppard [START_REF] Wildenschild | X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems[END_REF] used this technique for quantifying the pore-scale structure Chapter 1. Literature review and processes in subsurface porous medium systems. A pore-network reconstruction is shown in Fig. 1.10, illustrating the results of applying watershed partitioning and region merging on the pore and grain space of an unconsolidated sand. However, this method reveals a major weakness which is the use of spherical structure elements that breaks down when the grain shapes are very oblate or prolate [START_REF] Wildenschild | X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems[END_REF]. Farber et al. [START_REF] Farber | Use of x-ray tomography to study the porosity and morphology of granules[END_REF] studied the porosity, pore size distribution, and geometric structure of pores in granulates by using dierent conditions and materials. In addition, the mor-phological information such as pore shape, spatial distribution, and connectivity were also observed by combining with mercury porosimetry. Kim et al. [START_REF] Kim | Water distribution variation in partially saturated granular materials using neutron imaging[END_REF] reported the effect of particle morphology on water distribution in partially saturated granular porous media. They also investigated the inuence of water content in order to evaluate the water phase-distribution spatially for compacted sand specimens in a aluminum cylinder.

Riedel et al. [START_REF] Riedel | Water retention behaviour explored by x-ray ct analysis[END_REF] characterized experimentally the link between partial water saturation and suction in a sand sample at the micro scale. Four dierent levels of imposed suction of a sample (with Hostun sand) obtained from trinarised 3D images, as shown in Fig. 1.11 [START_REF] Riedel | Water retention behaviour explored by x-ray ct analysis[END_REF]. The vertical location of the slices is chosen to best highlight dierent water retention states. The two slices on the top row show liquid phase continuity for both saturations. At higher suction, the pendular domain is reached in some parts of the sample (in bottom left of Fig. In a more recent study, Menon et al. [START_REF] Menon | Analysing the impact of compaction of soil aggregates using X-ray microtomography and water ow simulations[END_REF] used XRCT combined with numerical sim-Chapter 1. Literature review ulation, to quantify the impact of compaction on macroaggregates, mainly on the pore size distribution and water ow. The impact of compaction on interpores and intrapores is presented in the proportion of inter-and intrapores, as shown in Fig. 1.12. Note that interpores or inter-aggregate pores, are the pores between soil aggregates; intrapores or intra-aggregates pores, are within soil aggregates (within the solid matrix of soil aggregates). Interpores dominated the total pores volume in comparison to the intrapores, representing > 90 % of the total pore volume before compaction in pre-compacted samples. However, after compaction, there was an increase in intrapores proportion in all cases. In addition, the total surface area of interpores and intrapores were also calculated. They concluded that in all cases, the decrease in interpores proportion produced a corresponding increase in intrapores fraction [START_REF] Menon | Analysing the impact of compaction of soil aggregates using X-ray microtomography and water ow simulations[END_REF]. Moreover, the size distribution of interpores volumes (i.e. volume of individual interpore) before and after compaction was also reported. They showed that the soil compaction and moisture content increased the total number of individual (very small) interpores (< 0.001 mm 3 ) although the total volume of interpores decreased sharply [START_REF] Menon | Analysing the impact of compaction of soil aggregates using X-ray microtomography and water ow simulations[END_REF]. methods have been widely used to investigate the macroscopic properties as well as the internal states of granular systems. In other words, these methods allow to investigate the behavior of granular systems from the macroscale to the microscale (and even the nanoscale); this is hardly made possible by using experimental methods.

The most common computer simulation method used to determine uid ow and analyze transport properties through porous media is the LBM [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF][START_REF] Sukop | Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers[END_REF]. LBM may be regarded as a digital method equivalent to the traditional Computational Fluid Dynamics.

Several authors used the LBM to investigate uid ow and permeability in porous materials [START_REF] Coles | Pore level imaging of uid transport using synchrotron X-ray microtomography[END_REF][START_REF] Harting | Large-scale grid-enabled lattice Boltzmann simulations of complex uid ow in porous media and under shear[END_REF][START_REF] Selomulya | Direct Prediction of Structure and Permeability of Flocculated Structures and Sediments Using 3D Tomographic Imaging[END_REF][START_REF] Wang | Domain-decomposition method for parallel lattice Boltzmann simulation of incompressible ow in porous media[END_REF][START_REF] Selomulya | An integrated methodology to evaluate permeability from measured microstructures[END_REF][START_REF] Sukop | Distribution of multiphase uids in porous media: Comparison between lattice Boltzmann modeling and micro-x-ray tomography[END_REF]. The permeability and tortuosity of porous materials were also studied by the RWS [START_REF] Nakashima | Estimate of transport properties of porous media by microfocus X-ray computed tomography and random walk simulation[END_REF][START_REF] Nakashima | Threedimensional diusion of non-sorbing species in porous sandstone: Computer simulation based on X-ray microtomography using synchrotron radiation[END_REF]. SPH is a meshless Lagrangian computational technique that was introduced in the astrophysics community to simulate the movement of masses of material in an unbounded three-dimensional space [START_REF] Potapov | Liquid-solid ows using smoothed particle hydrodynamics and the discrete element method[END_REF][START_REF] Cleary | Discrete-element modelling and smoothed particle hydrodynamics: potential in the environmental sciences[END_REF][START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF][START_REF] Monaghan | Smoothed Particle Hydrodynamics and Its Diverse Applications[END_REF][START_REF] Chen | Numerical simulations for large deformation of granular materials ssing Smoothed Particle Hydrodynamics Method[END_REF].

SPH is used to obtain approximate numerical solutions of the equations of uid dynamics by replacing the uid by a set of particles. MPM, known as the particle-in-cell method, combines the strengths of Eulerian and Lagrangian descriptions of the material. The Lagrangian description is provided by discretizing each body by a collection of material points, and the Eulerian description is based on a background computational mesh. The MPM extends these capabilities to materials modeling [START_REF] Bardenhagen | The material-point method for granular materials[END_REF][START_REF] Bardenhagen | The generalized interpolation material point method[END_REF][START_REF] Guilkey | Implicit time integration for the material point method: Quantitative and algorithmic comparisons with the nite element method[END_REF][START_REF] Eckowski | The material point method in large strain engineering problems[END_REF]. PFEM uses a Lagrangian description to model the motion of nodes (particles) in both the uid and the structure domains. Nodes are thus viewed as particles which can freely move and even separate from the main analysis domain representing, for instance, the eect of water drops. A mesh connects the nodes dening the discretized domain where the governing equations, expressed in an integral form, are solved as in the standard Finite Element Method (FEM) [START_REF] Oñate | The Particle Finite Element Method -an Overview[END_REF][START_REF] Idelsohn | The particle nite element method: A powerful tool to solve incompressible ows with free-surfaces and breaking waves[END_REF].

The most popular computer simulation methods used to calculate the contact and the displacement of particles in granular materials is the discrete element type. This method was rstly introduced in 1979 [START_REF] Cundall | A disrete numerical model for granular assemblies[END_REF] and has proved to be an ecient and widespread tool to observe the micro-structure and micro-mechanical behavior of granular systems.

Two versions of the molecular dynamics simulation method are introduced. The rst is the so-called soft sphere molecular dynamic (MD=DEM). It is a straightforward implementation to solve the equations of motion for a system of many interacting particles [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF]. For DEM, both normal and tangential interactions, like friction, are discussed for spherical spheres. The integration of the dynamic equation is illustrated in Fig. 1.13. The second method is the so-called Event-Driven (ED) or Contact Dynamics (see Jean in [START_REF] Radjaï | Discrete-element Modeling of Granular Materials[END_REF]) simulation, which is conceptually dierent from DEM, since Chapter 1. Literature review collisions are dealt with using a collision matrix that determines the momentum change on physical grounds. For the sake of brevity, the ED method (see McNamara in [START_REF] Radjaï | Discrete-element Modeling of Granular Materials[END_REF]) is only discussed for smooth spherical particles. In addition, the quasi-static method is developed based on the stiness matrices of the granular assemblages. This method is appropriate to analyze the granular assemblages at the equilibrium or in the quasi-static deformation (see Roux and Combe in [START_REF] Radjaï | Discrete-element Modeling of Granular Materials[END_REF]). In some cases, a computational method coupling two methods is necessarily applied to study simultaneously two objects, for instance the uid and particles ows. Hence, a coupling of the SPH and the DEM was proposed in order to simulate the liquid-solid ows in granular materials [START_REF] Potapov | Liquid-solid ows using smoothed particle hydrodynamics and the discrete element method[END_REF][START_REF] Cleary | Discrete-element modelling and smoothed particle hydrodynamics: potential in the environmental sciences[END_REF]. Another coupling of two methods, used in mechanics of granular media and uid dynamics, are the LBM and the DEM, respectively [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent uid ows: Computational issues[END_REF][START_REF] Han | Coupled lattice Boltzmann and discrete element modelling of uid-particle interaction problems[END_REF][START_REF] Feng | Combined three-dimensional lattice Boltzmann method and discrete element method for modelling uid-particle interactions with experimental assessment[END_REF][START_REF] Lominé | Modeling of uid-solid interaction in granular media with coupled lattice Boltzmann / discrete element methods: application to piping erosion[END_REF][START_REF] Galindo-Torres | A coupled Discrete Element Lattice Boltzmann Method for the simulation of uid-solid interaction with particles of general shapes[END_REF][START_REF] Aminpour | Micro-scale characterization of uid ow in a uniform particle pack using coupled discrete element-lattice Boltzmann method[END_REF]. Besides, several authors proposed a coupling of the DEM and the FEM to tackle the behavior of granular assemblies from the microscale to the macroscale.

The concept of Representative Elementary Volume (REV) is applied to bridge the microand macro-scales. This computational homogenization (FEM × DEM) has been use to simulate biaxial tests [START_REF] Nitka | Two-scale modeling of granular materials: a DEM-FEM approach[END_REF][START_REF] Nguyen | FEM x DEM modelling of cohesive granular materials : numerical homogenisation and multi-scale simulation[END_REF][START_REF] Kh | Numerical simulation of particle breakage of angular particles using combined DEM and FEM[END_REF], triaxial tests [START_REF] Wellmann | A two-scale model of granular materials[END_REF], and for the dynamic analysis of geomechanics problems [START_REF] Oñate | Combination of discrete element and nite element methods for dynamic analysis of geomechanics problems[END_REF]. In addition, a coupling between the MPM and the FEM was implemented to demonstrate a good agreement between explicit and implicit solutions to the problem of stress wave propagation [START_REF] Guilkey | Implicit time integration for the material point method: Quantitative and algorithmic comparisons with the nite element method[END_REF].

Typical observations

The DEM or MD simulation are getting more and more popular in granular material science, for 2D and 3D applications. However, studies of cohesive granular materials, especially in loose states, have been investigated less frequently than dense or noncohesive granular materials. For instance, Onoda and Liniger [START_REF] Onoda | Random loose packings of uniform spheres and the dilatancy onset[END_REF] studied the randomloose-packing fraction of uniform spheres at the limit of zero gravitational force. Dong et al. [START_REF] Dong | Role of interparticle forces in the formation of random loose packing[END_REF] then reported on the loose structures formed by particle packed under gravity and stabilized by adhesion, as illustrated in Fig. 1.14. It can be seen from Fig. 1.14 that the DEM can reasonably reproduce the experimental results for all sized particles, conrming the validity of the proposed simulation technique [START_REF] Dong | Role of interparticle forces in the formation of random loose packing[END_REF]. square, d = 250 µm; lozenge, d = 110 µm; and ×, result of [START_REF] Onoda | Random loose packings of uniform spheres and the dilatancy onset[END_REF]. Lines are the simulated results (from Dong et al. [START_REF] Dong | Role of interparticle forces in the formation of random loose packing[END_REF]).

Other authors investigated loose packing stabilized by cohesion and the collapse of such materials under growing loads during compression tests [START_REF] Kadau | Contact dynamics simulations of compacting cohesive granular systems[END_REF][START_REF] Kadau | Pore stabilization in cohesive granular systems[END_REF][START_REF] Wolf | Compaction of cohesive powders[END_REF][START_REF] Bartels | The eect of contact torques on porosity of cohesive powders[END_REF]. These authors studied a dynamical compression regime, and observed a shock wave propagating through the samples. Fig. 1.15 [START_REF] Kadau | Pore stabilization in cohesive granular systems[END_REF] shows a simulation of a ballistic deposit containing 2746 monodisperse spherical particles, the top gure is the initial arrangement, the bottom gure is the nal conguration. The contact between particles includes Coulomb friction, rolling friction and cohesion. They also compared with the case using only

Coulomb and rolling friction. This comparison conrms that the eect of cohesion on the porosity is higher if one uses Coulomb and rolling friction. Moreover, the low cohesive forces lead to similar packings as without cohesion [START_REF] Kadau | Pore stabilization in cohesive granular systems[END_REF].

Gilabert et al. performed a 2D investigation into the microstructure and mechanical

properties of model cohesive powders [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF], and also the quasi-static behavior of the cohesive powders under an isotropic compression test [START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF]. Fig. 1.16 [START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF] shows how this initial inuence aects the beginning of compression curves and fades out later on. (from Kadau et al. [START_REF] Kadau | Pore stabilization in cohesive granular systems[END_REF]).

Compression curves are shown in

low pressure P * , as does a decrease of agitation level (V 0 /V * ). However, such a change in material properties aects not only the initial stage (rst plateau of compaction curve) but also the macroscopic mechanical behavior at larger densities. As a consequence, the slope of the compaction curve is lower for larger rolling resistances [START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF].

Figure 1.16: Compression and decompression curves: eect of initial agitation level in aggregation stage and inuence of rolling resistance parameter (from Gilabert et al. [START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF]).

In addition, several studies focused on the fracture of bound particle assemblies in static [START_REF] Delenne | Mechanical behaviour and rupture of cohesive granular media[END_REF][START_REF] Delenne | Mechanical behaviour and failure of cohesive granular materials[END_REF] or dynamic conditions [START_REF] Thornton | How do agglomerates break?[END_REF]. The others studied wet beads packs in 1.3. Numerical investigations which neighboring particles are joined by liquid bridges, and investigated the structure of poured samples [START_REF] Yang | Numerical study of the packing of wet coarse uniform spheres[END_REF], or shear strength of this material [START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF]. Besides, the shear ows of cohesive granular materials have been also simulated [START_REF] Brewster | Plug ow and the breakdown of bagnold scaling in cohesive granular ows[END_REF][START_REF] Aarons | Shear ow of assemblies of cohesive and non-cohesive granular materials[END_REF][START_REF] Rognon | Rheophysics of cohesive granular materials[END_REF][START_REF] Luding | Shear ow modeling of cohesive and frictional ne powder[END_REF][START_REF] Luding | The eect of friction on wide shear bands[END_REF]. This means that with the DEM, it is easily to simulate all the tests with any material.

In the microscopic point of view, the DEM simulation give easily access to the microstructure of wet granular materials, thus allowing the measurement of intergranular forces between grains, contact networks, etc. For instance, Richefeu et al. [START_REF] Richefeu | Stress transmission in wet granular materials[END_REF] analyzed stress transmission in wet granular media in the pendular state by means of 3D molecular dynamics simulations. Fig. 1.17 shows the force network in a narrow slice whose thickness is nearly three particle diameters. The tensile and compressive forces are represented by segments of dierent colors joining particle centers. The line thickness is proportional to the force. It is remarkable that tensile and compressive force chains could be obviously observed although the slice is quite narrow [START_REF] Richefeu | Stress transmission in wet granular materials[END_REF]. Besides, the distribution of forces and the connectivity of the bond network are the main properties, that most often are characterized in studies using DEM simulations. The appearance of a small quantity of water in the intergranular voids, will preferentially accumulate near the contacts or the narrow interstices, forming liquid bridges that join contacting grains or near neighbors. These bridges are formed by assuming the shape of the meniscus as illustrated in Fig. 1.18. This formation may be spontaneous if the liquid has exchanges with its vapor, present in the surrounding atmosphere. The pressure dierence between two sides of a curved liquid-air or liquid-solid interface is determined by the Laplace-Young equation which is the product of the surface tension Γ and principal radii of curvature r 1 and r 2 (as shown in Fig. 1.18(a)):

∆u = u a -u w = Γ 1 r 2 - 1 r 1 , (1.4)
in which u a -u w is the pressure dierence across the uid interface and Γ is the surface tension of the air-liquid interface.

(a) The capillary force F cap due to the liquid bridge can be also calculated according to the gorge method [START_REF] Hotta | The capillary binding force of a liquid bridge[END_REF]. This force involves eects of both the capillary pressure ∆u and the surface tension Γ over the bridge section, as

R 1 R 2 q j r 1 r 2 h (b) R q j r 2h 2a
F G cap = πr 2 2 ∆u + 2πr 2 Γ = πΓr 2 1 + r 2 r 1 .
(1.5)

Numerical investigations 35

For the monodisperse case, r 1 and r 2 , respectively denote the radius of the meridian prole and the radius at the neck (Fig. 1.18(b)), are given by

     r 1 = h/2 + R(1 -cos ϕ) cos(ϕ + θ) , r 2 = R sin ϕ + [1 -sin(ϕ + θ)]r 1 .
(1.6)

The lling angle ϕ cannot be calculated explicitly and an iterative procedure must be used to get it as a function of the other parameters. The contact angle θ is illustrated in the relationship between capillary force, suction, and the gorge's radius (see Fig. 1.

19).

The contact angle is null when the gorge's radius is larger than 0.01 to 0.1 mm, the capillary force varies from 78 µN to 57 µN, and the suction from 23 kPa to 0.3 kPa.

Therefore, the strong variation of suction slightly aects the capillary force. Moreover, when the contact angle increases, the capillary force decreases at the constant of gorge's radius, and the capillary force reaches a maximum value of 2πΓR cos θ. In experiments on granular materials, it is really dicult to determine this contact angle. Hence, for the simplication, the contact angle is usually considered null in many dierent models.

Capillary force (µN)

The gorge's radius (µm) Suction (kPa) Besides, Soulié et al. [START_REF] Soulié | Inuence of liquid bridges on the mechanical behaviour of polydisperse granular materials[END_REF] also proposed a numerical solution to calculate the capillary force, as below

F cap = πΓ R 1 R 2 c + exp a h R 2 + b , (1.7)
where R 1 and R 2 are the radii of particles and h is also the interparticle distance. The coecients a, b, c are functions of the volume of liquid bridge V m , the contact angle θ, Chapter 1. Literature review and R 2 :

               a = -1.1 V m R 3 2 -0.53 b = -0.148 ln V m R 3 2 -0.96 θ 2 -0.0082 ln V m R 3 2 + 0.48 c = 0.0018 ln V m R 3 2 + 0.078 (1.8)
In the case of a capillary bridge at the contact point between two equally-sized spheres with radius R (see Fig. 1.18(b)), assuming that the principal radii of curvature are denoted a and r, relation (1.4) can be rewritten as:

∆u = Γ 1 a - 1 r .
(1.9)

As the diameter a of the meniscus is much larger than its meridian radius of curvature r, ∆u is negative. The liquid is sucked towards the contact region. The angle θ is the equilibrium contact angle between the liquid and the grain surfaces. In the case of perfect wettable materials, this contact angle is equal to zero.

The capillary force between equal-sized spheres is approximately equal to the suction force πa 2 ∆u. If we assume a r, the pressure dierence is proportional to -Γ/r. For the very small grain size, assuming h R (2h is the distance between grain surfaces),

we nally obtain the capillary force, F cap 2πΓR cos θ.

(1.10)

The capillary force between two equal-sized spheres in contact is proportional to surface tension Γ and particle radius R. This relation is remarkably independent of the volume of meniscus V m . When the particles are small enough, the inuence of gravity on the curvature of the liquid bridge and on the capillary force is negligible [START_REF] Santamarina | Soil Behavior at the Microscale: Particle Forces[END_REF].

In order to calculate approximately the capillary force, Maugis [START_REF] Maugis | Contact, Adhesion and Rupture of Elastic Solids[END_REF] proposed a simpler formula but less accurate, also called Maugis approximation, by which the capillary force is expressed as a function of the meniscus volume V m and the interparticle distance h (also called separation distance) as below,

F cap = 2πΓR cos θ 1 - 1 1 + 2V m πRh 2 . (1.11)
This formula is chosen to calculate the capillary force of the monodisperse grains in the present works.

Fig. 1.20 shows a typical example of the evolution of the capillary force with the relative displacement between two interacting grains for a given capillary pressure. The capillary force acts exclusively in the axial direction of the interaction (the normal to the contact plane). F cap is maximum and remains constant for grains strictly in contact (h ≤ 0), assuming small deformations of the contact zone compared with the grain size. D r (rupture distance) is the debonding distance dened as the maximum distance h beyond which the Laplace-Young equation has so solution. A hydraulic hysteresis can be introduced by dening D c (creation distance) as the distance below which a liquid bridge can form. For numerical simulation, a possible choice is to allow a bridge to form when particles come strictly in contact D c = 0, hence neglecting the possible eect of adsorbed water on the grain surface. This assumption is used in our study. The rupture distance D r can be found in several studies existing. Lian et al. [START_REF] Lian | A Theoretical Study of the Liquid Brid Forces between Two Rigid Spherical Bodies[END_REF] proposed the following relation between the rupture distance D r , the liquid bridge volume V m and the contact angle θ:

D r = 1 + θ 2 V 1/3 m .
(1.12)

The formula proposed by Herminghaus [START_REF] Herminghaus | Dynamics of wet granular matter[END_REF] also related to meniscus volume V m and contact angle θ, as follows:

D r = 1 + θ 2 V 1/3 m + 0.1V 2/3 m . (1.13)
The approximate formula of rupture distance, D r V 1/3 m , was also proposed by Willett et al. [START_REF] Willett | Capillary bridges between two spherical bodies[END_REF], and Maeda et al. [START_REF] Maeda | Evaporation and instabilities of microscopic capillary bridges[END_REF]. This approximate equation has been widely used in the numerical simulations and it is applied in the present study. An experimental study by Pitois et al. [START_REF] Pitois | Liquid Bridge between Two Moving Spheres: An Experimental Study of Viscosity Eects[END_REF] shows that the rupture distance also varies with the separation velocity.

Loose conguration of adhesive particles

Adhesive particles tend to form tenuous aggregates, which may sustain the increasing conning stress without rearranging, as shown for a 2D model [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF] in Fig. 1.21. The authors observed that the equilibrated congurations of cohesive packing are sensitive to the applied pressure. The structures can be stabilized for very small applied pressure, but then collapse when increasing the pressure in some pressure range. This can be found in the rst regime of compression curves in 2D models of loose systems [START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF], and in 3D ones [START_REF] Than | Basic mechanical properties of wet granular materials: A DEM study[END_REF], as shown in Fig. 1.22. Indeed, the rst part of the compression curves (see rectangle boxes) conrms that the initial structure is able to sustain the increasing pressure without rearranging. Note that the compression curve of Gilabert et al. [START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF] is expressed as the relationship between solid fraction Φ and reduced pressure P * instead of the e -log(P * ) curve conventionally used in soil mechanics. Void ratio, e (-)

Reduced pressure, P * (-)

Experimental investigation versus numerical investigation 1.4.1 Existing studies combining experimental and numerical methods

It is well recognized that experimental methods hardly enable to get insight into the microstructure of specimens, especially to measure precisely intergranular forces between particles subjected to growing external forces. However, with the advanced techniques, we can count accurately the capillary bridges between particles (using FM for instance), or capture and distinguish clearly convex and concave menisci (using ESEM for instance), or access and reconstruct completely the microstructure of specimens at dimensions from millimeters to micrometers (using XRCT for instance). In practice, some factors can considerably and directly aect the results. Among other factors, one can cite the technique itself, the quality and the inherent noise of the equipment, or the environment, which is of critical importance in the case of in-situ observations. Therefore, we sometimes need to perform several tests for measuring a given quantity. Generally, it is really dicult to overcome the inconvenience of external factors. Conversely, with the numerical simulations, we can idealize the materials and the testing conditions, with minimal disturbance from external factors. Alternatively, we can simulate in a realistic manner the imperfect conditions applied or observed while testing physically the specimens. Therefore, the numerical results will be helpful to get information on unmeasurable quantities or to interpret non-ideal experiments. Each method, either numerical or experimental investigations, has its own strengths as well as inherent limitations. Both could provide useful observations to improve the general knowledge. As a consequence, the combination of numerical and experimental observations is not only a way to overcome the shortcomings of each method but also 40 Chapter 1. Literature review an opportunity to complete and validate each other. A summary of the main applications that have utilized the combination of experimental observations and computer simulations to study the behavior of porous or granular materials is provided in Table 1.3 (partially succeeding from Moreno-Atanasio et al. [START_REF] Moreno-Atanasio | Combining X-ray microtomography with computer simulation for analysis of granular and porous materials[END_REF]). Selected examples (typical studies) found in the literature are shown in Table 1.3, where are indicated the application, material and parameter used to link experiment and computer simulation, and appropriate referenced works.

The link between experiments and computer simulations is always established through the exchange of structural information between the two techniques. In Table 1.3, the link between experimental and numerical investigations mainly focuses on the following phenomena: uid ow and transport properties, mechanical loading of granules and granular beds, and particle packing. A brief summary of the most characteristic parameters obtained from XRCT and computer simulations is described below.

Firstly, the study on oil and water distribution within the beads, reported by Sukop et al. [START_REF] Sukop | Distribution of multiphase uids in porous media: Comparison between lattice Boltzmann modeling and micro-x-ray tomography[END_REF], has shown one of the most robust ways of coupling XRCT with computer simulations. In this study, the initial structure was transferred into the computer software.

The authors then compared the distribution of oil and water by comparing the relative content of each uid plane by plane. The distribution of the two phases within specic pores was then visualized. In addition, Selomulya et al. [START_REF] Selomulya | Direct Prediction of Structure and Permeability of Flocculated Structures and Sediments Using 3D Tomographic Imaging[END_REF] also presented the transport properties (permeability and tortuosity) of lter cakes of which XRCT digital images were transferred into a LBM software. The simulated relative permeabilities were calculated according to Darcy's law. Moreover, Wang et al. [START_REF] Wang | Domain-decomposition method for parallel lattice Boltzmann simulation of incompressible ow in porous media[END_REF], using the LBM simulations, obtained a similar value of permeability for exactly the same bead which was used in the experiments by Nakashima and Watanabe [START_REF] Nakashima | Estimate of transport properties of porous media by microfocus X-ray computed tomography and random walk simulation[END_REF].

Secondly, Golchert et al. [START_REF] Golchert | Application of X-ray microtomography to numerical simulations of agglomerate breakage by distinct element method[END_REF][START_REF] Golchert | Eect of granule morphology on breakage behaviour during compression[END_REF] studied the inuence of the structure of granules made of glass beads on the compressive behavior of the granule using XRCT and DEM.

They determined particle sizes and positions from the XRCT images and this information was used to create the computer simulated granule. Besides, Stock [START_REF] Stock | Recent advances in X-ray microtomography applied to materials[END_REF] suggested that the displacement of marker particles in powder beads allows the determination of the eld displacement. For instance, Fu et al. [START_REF] Fu | Application of x-ray microtomography and image processing to the investigation of a compacted granular system[END_REF] used marker particles in the study of the compressive behavior of a bead of sugar particles mixed with some glass particles. They observed that the markers at the bottom of the bead were displaced less than the markers near the top of the bead. The authors stated that their results were in agreement with computer simulation results obtained by other authors using FEM [START_REF] Wu | Modelling the mechanical behaviour of pharmaceutical powders during compaction[END_REF]. Therefore, the comparison was purely qualitative.

Finally, the combination of XRCT and computer simulations allows a direct and easy comparison and validation of the particle packing predictions by computer models [START_REF] Fu | Investigation of particle packing in model pharmaceutical powders using x-ray microtomography and discrete element method[END_REF][START_REF] Fu | Application of x-ray microtomography and image processing to the investigation of a compacted granular system[END_REF]. The investigated materials include the monodisperse glass spheres, polydisperse glass spheres, and irregular polydisperse microcrystalline cellulose particles. By obtaining dierent cross sections from the XRCT, the particle size distribution was ob-tained, compared to SEM results and provided to DEM simulations. They found a good agreement for the radial distribution function between experiment and simulation, except the case for cellulose particles. This work shows a typical case in which limitations in the computer models are the origin of the dierences between the data obtained from XRCT and the predicted simulation data. In contrast to Fu et al. [START_REF] Fu | Investigation of particle packing in model pharmaceutical powders using x-ray microtomography and discrete element method[END_REF][START_REF] Fu | Application of x-ray microtomography and image processing to the investigation of a compacted granular system[END_REF], Jia et al. [START_REF] Jia | Validation of a digital packing algorithm in predicting powder packing densities[END_REF][START_REF] Jia | A novel approach to predicting the behavior of arbitrary particulate mixtures under vibration[END_REF] studied the predicted packing of irregular particles by using the real particle shape as an input for the computer simulations. Therefore, the link between simulation and experiments presented by these authors is a combination of parametric and imaging link.

Furthermore, a sophisticated algorithm able to produce particulate beads with a specic value of coordination number and packing fraction distribution has been developed [START_REF] Al-Raoush | Simulation of random packing of polydisperse particles[END_REF].

The results were compared to the values obtained from XRCT for a bead made of glass beads [START_REF] Al-Raoush | Extraction of physically realistic pore network properties from three-dimensional synchrotron X-ray microtomography images of unconsolidated porous media systems[END_REF]. They also found a good agreement between simulations and experiments in all cases of monodisperse and polydisperse spheres. [START_REF] Al-Raoush | Simulation of random packing of polydisperse particles[END_REF][START_REF] Al-Raoush | Microstructure characterization of granular materials[END_REF] Bed structure, packing, bed com-pression Glass mono, glass poly and MCC DEM XRCT Input structure, pair correlation function Fu et al. [START_REF] Fu | Investigation of particle packing in model pharmaceutical powders using x-ray microtomography and discrete element method[END_REF][START_REF] Fu | Application of x-ray microtomography and image processing to the investigation of a compacted granular system[END_REF] Bed structure, packing Industrial powders Digital model XRCT Comparison of porosity and structure between real and simulated Jia et al. [START_REF] Jia | Validation of a digital packing algorithm in predicting powder packing densities[END_REF][START_REF] Jia | A novel approach to predicting the behavior of arbitrary particulate mixtures under vibration[END_REF] Granule compression Glass beads DEM XRCT Input structure Golchert et al. [START_REF] Golchert | Application of X-ray microtomography to numerical simulations of agglomerate breakage by distinct element method[END_REF][START_REF] Golchert | Eect of granule morphology on breakage behaviour during compression[END_REF] Stress-dilatancy and strength

Comparison

Feng et al. [START_REF] Feng | Combined three-dimensional lattice Boltzmann method and discrete element method for modelling uid-particle interactions with experimental assessment[END_REF] Note: EXP, Experiment; SIM, Simulation; MCC, MicroCrystalline Cellulose; CT, Computed Tomography 1.4.2 Combination between the XRCT and the DEM One of the most suitable simulation techniques to couple with the XRCT in the study of mechanical properties of granular materials is the DEM. DEM considers particles as individual entities with physical and geometrical properties (elastic modulus, density, particle size distribution and dierent shapes). The DEM simulation could provide a direct correlation between external load and internal granular behavior with realistic contact deformation [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF]. Therefore, we present briey here some of the most important parameters that are obtained by combining XRCT and DEM, such as particle shape, porosity distribution and structures of pores, characterization of contact network, radial distribution function, and particle displacements.

Particle shape

The most common way to simplify particle shape is to consider irregular particles as ellipsoids. Wang et al. [START_REF] Wang | Representation of real particles for DEM simulation using X-ray tomography[END_REF] applied this technique to study the irregular particles (limestone and sandstone), whose shapes were obtained from XRCT. The model particles were then characterized by having the same mass and inertia moments as the original particles. This method was described in terms of a spherical harmonic expansion [START_REF] Garboczi | Three-dimensional mathematical analysis of particle shape using X-ray tomography and spherical harmonics: Application to aggregates used in concrete[END_REF][START_REF] Taylor | Some properties of irregular 3-D particles[END_REF].

By applying the volume of the equivalent ellipsoids, Wang et al. [START_REF] Wang | Representation of real particles for DEM simulation using X-ray tomography[END_REF] shown that the ellipsoid volume overestimated the real particle volume meanwhile the reconstructed images of limestone and sandstone were reproduced qualitatively well by the equivalent ellipsoidal simulation. Therefore, this combination is extremely useful because it allows an easy manipulation of irregular particle without making the codes very complex. Several algorithms used to describe dierent particle shape properties (such as roundness and inertia moments) were also proposed by Al-Raoush [START_REF] Al-Raoush | Microstructure characterization of granular materials[END_REF] and Lin and Miller [START_REF] Lin | 3D characterization and analysis of particle shape using X-ray microtomography (XMT)[END_REF].

Wang et al. [START_REF] Wang | Representation of real particles for DEM simulation using X-ray tomography[END_REF] also proposed the clumping of spheres method, using XRCT images to simplify particle shapes, and they can be easily handled by simulations. Fig. 1.23 shows two schematic representations of the clumping of spheres method in which the tomographic image of a real particle was used to create a computational equivalent particle. Fig. 1.23(a) illustrates the equivalent particle using the non-overlapping technique.

The overlapping method of clumping the spheres is presented in Fig. 1.23(b). In the rst method, the shape and volume of the nal modeled particle are very dierent from the real one because of the presence of pores between the clumped spheres, unlike the real object. Therefore, the nal structure is made of clusters of spheres which do not show a relative motion between them, in principle [START_REF] Wang | Representation of real particles for DEM simulation using X-ray tomography[END_REF]. However, with the second method, it allows the overlap of the individual spheres, and therefore avoiding the interstitial articial pores. Nevertheless, this method may be computationally very expensive depending on the irregularity of the real particle and on the requirement of accuracy of the simulated particle [START_REF] Moreno-Atanasio | Combining X-ray microtomography with computer simulation for analysis of granular and porous materials[END_REF].

Other authors proposed a simple solution and ignored the non-sphericity of particles if the particle shape was not strongly elongated [START_REF] Golchert | Application of X-ray microtomography to numerical simulations of agglomerate breakage by distinct element method[END_REF][START_REF] Golchert | Eect of granule morphology on breakage behaviour during compression[END_REF]. In these cases, the centers and radii of particles can be obtained directly from XRCT and these parameters are straightforwardly transferred to computer simulation.

Porosity, porosity distribution and structures of pores

Porosity and distribution of porosity in a sample is probably the most common parameter used to link XRCT and DEM. This parameter strongly aects the mechanical and physical behavior of materials, such as mechanical strength, ow-ability, or permeability. Therefore, the porosity predicted by simulations necessarily has to be compared and validated with the values obtained from XRCT. Some typical works reported here have not been linked with simulations but these works provide the possibility to relate XRCT and numerical simulations in future works.

The analysis of the bulk or local porosity of the samples exposes some problems: the scanned areas may not be representative of the whole specimen, the identied algorithm for dierentiation between grains and pores (namely, the segmentation process) may be inaccurate, and the resolution of XRCT images may not be high enough. Therefore, we can consider dierent local scans from dierent parts of the specimen as a representative value of porosity. We use this value to compare with the bulk porosity of whole sample.

Selumulya et al. [START_REF] Selomulya | Direct Prediction of Structure and Permeability of Flocculated Structures and Sediments Using 3D Tomographic Imaging[END_REF] performed six scans of dierent sample region in order to obtain an appropriate value of bulk porosity of the material, and compared with the porosity of the real sample. The authors analyzed the porosity of the sample versus height. This is the most common way of characterizing the porosity of materials. Interesting studies regarding the analysis of porosity distribution with height within a cylindrical specimen has been performed [START_REF] Suzuki | Measurement of voidage distribution in particle packed bed using X-ray micro computed tomography[END_REF][START_REF] Than | Macromicroscopic one-dimensional compression of wet granular soils by experimental investigation[END_REF].

The characterization of the porosity of samples was also analyzed by the pore network structure in term of the number of pore bodies, pore throat, and the coordination number of pores [START_REF] Al-Raoush | Extraction of physically realistic pore network properties from three-dimensional synchrotron X-ray microtomography images of unconsolidated porous media systems[END_REF][START_REF] Taud | Porosity estimation method by X-ray computed tomography[END_REF]. Al-Raoush and Willson [START_REF] Al-Raoush | Extraction of physically realistic pore network properties from three-dimensional synchrotron X-ray microtomography images of unconsolidated porous media systems[END_REF] considered eight dierent systems mainly made of glass beads although some of them were made of natural marine sand. They found that the dierences in the number of pores, body pore sizes and pore coordination number is approximately 50 % despite the fact that the macroscopic value of porosity was similar in all samples. Interestingly, they also compared the results with low and high resolution XRCT images and concluded that a good resolution is needed to minimize the error of the estimation. This is one of the most detailed study in literature which analyzes clearly pore properties.

Other studies using XRCT, with no comparison with simulation, were performed to determine parameters such as the average pore diameter, number of pores per cross sectional area [START_REF] Farber | Use of x-ray tomography to study the porosity and morphology of granules[END_REF][START_REF] Vladisavljevi¢ | Permeability of hydrophilic and hydrophobic Shirasu-porous-glass (SPG) membranes to pure liquids and its microstructure[END_REF]. XRCT results are also compared with the ones obtained from mercury intrusion porosimetry. One of the most interesting parameter is the pore size distribution analyzed directly from the tomographic images, as illustrated in Fig. 1.24 [START_REF] Farber | Use of x-ray tomography to study the porosity and morphology of granules[END_REF].

By analyzing dierent samples with the two techniques, the authors concluded that XRCT and mercury intrusion porosimetry yield very similar overall porosities; noteworthy is the fact that only the largest pores (so-called macropores) were analyzed [START_REF] Farber | Use of x-ray tomography to study the porosity and morphology of granules[END_REF]. Finally, the porosity of granular materials can be characterized in many possible ways and using various metrics, such as bulk and local void ratio, spatial distribution of density, number and size of pores, etc. These quantities do not only complement each other but also can be favourably injected into computer simulations in general, and the DEM in particular in order to provide a detailed comparison with XRCT.

Microstructure characterization

The characterization of the contact network of granular materials is very important to understand how the material behaves under an external applied force. One of the most important characteristic of the contact network is the coordination number of the system.

The coordination number is dened as the number of contacts per grain. A contact network can be characterized not only by the total coordination number of the whole structure but also by the local values of the number of contacts. In fact, the coordination number has been characterized in many studies using XRCT [START_REF] Seidler | Granule-by-granule reconstruction of a sandpile from X-ray microtomography data[END_REF][START_REF] Al-Raoush | Extraction of physically realistic pore network properties from three-dimensional synchrotron X-ray microtomography images of unconsolidated porous media systems[END_REF][START_REF] Neethirajan | Investigation of 3D geometry of bulk wheat and pea pores using X-ray computed tomography images[END_REF][START_REF] Bhattad | Eect of Network Structure on Characterization and Flow Modeling Using X-ray Micro-Tomography Images of Granular and Fibrous Porous Media[END_REF][START_REF] Andò | Experimental micromechanics: grainscale observation of sand deformation[END_REF], in studies using XRCT and optimization methods [START_REF] Al-Raoush | Simulation of random packing of polydisperse particles[END_REF][START_REF] Al-Raoush | Microstructure characterization of granular materials[END_REF], and many studies using the DEM. However, one of the most interesting combination of XRCT and DEM for analyzing the coordination number can be found in Marmottant et al. ). This discrepancy is due to the digitization procedure, which articially increases the measured area (because of the nite voxel volume), causing for example a nonzero area for a point-contact between two particles [START_REF] Marmottant | Coordination measurements in compacted NaCl irregular powders using X-ray microtomography[END_REF].

The determination of the coordination number in granular media is also a typical subject in stereology and algorithms for extracting this parameter from digital images [START_REF] Jouannot-Chesney | Practical Determination of the Coordination Number in Granular Media[END_REF].

Nevertheless, the estimation of the coordination number and the characterization of the contact network have not been extensively investigated in engineering applications in spite of the key role of these parameters to understand comprehensively the behavior of granular materials under dierent loading conditions.

Radial distribution function

Radial distribution function (RDF) or pair correlation function is also a common way of comparing the structures obtained from simulations with those from experiments. The RDF is the probability of nding a particle center in a spherical shell. Many examples of the use of RDF can be found in the literature [START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF]. The detailed information of the RDF will be also described in Chapters 2 and 3. In this section, we only present a typical study using both of XRCT and DEM for the RDF.

The inuence of particle shape on the local packing of particles is studied by Fu et al. [START_REF] Fu | Investigation of particle packing in model pharmaceutical powders using x-ray microtomography and discrete element method[END_REF][START_REF] Fu | Application of x-ray microtomography and image processing to the investigation of a compacted granular system[END_REF]. These authors observed that the rst peak of the RDF, obtained from XRCT, was lower than the one obtained from DEM simulations, as illustrated in Fig. 1.26.

The oscillations in the RDF extended through a longer distance in densely packed states and were quickly damped in the loosely packed states. This fact was attributed to the presence of large disordered structures when the system was loosely packed. Other studies also observed the RDF, such as the study of bead compression [START_REF] Richard | Analysis by x-ray microtomography of a granular packing undergoing compaction[END_REF], characterization of packing structure [START_REF] Seidler | Granule-by-granule reconstruction of a sandpile from X-ray microtomography data[END_REF], and the analysis of the inuence of particle Chapter 1. Literature review shape [START_REF] Al-Raoush | Microstructure characterization of granular materials[END_REF]. In summary, the RDF obtained from XRCT can provide extensive information for the corroboration of the DEM.

Particle displacements

The determination of particle displacements is usually applied in XRCT and DEM to visualize the changes in the structure of a granular material. Although in most cases, In many studies using XRCT, the 2D/3D Digital Image Correlation (DIC) method has been extensively applied to analyze the micro-macroscopic displacements of particles under dierent loading conditions [START_REF] Hall | Discrete and continuum analysis of localised deformation in sand using X-ray µCT and volumetric digital image correlation[END_REF][START_REF] Watanabe | Displacement in sand under triaxial compression by tracking soil particles on x-ray ct data[END_REF][START_REF] Andò | Experimental micromechanics: grainscale observation of sand deformation[END_REF][START_REF] Andò | Grain-scale experimental investigation of localised deformation in sand: A discrete particle tracking approach[END_REF][START_REF] Viggiani | Identifying and following particle-to-particle contacts in real granular media: An experimental challenge[END_REF]. This approach provides an ideal framework to compare experiments on granular materials and DEM simulations at both macroscopic and microscopic scales.

In our study, the DEM and the XRCT are chosen to observe the mechanical behavior of wet granular materials from macroscopic level to microscopic level. We performed oedometric compression tests in order to investigate the dierences between the irreversible behavior observed experimentally and the one predicted by numerical simulations.

Conclusion

In this chapter, we have reviewed the mechanical behavior of bonded granular materials which were observed by using both experimental and numerical investigations, with a particular focus put on wet granular materials. A brief description of the behavior of granular materials is also given. Three common non-destructive experimental methods and several computer simulation methods which were usually applied to approach the behavior and the structure of granular materials, from macroscopic to microscopic levels, are succinctly presented. Among various numerical methods, the DEM is one of the most popular applied to investigate the behavior of wet granular materials. We have introduced the following remarks: (i) adhesive granular materials have been observed less frequently than cohesionless ones; (ii) dense states of granular materials were studied signicantly more than loose or very loose states; (iii) the association of the DEM and the XRCT is still limited in the study of granular materials in general and wet granular materials in particular. Therefore, the work presented here will combine DEM and XRCT to investigate wet granular materials, in very loose states stabilized by the presence of a very small amount of interstitial liquid. In the following chapters, we will present experimental and numerical evidences, which will be analyzed at both macroscopic and microscopic scales.

CHAPTER 2 Experimental investigations

In this chapter, we present the observations of grain-scale behavior of a wet granular material in very loose state under one-dimensional compression tests, combined with X-ray computed tomography. We rst describe the global compaction behavior of this material. A procedure to process images is then presented. It involves an algorithm to detect spherical structures from 3D tomography images, and the analysis of microstructure properties. A link between micro and macroscopic behavior is nally discussed to get a better insight into the macroscale behavior based on the microscale one.

Introduction

So far, many authors investigated the macro-and micro-mechanical behavior of wet granular materials [START_REF] Pierrat | Tensile strength of wet granula materials[END_REF][START_REF] Mason | Critical angle of wet sandpiles[END_REF][START_REF] Fournier | Mechanical properties of wet granular materials[END_REF][START_REF] Mitarai | Wet granular materials[END_REF][START_REF] Jiang | Summary of collapsible behaviour of articially structured loess in oedometer and triaxial wetting tests[END_REF][START_REF] Bruchon | Full 3d investigation and characterisation of capillary collapse of a loose unsaturated sand using x-ray ct[END_REF] (among others). These works evidenced the link between the macroscopic mechanical behavior and the changes of microstructure.

However, few works have focused on the grain-scale behavior.

Moreover, among various non-destructive experimental techniques, the XRCT is a high-resolution non-destructive 3D observation method, that does not require any sample pre-treatment. This technique can advantageously be combined with mechanical loading tests. Such a combination constitutes a clear strength of this method. Therefore, several recent studies used this technique to investigate the microstructure of granular materials with / without combining mechanical loading or any other test [START_REF] Farber | Use of x-ray tomography to study the porosity and morphology of granules[END_REF][START_REF] Aste | Investigating the geometrical structure of disordered sphere packings[END_REF][START_REF] Aste | Variations around disordered close packing[END_REF][START_REF] Aste | Geometrical structure of disordered sphere packings[END_REF][START_REF] Al-Raoush | Simulation of random packing of polydisperse particles[END_REF][START_REF] Al-Raoush | Microstructure characterization of granular materials[END_REF][START_REF] Andò | Experimental micromechanics: grainscale observation of sand deformation[END_REF][START_REF] Andò | Grain-scale experimental investigation of localised deformation in sand: A discrete particle tracking approach[END_REF][START_REF] Bruchon | Full 3d investigation and characterisation of capillary collapse of a loose unsaturated sand using x-ray ct[END_REF]. However, most of these works focus on the analysis of dense packing of cohesionless particles of large size.

In this chapter, we present an experimental investigation of a wet granular material assembled in a very loose state by means of XRCT. A cylindrical sample constituted of thousands of glass beads (average diameter: 100 micron) with liquid bonds in the pendular state was submitted to an oedometric compression test, which consisted in various compression steps. At every step, a XRCT scan was performed.

Material and methods

Material

The studied material is composed by spherical glass beads, an industrial product. Beads with diameters ranging from 80 µm to 104 µm were chosen. With these diameters, particle weight looses relevance with respect to capillary forces [START_REF] Santamarina | Soil Behavior at the Microscale: Particle Forces[END_REF]. Furthermore, the • The glass beads were mixed with water at an initial water content of w 0 = 3.5 %,

to ensure that the specimens remained in the pendular regime at any time.

• The oedometric cell was placed on the support system (see Fig. 2.4). It was then lled with the wetted glass beads by using a controlled pluviation method, described hereafter. The lling was performed by placing the wetted glass beads on a xed sieve with a 200 µm opening. The sieve was vibrated vertically with a vibration machine. Thanks to the vibration, clusters of glass beads passed sequentially through the sieve to nally fall into the cell. The height of free fall (h f ) from the sieve bottom to the cell bottom was varied between 80 mm and 260 mm.

• All specimens were prepared in an air-conditioned room. The preparation of each specimen took approximately from 30 to 45 minutes.

• The top surface of the specimen was then perfected by using a small iron wire to cut the excess soil (see Fig. • After the compaction test, the water content of the specimens was measured by oven-drying (w ap ).

The initial parameters of specimens after preparation are listed in after preparation is strongly aected by the height of free fall (h f ). For instance, the water content after testing is w ap = 1.60 % for the highest height (h f = 260 mm), and w ap = 2.71 % for the smallest height (h f = 80 mm). Indeed, with same initial conditions, the evaporation of samples is proportional to the height of free fall. Clear evidences are given in Table 2.1: the maximal evaporation amounts to 1.90 % of water content (relative reduction of ≈ 54.29 %) for the highest height and the minimal evaporation corresponds to 0.79 % of water content (relative reduction of ≈ 22.57 %) for to smallest height.

Although the water contents after preparation were quite small, the specimens were still stable in very loose states. The calculated degree of saturation (S < 5 %) conrms that the specimens are always maintained in the pendular state. This conrms that the capillary forces between grains play a key role to maintain the stability of specimen's structure.

Moreover, the relationship between initial void ratio e 0 , initial porosity n and the height of free fall h f is illustrated in Fig. 2.5. The circle-red points represent the specimens with cell φ50, the square-blue points represent the specimens with cell φ20, and the triangle-green points represent other specimens which were not used for further analysis. Fig. 2.5 conrms that the higher is the height of free fall, the smaller the initial void ratios are.

2.2. Material and methods 57

Oedometric compression tests

After its preparation, the specimen was carefully placed on the loading platform of the compression system, as shown in Fig. 2.6. In this system, the piston was xed to the load cell which has a maximum capacity of 1,000 N. A linear variable dierential transformer (LVDT) was installed to measure the displacement of the loading platform. Before Chapter 2. Experimental investigations recording the signals, the cell was moved upward in order to approach the piston. A rate of 0.5 mm/min was set up for the loading cycle. The maximal applied force was limited approximately to 700 N. All signals of load cell and LVDT were received by the data logging system using the software LabVIEW. After compression, all specimens were taken out from the cell (see Figs. The XRCT images are reconstructed by using the GPU-based reconstruction software Xact developed by RX-Solutions. Some lters may be used to correct some artifacts as ring artifacts.

Cell

Apparatus

Apart from the cell φ20, several gaskets and pistons were also designed for performing XRCT scans, as shown in Fig. 2.9. The height of pistons were manufactured with dierent values in millimeter, as written in Fig. 2.9(a). Several wedges in paper were

Scanning procedure

The scanning procedure combined with an in-situ compression test for the sample CT03 (mentioned in Table 2.1) is now described. The specimen was scanned at various stages of compression; the order of scanned compressed heights is depicted in Fig. 2.10. The detailed scanning program is given in Table 2.2, in which the void ratio was calculated from the compressed height and initial parameters of sample CT03. The height of pistons δ is given in Fig. 2.9.

At each level of the compressed height, a 3D image was scannned at the scale of the whole specimen (scans hereafter referred to as global scan or GS). Besides, at the levels corresponding to heights h 1 , h 3 , h 5 , and h 7 , 3D images of the central part of the specimen (see the 10-mm diameter dashed-region in Fig. 2.10) were also scanned at a ner scale for a more detailed view of the grains and the pores (scans hereafter referred to as local scan or LS). For the global scans, X-ray source voltage and current were respectively 70 kV and 130 µA; the voxel size was 13 µm and a single scan lasted 80 min. For the local scans, X-ray source voltage and current were respectively 70 kV and 95 µA; the voxel size was 6.5 µm and a single scan lasted 240 min. 3D images were saved in 16 bits and subsequently converted into 8 bits in order to reduce memory usage for subsequent image analyses. Void ratio, e (-)

h 0 = 10 mm h 0 h 1 h 2 h 3 h 4 h 5 h 6 h 7 3 mm d
Vertical stress, σ v (kPa) Based on the shape of these curves and the range of vertical stress, three stages can be identied. During the rst stage, under low stresses (< 0.1 kPa), void ratio e slightly decreased. Void ratio then sharply decreased in the intermediate stage. Finally, in the third stage, void ratio gradually decreased for vertical stresses larger than 10 kPa. In particular, all curves (except those for cell φ20) are very close in this stage, irrespectively of the initial void ratio. In other words, when almost no further rearrangement is possible, the decrease of void ratio was very small and almost limited to the elastic deformation of grains at the contact points.

OD01 OD02 OD03 OD04 OD05 OD06 OD07 OD08 CT01 CT02
For the cell φ20, upon growing of the vertical stress, the decrease of void ratio is slower than for the cell φ50. This can be explained by the eect of the diameter (20 mm against 50 mm) and the height (10 mm against 20 mm) of the cell, which tend to increase lateral friction.

It should recognized that the macroscopic behavior of the studied material under 1D compression tests is very special. Indeed, in soil mechanics [START_REF] Wood | Soil Behaviour and Critical State Soil Mechanics[END_REF][START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF][START_REF] Tang | Analysing the form of the conned uniaxial compression curve of various soils[END_REF][START_REF] Muñoz Castelblanco | Some aspects of the compression and collapse behaviour of an unsaturated natural loess[END_REF], only the rst two stages are observed in oedometric compression curves. Interestingly, the last stage 2.4. Global scan observations 63 observed in the compaction curve in the present study was also described in the Sshape curves obtained by Tang et al. [START_REF] Tang | Analysing the form of the conned uniaxial compression curve of various soils[END_REF]. Similarly, three stages are also observed in the numerical results of isotropic compression with the same material [START_REF] Than | Basic mechanical properties of wet granular materials: A DEM study[END_REF].

Global scan observations

Eight global scans (Table 2.2) were successively captured at various the compaction steps, and were named as follows: S01G, S02G, S04G, S05G, S07G, S078G, S10G and S11G. 3D image analyses were carried out by using the open source software Fiji.

3D reconstruction from the image stack

After converting to the 8-bit images, the 3D sample volume at dierent states of compaction was reconstructed, as shown in Fig. 2.12. These 3D images, with diameter 20 mm and height h i , correspond to the global scans given in Table 2.2. A quarter of these 3D images was removed to observe the inside of the specimen. Pores and grains can be distinguished thanks to grey levels. Lighter grey levels typically correspond to more absorbing materials, while darker grey levels correspond to less absorbing materials. There exist a lot of segmentation methods in the Fiji software. In the present study, the IsoData algorithm, proposed by Ridler and Calvard [START_REF] Ridler | Picture Thresholding Using an Iterative Slection Method[END_REF], was chosen for all analyses. The procedure divides the image into object and background by taking an initial threshold, then the averages of the pixels at or below the threshold and pixels above are computed. The averages of those two values are computed, the threshold is incremented and the process is repeated until the threshold is larger than the composite average. That is, threshold = (average background + average objects)/2. As shown in Fig. 2.12, it appears clearly that the density is not homogeneous within the cylindrical sample, especially at the periphery where more voids are observed. A larger density with less voids at the top and bottom of the specimen is also observed.

Global scan observations 67

There exist also several large holes in the center of the specimen in the rst compression steps (S01G, S02G, S04G, and S05G). This is due to the specimen preparation, and also to the boundary conditions imposed by the ring. In other words, during the specimen preparation, the friction between the external ring and the grains opposes to the displacement of other particles. The cohesion of grains also plays a remarkable role to limit the movement of particles.

At higher levels of compaction (Figs. 2.12(e, f, g, and h)), the structure is strongly rearranged in the central layers because the densest state is obtained in this region (horizontal band). Reversely, the structure is slightly aected at the bottom and top layers.

Moreover, the average grey level distributions of all global scans were also extracted to observe the change of density during compaction. Fig. 2.15 illustrates the average grey level according to the elevation (h) of all scans during compaction. In general, the specimen is not homogeneous. Initially, the specimen is fairly homogeneous from h = 1 mm to h = 7 mm (from bottom to top of specimen). After a slight compaction (S02G), the density of the upper part of specimen slightly increases while the one at the lower part (from 0 to 2 mm) is not changed. The upper part of specimen is continuously compressed at higher compaction (S04G and S05G) while the lower part is negligibly disturbed. For the highest compaction (S07G, S08G, S10G, and S11G), the central layers of the specimen are strongly settled, while the bottom layers are insignicantly rearranged.

Spatial distribution of void ratio

In For the diagram 1 (Fig. 2.16(a)), the void ratio distribution according to the elevation for dierent compaction steps is depicted in Fig. 2.17. At the initial state (S01G and S02G), from bottom to top, the void ratio is largest at the sixth layer (e max = 2.65)

because of the inuence of the wall in the specimen preparation process. Meanwhile the sample is denser at the bottom layers and densest at the top layers. This is caused by the end of preparation process when the grains were widely spread on the surface of the cell. The top layers are therefore more homogeneous than the other layers. These observations conrm that the void ratio in this region is the smallest initially. For the higher compaction level (S04G and S05G), the structure is notably compacted and rearranged. Void ratio in the central and top layers is reduced, especially the central layers of scan S05G. Up to scans S07G and S08G, the structure is signicantly perturbed.

h i R i 1 2 5 7 8 V1 V2 V3 V4 V5 X Y (a) (b) (c) (d) (e)
The void ratio in the central layers strongly decreases (from e = 2.25 to e = 1.15, ≈ 49.80 %). This reduction is larger than the one observed in the bottom layers and the top layers (from e = 1.51 to e = 1.15, ≈ 23.25 %). At the densest states (S10G and S11G), a strong decrease of void ratio is obtained in the central layers (e min = 0.7) while the structure in the bottom layers is slightly rearranged. Actually, when the specimen becomes denser in the center, the corresponding layers seem to prevent the displacement of the top layers. This explains why it is really dicult to compact this specimen at S11G although the total void ratio is decreased by approximately 50 % (from e = 2.23 to e = 1.11). This is already conrmed via the decrease in the total void ratio in Fig. 2.11 from the loosest state (S01G) to the densest state (S11G).

For the diagram 2 (Fig. 2.16(b)), Fig. 2.18 depicts the spatial distribution of void ratio which is investigated as a radial function, considering dierent concentric hollow cylinders having the same height for each scan. In general, the mean distribution of void is rather regular, except the center and periphery of the cylinders in the rst four scans.

Once again, this demonstrates the remarkable eect of the wall friction which is already mentioned above. In the last four scans, void ratio e variation with the distance to the center is negligible. This reects the homogeneity of the specimen under the higher For the diagram 3 (Fig. 2.16(c)), the void ratio distribution according to eight equal sectors for dierent compaction steps is illustrated in Fig. 2.19. In general, the distribution of void ratio is fairly regular at the initial state, except fro sectors S1 and S8 in rst four steps (S01G, S02G, S04G and S05G). In these sectors, there exist several large pores between grains and the wall (see Fig. 2.12). For the denser states (last four steps), the structure is really homogeneous. Besides, the distribution of pore volume is also analyzed. The algorithm thickness [START_REF] Dougherty | Computing local thickness of 3D structures with ImageJ[END_REF][START_REF] Hildebrand | A new method for the model-independent assessment of thickness in three-dimensional images[END_REF] in the plugin BoneJ of Fiji is applied to compute the local thickness of 3D structures in 3D images. It denes the thickness at a point as the diameter of the greatest sphere that ts within the structure and which contains the point. In other words, the grains and the pores are precisely separated by this algorithm. The distribution of pore volume in volume V2 for each of two adjacent scans is illustrated in Fig. 2.24. Numbers of pores mean the number of all the pores which have the same size (same colored voxels). In general, the shape of the distribution of pore volume in all cases is almost similar. More small pore volumes are created under growing of the compaction step. In each of two adjacent scans, the larger pore volumes are negligibly decreased while the smaller ones are signicantly reduced (see the continuous line and the dashed line in Fig. 2.24). For instance, in the initial stages (S01G and S02G red curves), the smallest pore volume varies from 3×10 -4 mm 3 to 2×10 -3 mm 3 , and in the nal stages (S10G and S11G green curves), the smallest one varies from 10 -6 mm 3 to 2 × 10 -4 mm 3 . Interestingly, in the larger compaction levels (S07G and S08G violet curves), several pore volumes from the scan S08G (dashed curve) are larger than the ones from the scan S07G. This reects a strong perturbation of the structure in these stages. The evidence in Fig. 2.17 also conrms this phenomenon. Image processing was also performed by using the open-source software Fiji. Besides, the author's codes (using Matlab software, the programming language Python and the open-source software ParaView) were mainly used for the image analysis. The algorithm used to detect spherical structures of 3D images, developed by Xie et al. [START_REF] Xie | Magnetic resonance histology of age-related nephropathy in the Sprague Dawley rat[END_REF] in 3D is applied with a suitable development. This algorithm was modied from the one of Peng et al. [START_REF] Peng | Algorithms for On-Line Monitoring of Micro Spheres in an Optical Tweezers-Based Assembly Cell[END_REF] in 2D and the Adaptive Hough Transform algorithm (original algorithm) of Illingworth and Kittler [START_REF] Illingworth | The Adaptive Hough Transform[END_REF]. Furthermore, a 3D skeletonization algorithm, proposed by Lee et al. [START_REF] Lee | Building skeleton models via 3D medial surface/axis thinning algorithms[END_REF] and then developed into a plugin of BoneJ in Fiji by Doube et al. [START_REF] Doube | BoneJ: Free and extensible bone image analysis in ImageJ[END_REF],

is also applied to plot the contact network of particles.

Algorithm to detect spherical structures

The algorithm to detect spherical structures from tomography images is illustrated in Fig. 2.27. A 3D image with a cubic shape, and formatted as a 3D array, was rst extracted from the local scans. The center's positions of spheres were then directly located within the 3D array. Afterward, the radii of the spheres (with identied central position) were identied. Finally, the 3D reconstruction process and analysis were carried out from the detected spheres.

The detailed description of the main algorithm (the dashed-red square in Fig. 2.27) will be given in the following subsections.

Create a 3D image in

Locating the center's position of spheres

The method to compute the center's position of a sphere in the 3D image is rst to identify the voxels that belong to that sphere, and then to compute the centroid of those voxels. The gradient eld of image intensity is rstly computed using the following equation:

∇I(i, j, k) = (V x , V y , V z )| (x,y,z) = I(i, j, k) -I(i, j, k -1), I(i, j, k) -I(i, j -1, k), I(i, j, k) -I(i -1, j, k) (2.1)
where (i, j, k) are the voxel indices, ∇I(i, j, k) is the gradient vector at voxel (i, j, k), which consists of x, y and z components, and I(i, j, k) is the image intensity at voxel (i, j, k).

In a gradient eld, the nonzero gradient vectors are either pointing toward the center of a sphere or away from it. So, a transform is dened to convert the gradient eld to an accumulation array, in which the voxel intensity corresponds to the probability of that voxel being the center of a sphere. The accumulation array has the same dimension as the gradient eld. In the accumulation image, the maximum intensity represents the center position of a sphere. The accumulation array is constructed by an adding process, as shown in Fig. 2.28. In the gradient eld, for each nonzero gradient vector ∇I(i, j, k), In order to determine the radii of spheres, a signature curve f (r) is dened over the gradient eld of the sphere's image, where r is the distance to the center position of the sphere and f (r) is the averaged dot product of the gradient vector and the radial vector. For a voxel (i, j, k) in the neighborhood of (i c , j c , k c ), let V (i, j, k) be the gradient vector of the voxel, (i, j, k) and q(i, j, k) be the vector joining (i c , j c , k c ) to the center of the voxel (i, j, k). A sequence of discrete r values with xed interval is selected. It represents distances to the sphere center (i c , j c , k c ). For each r value, f (r) is dened by the following equation:

a
f (r) = Σ || q(i,j,k)|-r|<∆r/2 V (i, j, k). q(i, j, k) | q(i, j, k)| Σ || q(i,j,k)|-r|<∆r/2 1 , (2.2)
where the sum is over all the voxels whose distance from (i c , j c , k c ) is within the range r -∆r/2, r + ∆r/2 .

Estimating the lost contacts after the detection process

After the detection process, the spheres which have a part located outside of the border of the cube were not identied (see Figs. 2.25(c) and 2.31). Hence, the contacts between these spheres and their neighbors inside the standard volume were lost. In order to estimate the lost contacts, a method is proposed as follows. A new cube which is larger than the SV (see denition in Sec. EVs having the same coordinates were rst identied. The grains, around the border of the EVs, which have contacts with the ones inside the SVs, were then distinguished.

The lost contacts N l C of the SVs were nally found. They will be used in the analyses of geometrical changes in the following sections. Hence, they cannot be seen in the cubes on the right-hand side in Fig. 2. 33 (b, d, f, andh). During the detection process, all the SVs were set up with the same input parameters.

3D reconstruction after detection process

The results show that the average number of particles at a denser state is higher: at the loosest state, S03L (e = 2.23), there are approximately 4000 particles in one cube, while at the densest state, S12L (e = 1.04), there are approximately 7500 particles. Besides, the macroscopic solid fraction and the average solid fraction of ten SVs are illustrated in Fig. 2.37. The triangle points are computed from the void ratio of the whole specimen (see Sec. 2.4.3). The circle points show the average values, which were directly calculated from the tomography images. In addition, after the detection process, the solid fraction is easily computed by the volume of all particles and the volume of the cube. The results show that the macroscopic solid fraction is dierent from that estimated at the center of the specimen. That can be partly explained by the heterogeneity of the microstructure of the sample. Besides, the solid fraction obtained after the detection process is smaller than that before the detection; the decrease equals to 0.05, 0.06, 0.07, and 0.1 for the scans S03L, S06L, S09L, and S12L, respectively. It seems that the detection process underestimate the radii of the spheres and this problem is more serious at a denser state.

Furthermore, the solid fraction after the detection process is also recalculated with a tolerance of 1 voxel (6.5 micron) added in the diameter of spheres (the diamond-shaped- scans before and after the detection process with a tolerance of 0 and 1 voxel. BD denote before detection, AD denotes after detection. magenta points). These values are very close to the ones before detection (the circle-red points). This conrms that the use of a tolerance of 1 voxel for all analyses of the 3D images after the detection process is quite suitable.

Geometrical structure analysis 2.5.7.1 Total coordination number

In granular materials, the coordination number z is dened as the average number of force-carrying contacts per grains. If N C is the number of force-carrying contacts between N grains, then the coordination number is z = 2N C /N as product zN counts each contact twice. In the present analysis, N C is redened as the number of pairs of contacts between N grains. Pairs of contacts mean two particles which have a distance between two centers smaller than or equal to the sum of its radii. In order to reduce the error of the rounded value of diameters in voxel during the detection process, the eect of the resolution of tomography images, a tolerance in voxel value is used when estimating the coordination number. This tolerance is dened as the distance between two surfaces of particles. Therefore, in the present part, z is called the total coordination numbers. As mentioned in Sec. 2.5.2.3, the diameters after the detection process are smaller than the real diameters of material. Hence, a tolerance of 1 voxel is chosen for the reference error of all calculations. stress. This is not reasonable but this conrms once again the error mentioned above. z then increases twofold and nearly threefold when the tolerance changes from 0 to 1 voxel. However, from 1 to 2 voxels, z is gradually raised from a half to nearly one and a half.

The error bars are wider at a denser state or with a higher tolerance.

The new z (with lost contacts) of four steps is also plotted in Fig. 2.38 for the tolerance of 1 voxel (the triangle-magenta points). The results show that this correction (i.e. using lost contacts) increases the total coordination number by approximately 10 % for the four cases. For the following analysis, the lost contacts were not considered for the sake of simplication.

Coordination number of close neighbors

The coordination number of close neighbors, z(h), is the average number of neighbors of one particle separated by an interstice narrower than h. z(0) is the usual total coordination number z with tolerance 0. (

2.3)

This distance is equivalently from 0 to 2 voxels.

For each scan, there is a threshold of h below which z(h) increases quickly when h increases; above this threshold, the z(h) increases with a lower rate. Based on the shape of these curves, we can suggest that this threshold is around 1 voxel (equivalently h/ d ≈ 0.07). At this point, four curves start to diverge and the z(h) of the lower compaction level is smaller than the one of the higher compaction level. This threshold, applied for all the scans, corresponds to the error induced by the detection process. Below the threshold, the rate at which z(h) increases with h is independent of the compaction state because it is related to the error of the detection process. Above the threshold, the rate at z(h) increases with h is higher for a denser state. Actually, this rate depends on the structure of the samples. The denser the states, the higher the numbers of close neighbors.

Number of contacts per grain

Fig. 2.40 shows the number of contacts per grain for the three tolerances 0, 1, and 2 voxels. In general, the grains have from one to eight contacts. Firstly, with a tolerance of 0 voxel, the grains having zero contact occupy the highest quantities, from 30 % to nearly 40 %. The denser the stages, the higher the percentage of grains without contact. This is due to the error of the detection algorithm as already mentioned above.

For the tolerance of 1 voxel, the grains having two and three contacts occupy the most quantities, about 27.5 % for each series. Meanwhile, the quantities of the grains without contact decrease from threefold to fourfold corresponding to growing of compaction. For the denser state, the probability to have a high value of N C (3, 4, 5, and 6) is higher than for the other states. Reversely, the probability to have low values of N C (0, 1, and 2) is lower for denser states than for looser states. The probability to have N C = 0 is approximately 6 % to 8 %. That can be attributed to the lost contacts mentioned above, which are not considered in this analysis, where the red spheres have N C = 0 (see Figs. 2.35(e) and (g)). They are mainly located close to the border of the cube. N C = 7 is the maximal value and it occupies a minimal probability, nearly 0.2 %. For the tolerance of 2 voxels, the probabilities P (N C ) of four states are signicantly changed for the ones of tolerance of 1 voxel, in which the P (N C ) at N C = 3 is negligibly changed. For the denser state, low value of N C (0, 1, and 2), the P (N C ) is notably decreased while it is considerably increased at high value of N C (4, 5, and 6). The grains having nine contact exist but the corresponding quantity is very small, approximately 0.004 %.

Several authors reported the distribution of number of contacts (i.e. the coordination number) [START_REF] Al-Raoush | Microstructure characterization of granular materials[END_REF][START_REF] Al-Raoush | Extraction of physically realistic pore network properties from three-dimensional synchrotron X-ray microtomography images of unconsolidated porous media systems[END_REF][START_REF] Aste | Local and global relations between the number of contacts and density in monodisperse sphere packs[END_REF][START_REF] Marmottant | Coordination measurements in compacted NaCl irregular powders using X-ray microtomography[END_REF][START_REF] Moreno-Atanasio | Combining X-ray microtomography with computer simulation for analysis of granular and porous materials[END_REF] by using dierent methods to detect the contacts between the spheres. However, calculating the number of contacts is dicult due to the uncertainty (due to the diraction of X-rays at the boundary) when calculating the accurate position of the center of the particle, especially with the small size of grains and low resolution of the tomography images. That is why the tolerance of 1 voxel is chosen for all calculations in the present study. This tolerance is close to the interval tolerance by 1-2 voxels 1/3 of Aste et al. [START_REF] Aste | Geometrical structure of disordered sphere packings[END_REF][START_REF] Aste | Local and global relations between the number of contacts and density in monodisperse sphere packs[END_REF] while their size of particles is bigger than our grains, and the resolution is not too much dierent. In several packing studies [START_REF] Al-Raoush | Simulation of random packing of polydisperse particles[END_REF][START_REF] Aste | Variations around disordered close packing[END_REF][START_REF] Aste | Local and global relations between the number of contacts and density in monodisperse sphere packs[END_REF],

the spheres are at maximum in contact with 11 particles and at minimum with 2. Meanwhile, in our study, the maximum number of contacts is 8 with an average percentage of 0.01 %, and there still exists a signicant percentage (around 7 %) of spheres that have no contact. We think that this can be explained by two reasons. Firstly, during the specimen preparation, the spheres are automatically stuck in groups of two or more grains after passing through the sieve. When falling into the cell, the cohesion of spheres tends to resist the interpenetration of the other groups, especially in some closed pores.

Hence the structure of the specimen is not homogeneous. Secondly, as mentioned in Sec. 2.5.2, the spheres located close to the borders of the analyzed volumes are partly outside the analyzed cubes and thus not detected. Hence some spheres which are in contact with incomplete spheres will be isolated after the detection process.

Interestingly, the sphere having one contact can exist in this structure, as illustrated in Fig. 2.41. This is a small cube extracted from the Fig. 2.35(a). Furthermore, from the loosest state (S03L) to the densest state (S12L), N C = 1 occupies a signicant amount, approximately 30 % (tolerance of 0 voxel), 17 % (tolerance of 1 voxel), and 11 % (tolerance of 2 voxels). These spheres (in green) can be clearly seen in Figs. 2.35(e) and (g). This is explained by the fact that the capillary force between spheres is higher than the weight of the spheres. Indeed, with the size of diameters is 100 µm, the particle weight looses relevance with respect to capillary forces [START_REF] Santamarina | Soil Behavior at the Microscale: Particle Forces[END_REF].

Radial distribution function

The Radial Distribution Function (hereafter referred to as RDF), g(r), is the probability of nding a particle center in a spherical shell of a radius r, given that there is another particle at the origin of the spherical shell [START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF]. The RDF is related to the average number of particles by p(r 2 ) -p(r 1 ) = r 2 r 1 g(r)4πr 2 dr, (2.4) where p(r) is the average number density of particles in a spherical shell of radius r. In this formula, the radial distances will be normalized by the average diameter d (given in Eq. 2.3). The RDF is a useful tool to describe the structure of a system, and it has been used widely for granular materials, especially in the numerical simulations. the RDF must be zero, because two particles cannot occupy the same space.

The RDF for S03L is rougher than for the other ones. The functions start from maximum probability of the contact particles at rst peak r = d with dierent amplitudes. The functions then decrease to minimum probability at r = 1.2 d , r = 1.3 d , r = 1.35 d , and r = 1.4 d upon growing of compaction S03L, S06L, S09L, and S12L, respectively. These functions increase subsequently to another local peak at r ≈ 1.9 d . From this peak, these functions continue to vibrate slightly until they approach a constant value at r = 3.0 d . It is clear that g(r) should go to 1 for large r. However, these RDFs in Fig. 2.42 go asymptotically to 1 for large r. This can explain why spherical shell cannot include the spheres which are located in the corner of the cube. Besides, due to the heterogeneity of the SVs, these RDFs slightly oscillate until their value is close to unity.

Chapter 2. Experimental investigations 2.5.8 Contact networks Fig. 2.43 illustrates the contact networks of a typical SV, represented in a narrow slice of nearly four average particle diameters (4 d ≈ 50 voxels) thick before detection (on the left-hand side) and after detection (on the right-hand side). Each pair (before and after detection) from (a and b) to (g and h) corresponds to growing of compaction, from S03L to S12L. The contact between particles are represented by segments of dierent colors (only for after detection) joining particle centers.

Before detection, the contact network is directly analyzed from the 3D binarized images by joining the particle centroids. On the left-hand side of Figs. 2.43(a, c, e, and g), there are still a lot of small blocks. This is due to some small pieces of incomplete spheres at the borders that are not identied. After detection, the contacts are plotted in two types (on the right-hand side of Fig. 2.43). The rst one (red lines) is the contacted particles including one voxel of diameter (tolerance of 1 voxel). The second one (blue lines) is the neighbors with their interstitial distance of 1 voxel.

In general, the contacts obtained before the detection process are more than the ones after detection. This is due to the grains loss during the detection process, already mentioned in the above sections. Moreover, the lines, joining two particle centroids on the left-hand side, are spatially connected from many voxels so the contact networks are rather complex for observing. This conrms that these lines are not straight. Furthermore, on the right-hand side of Fig. 2.43, the contact between neighbors with their interstice of 1 voxel are also drawn in order to reduce the eect of detection process. Indeed, the general shape of the contact networks before and after detection are fairly similar.

A supplemental illustration of contact networks between grains is shown in Fig. 2.44.

A typical slice (4 d ) of the typical SV is extracted. It has the same position of the thickness in Fig. 2.43. For depicting the number of contacts per grains, a tolerance of 1 voxel is computed to the diameters of particles. The number of contacts per grains are illustrated by dierent colors. The red, green, blue, dark-yellow, magenta, dark-cyan, dark-green, dark-purple, black and dark-red represent the spheres having zero, one, two, three, four, ve, six, seven, eight, and nine contacts, respectively. In these slices, the dark-purple, black, and dark-red spheres do not exist but the dark-green spheres having 6 contacts can be easily found. Nevertheless, the red spheres still exist signicantly. Therefore, we can conclude that the number of contact is underestimated because of the signicant inuence of the lost contacts.

Links between macro and microscopic behaviors

The links between macro and microscopic behavior of specimen under growing of compaction can be observed via the compression curves, the global and local observations.

In other words, these links are successively depicted by the Figs. In the stage 1 of the compaction curve (Fig. 2.11), the void ratio remains nearly constant and slightly decreases at the end of this regime (equivalently the position of the scan S02G). This is depicted well by the distribution of void ratio of scans S01G and S02G in Fig. 2.17.

Under growing of stress in stage 2, the structure of specimen quickly collapses. This is evidently reected by a twofold decrease of void ratio in the central layers of specimen (see Fig. 2.17). At these regions, there are several big pores which are located at the periphery and the central zones of the specimen (see Figs. ). Furthermore, the highest density of grains at the central layers of the specimen is also a main factor to sustain the settlement of the top layers. This is conrmed by the smallest void ratio of S11G in Fig. 2.17 and a band in specimen of the scans S10G and S11G (see Fig. 2.12).

Conclusion

In this chapter, we have performed two observations of wetted glass beads: the rst one is the 1D compression test, the second one is the XRCT test combined with in-situ compression test. Two dierent scales of scans were performed in order to quantify the degree of homogeneity and to observe the microstructural changes of specimen under growing of compaction levels.

Firstly, the results show that a very loose yet stable state of wet granular soils can be made in laboratory for tests, even using very small water content. This is one the most interesting aspects of this work as compared to the past studies which were focused on denser states and/or higher water contents (e.g. Bruchon et al. [START_REF] Bruchon | Full 3d investigation and characterisation of capillary collapse of a loose unsaturated sand using x-ray ct[END_REF]). Evidently, the structure of the specimen is heterogeneous but it seems very dicult to make the specimen more homogeneous with so low densities.

Secondly, the macroscopic behavior of wet beads at loose state obtained from oedometric compression tests has usually three regimes: (i) at stresses lower than the preconsoldiation stress, the initial structure can sustain these stresses with a very small rearrangement; (ii) when the stresses exceed the preconsoldiation stress, the loose structures collapse and restructure; and (iii) at higher stresses, the elastic deformation of grains at contact points occurs.

Thirdly, the XRCT investigations of the microstructure allowed quantifying the degree of homogeneity, the distribution of porosity in the specimen as well as the pore-scale behavior under growing of compaction level.

Finally, a new approach for detecting the spherical structures from tomography images is proposed in order to characterize the microstructures of loose systems. Several fundamental microstructural properties are obtained such as total coordination number, number of close neighbors, and distribution of number of contacts, contact network, and radial distribution function. Although our results show an underestimated prediction of several microstructural properties for the other studies, the method presented in this work provides an interesting characterization of the grain-scale compression behavior of wet granular soils (with type of spherical shape) at loose state.

CHAPTER 3 Isotropic compression: DEM study

This chapter reports on a 3D numerical study, using the discrete element method, to investigate the quasistatic response of very loose assemblies of frictional spherical grains to an isotropic compression. The grains are mixed with a small amount of an interstitial wetting liquid, which forms capillary menisci causing an attractive intergranular forces. Under growing pressure, the inuence of the initial assembling process and of various micromechanical parameters on the plastic compression curves are rstly studied. The plastic responses along those compression curves are then shown in particular without and with the appearance of rolling and pivoting resistances in contacts and also the eect of size polydispersity. The evolution of microstructure and force transmission along the compression curve are nally characterized.

Introduction

This chapter presents the results of a discrete numerical simulation study of wet beads under isotropic compression. Capillary cohesion due to the liquid bridges beween neighboring grains stabilizes very loose structures which tend to collapse under growing pressure. The material behavior in compression is thus characterized by strongly irreversible plastic eects, which are investigated at the macroscopic and microscopic levels. This study is similar to the one carried out by Gilabert et al. [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF][START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF] in 2D. However, the numerical model is designed to quantitatively simulate a real material and enables direct comparisons with laboratory experiments.

The model material, which is the same as in the numerical simulations of Khamseh et al. [START_REF] Khamseh | Flow of wet granular materials: a numerical study[END_REF], is described in Sec. 3.2. Some basic properties and important relations in the DEM approach to the mechanical behavior of granular materials are recalled in Sec. 3.3. The numerical procedure and the choice of parameters use to simulate isotropic compression are presented in Sec. 3.4. Then the results both for the macroscopic behavior and the microstructure evolution are given in Sec. 3.5 in which the inuence of various parameters is investigated. The material in which the contacts are endowed with resistance to rolling and pivoting (RPR), a micromechanical feature which is known to strongly aect the mechanical behavior, is studied in Sec. 3.6.

Model material and interaction laws

A granular material is dened by the grain geometry and the intergranular forces. The contact law in a granular material is the relationship between the relative motion of two contacting grains and the contact force. As we deal with particles that may attract one another at short distance without touching, the law governing intergranular forces and moments is best referred to simply as the interaction law. Such elements of the numerical model, designed to simulate wet solid beads, are given here.

Particle shape and size distribution

Eects of angularity on internal friction are known to be important, and the morphology of the liquid phase is also expected to be strongly aected by the geometry of the solid grains. Nevertheless this study focusses on spherical grains, as the simplest model material for which laboratory studies are possible.

One possible way to represent non-spherical particle shape is to implement rolling and pivoting resistance (RPR), as dealt with in Sec. 3.6, Ref. [START_REF] Estrada | Identication of rolling resistance as a shape parameter in sheared granular media[END_REF] established a correspondence between rolling resistance of spherical grains and a certain shape parameter for non-spherical grains. The material behavior reported in Sec. 3.6 might thus be expected to be similar to the one of assemblies of grains departing from the spherical shape, although the change in the geometry of liquid menisci is not accounted for.

Ẽ = E/1 -ν 2 and d ij = 2d i d j /(d i + d j )
, the Hertz law, as implemented in a number of previous numerical studies [START_REF] Agnolin | Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks[END_REF][START_REF] Khamseh | Flow of wet granular materials: a numerical study[END_REF], adequately describes the normal elastic force, F eN ij , in their contact, depending on deection h ij , as

F N,e ij = Ẽ 3 d ij |h ij | 3/2 . (3.4)
This corresponds to a deection-dependent (or force-dependent) normal stiness

K ij N = dF eN ij d|h ij | = Ẽ 2 d ij |h ij | 1/2 .
d i = d j = d ij = d)
or polydisperse ones with the diameter distribution given in (3.1).

Viscous damping

Just like in simulations of dry granular materials [START_REF] Agnolin | Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks[END_REF][START_REF] Peyneau | Frictionless bead packs have macroscopic friction, but no dilatancy[END_REF][START_REF] Peyneau | Solidlike behavior and anisotropy in rigid frictionless bead assemblies[END_REF] we introduce a normal viscous force in intergranular contacts. The viscous normal force opposes the relative normal velocity δV ij N = ḣij of contacting beads, and is chosen as a constant fraction ζ of the critical damping coecient in a linear oscillator with mass m * ij /2 (dening m * ij as the harmonic mean 2m i m j /(m i + m j ) of masses m i and m j of the contacting grains), and stiness constant K ij N as written in (3.5):

F ij N,v = 2ζ m * ij K ij N δV ij N = ζ(m * ij Ẽ) 1/2 d ij h 1/4 ij δV N . (3.6)
This model corresponds to a velocity-independent restitution coecient in a binary collision [START_REF] Peyneau | Frictionless bead packs have macroscopic friction, but no dilatancy[END_REF]. Such a form of the viscous force is more justied by its computational convenience than by its physical relevance. In particular, it is not related to the viscosity of the liquid forming the menisci. The viscous forces, anyway, are quite irrelevant for the quasistatic mechanical properties of the wet bead assemblies, which are studied here. Their role, in the simulation practice, is to ease the obtaining of equilibrium states by dissipating the kinetic energy.

Capillary forces

The appearance of a small quantity of an interstitial wetting liquid introduces additional capillary forces, transmitted between contacting or neighboring grains by a liquid bridge, or meniscus, as depicted in the inset gure of Fig. 3.2. Our work is limited to the pendular state of isolated menisci [START_REF] Mitarai | Wet granular materials[END_REF], and the capillary forces are pairwise additive, acting on grains joined by such menisci. These forces are only present if grains have been in contact [START_REF] Kohonen | On capillary bridges in wet granular materials[END_REF] and the meniscus has not broken since. A meniscus of volume V m remains until the separation distance reaches its rupture value D r V 1/3 m , as a result of calculations by Willett et al. [START_REF] Willett | Capillary bridges between two spherical bodies[END_REF], Maeda et al. [START_REF] Maeda | Evaporation and instabilities of microscopic capillary bridges[END_REF], and Herminghaus et al. [START_REF] Herminghaus | Dynamics of wet granular matter[END_REF].

In the case of monodisperse grains, the Maugis approximation [START_REF] Maugis | Contact, Adhesion and Rupture of Elastic Solids[END_REF] is chosen to calculate the attractive force which varies as a function of the distance h between particle surfaces: Two particles moving away from each other after a collision will only separate if their receding relative velocity is large enough to overcome the capillary attraction. The condition that the kinetic energy of the relative normal motion should be larger than the potential energy of the capillary force, i.e. sets a threshold relative velocity V th proportional to V * [13], with

F cap =                        -F 0 h ≤ 0 -F 0 1 - 1 1 + 4V m πdh 2 0 < h ≤ D r (unbroken meniscus since last contact) 0 0 < h ≤ D r (no
V * = F 0 D r m . (3.8)
More precisely, one has

- Dr 0 F cap (h)dh = F 0 Dr 0 1 - 1 1 + 4V m πdh 2
dh.

(3.9)

V * is thus equal

V * = kF 0 V 1/3 m m , (3.10) with D r V 1/3 m
and k is a factor from the result of Eq. 3.9, as follows

k = 1 + 4V 1/3 m πd -1 + 4V 1/3 m πd . (3.11)
For polydisperse beads we resort to the Soulié formula [START_REF] Soulié | Inuence of liquid bridges on the mechanical behaviour of polydisperse granular materials[END_REF], which is a t to a numerical computation of the capillary force between particles with dierent diameters.

Eq. (1.7) in Sec. 1.3.3.1 of Chapter 1 is rewritten, for two beads of radii R 1 and R 2 , as

F cap = F 0 R 1 R 2 c + exp a h R 2 + b , (3.12) with          a = -1.1 V * m -0.53 b = -0.148 ln V * m -0.96 θ 2 -0.0082 ln V * m + 0.48 c = 0.0018 ln V * m + 0.078, (3.13)
in which V * m is the dimensionless ratio of the meniscus volume (proportional to the surface area of the grain pair) to the cube of the maximum R max = max(R 1 , R 2 ):

V * m = V m R 3 max = V (0) m (R 2 i + R 2 j ) 2R 2 min . (3.14)
In all our simulations the volume meniscus is attributed a xed value V m : a meniscus is assumed to form, from the liquid layer coating isolated grains, as soon as a contact appears. The assumption of a constant volume, while not really justied on physical grounds, enables calculations to be carried out without having to model the liquid transport. While seemingly a gross oversimplication, it will be shown to have negligible inuence on the results. In polydisperse systems, V m varies according to the radii of the grains, and is proportional to the available surface area in the pair, as implied in (3.14), in which V (0) m denotes the meniscus volume for two beads with the minimum diameter of the distribution.

to the normal force F N ij in the right-hand-side. As shown in Fig. 3.3, this entails that the Coulomb cone corresponding to total static contact forces in the space spanned by force components F N , F T is shifted to the left , with its tip at F N = -F 0 , the maximum tensile force. As a consequence, in the contact between two grains of an isolated pair in equilibrium, the total normal force vanishes, but a sliding motion might be opposed by a tangential force as large as µF 0 .

In all the numerical computations of the present study the friction coecient is set to µ = 0.3.

Resistance to rolling and pivoting (RPR)

The inuence of rolling and pivoting resistances (RPR) at contacts, as investigated in 2D by Gilabert et al. [START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF], is also studied. For simplicity, we only implement this feature with linear contact elasticity. Therefore, the linear elastic repulsive force, F N,e , is dened as:

F N,e = K N |h|, (3.18)
where |h| is the normal deection in the contact as the spheres are pressed against each other, K N is the normal stiness coecient, related to the elastic moduli of the material that the grains are made of. In the present study, K N is chosen such that K N d = 4.10 8 Pa.m. Likewise, tangential stiness K T is deection-independent in the linear elastic contact model, and chosen here equal to K N .

The moment that sphere i exerts onto its contacting neighbour j, of radius R j , at its centre, is denoted as Γ ij . It is rst due to the tangential contact force, then to a possible moment Γ r ij of the force density distribution within the contact region, in its centre. One writes

Γ ij = F T ij × R j n ij + Γ r ij . (3.19)
Γ r ij is most often neglected on dealing with smooth, convex particle shapes, because the contact region is very small on the scale of the particle radius.

However, non-negligible moments at the contact point Γ r ij might result from particle surface roughness, as contact regions become larger than the ones deduced from contact elasticity. Separating tangential and normal components of Γ r , one gets the rolling torque, Γ R (a tangential vector), and the pivoting torque, scalar Γ P = n • Γ r . Those moment components respectively work in rolling and pivoting relative motions. The rolling velocity δΩ R ij in contact ij is the tangential part of the dierence in rotation velocities, δΩ ij = Ω i -Ω j , while the pivoting velocity δΩ P ij is its coordinate along normal vector n ij .

The contact law for contact moments, expressing RPR, is similar to the tangential contact law, it involves spring constants K R , K P and Coulomb-like thresholds with friction coecients µ R , µ P . Just like the tangent elastic force, Γ R and Γ P are incrementally 106 Chapter 3. Isotropic compression: DEM study updated in the bead motions, as ∆Γ R = K R δΩ R ∆t ∆Γ P = K P δΩ P ∆t (3.20) in each time step ∆t, while the following Coulomb inequalities are enforced whenever necessary by reducing the norm of vector Γ R and the value of scalar Γ P :

||Γ R || ≤ µ R F N,e |Γ P | ≤ µ P F N,e . (3.21)
Friction coecients µ R and µ P have the dimension of a length. They are of the order of magnitude of the lateral dimension of the contact region, due to surface irregularities.

In practice, they should be equal to a small fraction of the particle radii.

Stiness constants K R , K P have the dimension of a force multiplied by a length, or of an ordinary linear spring stiness multiplied by the square of a length. K R and K P should be such that equality in conditions (3.21) should be reached for small relative rotation angles (otherwise the integration of rotations should not be carried out as vectors).

Values of order µ R RK N (R being the particle radius) are physically reasonable in this respect.

Just like sliding friction (with Coulomb condition (3.17)), rolling and pivoting friction involve the normal elastic force in the right-hand-side of inequalities (3.21). One could thus represent possible contact moment values versus this total static normal force as a Coulomb cone, with its tip shifted to negative values, at -F 0 , similar to Fig. 3.3. A contact between two beads in which the total normal force vanishes may transmit torque components as large as µ R F 0 in a tangential direction or µ P F 0 in the normal direction.

Basic denitions and properties

Basic aspects of DEM simulations are recalled here, which prove useful in the design of the numerical tests and in the analysis of their results. More precisely, F ij and Γ ij are exerted by i onto j at the contact point. Simulations do not model material deformation in a contact region, but consider overlapping particles.

DEM equations, boundary conditions and stress control

The contact point is dened as the center of the intersecting region of the two spheres. In the case of the capillary interaction between spherical beads without contact, the force will be normal to the surfaces at the points of nearest approach and therefore carried by the line of centers.

Let r i and r j , respectively, denote the position of the centre of spheres i and j, as shown in Fig. 3.4. r ij = r j -r i is the vector joining the centres of i and j and

h ij = | r ij | -(d i + d j )/2
is the overlap distance. The degrees of freedom, in addition to the positions r i are the angles of rotation θ i , velocities v i , angular velocities ω i of the grains (1 ≤ i ≤ N ), the dimensions (L α ) α=1,3 of the cell containing the grains and their time derivatives, through the strain rates εα = -Lα L 0 α , (3.22) in which L 0 α denotes the initial size for the corresponding compression process. The time evolution of the particle degrees of freedom are usually governed by the following equations: Stresses are measured using the standard formula: for all coordinate index pairs α, β,

m i d 2 r i dt 2 = N j=1 F ij I i dω i dt = N j=1 Γ ij .
σ αβ = 1 Ω N i=1 mv α i v β i + i<j F α ij r β ij .
(3.24)

Here, Ω = L 1 L 2 L 3 is the sample volume. Cell lengths (L α ) α=1,3 vary simultaneously with the grain positions. We use stress-controlled dynamics, requesting cell dimensions and particle positions to adjust until a mechanical equilibrium state is obtained in which the diagonal elements of the stress tensor take their prescribed values Σ α , 1 ≤ α ≤ 3.

Equations of motion for cell dimensions L α are written as [START_REF] Agnolin | Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks[END_REF].

M Lα = 1 L α L 2 α i m i r(α) i 2 + i<j F (α) ij r (α) ij - Ω L α Σ α . (3.25)
The right-hand-side of (3.25) is proportional to σ αα -Σ α . M denotes a generalized mass associated with changes in the cell dimensions. It is usually chosen, in compression tests, of the order of the total mass of all particles in the system. Eq. (3.25) requests the cell to expand if the stress (evaluated with the sign convention of soil mechanics: tensile stresses are negative) is too large compared to the prescribed value, and to shrink if it is too small.

In all equations written here, the periodicity of the simulation cell has to be taken into account: particles exiting the cell through any boundary are being replaced by their image (or copy ) entering the cell through the opposite boundary. Sums over interacting pairs involving branch vectors (centre-to-centre vectors) r ij , as in (3.24), are to be understood with the nearest-image convention: vector r ij actually points from the centre of grain i to the centre of the nearest copy of j within the periodic set of replica of grain j.

Equilibrium condition

Although numerical simulations of the quasistatic response of granular materials requires mechanical equilibrium to be reached, equilibrium criteria are sometimes left unspecied or quite vaguely stated in the literature. However, in order to report on usual quantities like the coordination number or the force distribution, it is essential to know which pairs of grains are in contact and which are not. Due to the frequent occurrence of small contact force values, this requires forces to balance with sucient accuracy. We found that the following criteria allowed us to identify the force-carrying structure clearly enough.

Basic denitions and properties 109

A typical intergranular force value F 1 =max(F 0 , P d 2 ) sets the tolerance levels for individual grain equilibrium, where P is the applied pressure. A conguration is deemed equilibrated when the following conditions are simultaneously satised: (i) the net force (the total force) on each spherical grain is less than 10 -4 F 1 ; (ii) the total moment on each sphere is lower than 10 -4 F 1 d; (iii) the dierence between imposed and measured stresses is less than 10 -4 F 1 /d 2 ; and (iv) the kinetic energy per grain is less than 10 -8 F 1 d.

Coordination numbers

In static granular assemblies, the coordination number z is, after solid fraction Φ (which is the volume occupied by the solid grains divided by the total system volume), the most important and well-studied variable characterising the material state. It is dened as the average number of force-carrying contacts per grain. If N c is the number of force-carrying contacts between N grains, then the coordination number is z = 2N c /N as the product zN counts each contact twice. In general disordered systems, the coordination number has an upper bound well below the values it may reach in ideal, regular arrangements like perfect crystals. Thus, disordered packings of spherical beads, in the limit of small contact deections (rigid contact limit), satisfy z ≤ 6 [START_REF] Roux | Geometric origin of mechanical properties of granular materials[END_REF][START_REF] O'hern | Jamming at zero temperature and zero applied stress: The epitome of disorder[END_REF][START_REF] Donev | Pair correlation function characteristics of nearly jammed disordered and ordered hard-sphere packings[END_REF][START_REF] Agnolin | Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks[END_REF]. This upper bound to z is due to the absence, in generically disordered grain packings, of contact networks able to carry self-balanced normal forces (the degree of force indeterminacy or degree of hyperstaticity of the contact network is equal to zero if contacts are regarded as frictionless).

The coordination number z is split into two terms: the average number of contacts per grain, z c , and the average number of the pairs interacting though a capillary bridge per grain, z d , so z = z c + z d . Hence, z is also called the total coordination number.

Furthermore, it is interesting to dene and to measure a coordination number of close neighbors, z(h), as the average number of neighbors of one sphere separated by an interstice ≤ h. z(0) is the usual contact coordination number z c , and z(h) characterizes the packing geometry in terms of position correlations.

A decomposition of the coordination number of contacts z c is:

z c = z + + z -, (3.26)
in which z + is the coordination number of compressive bonds (positive forces) and z -is the coordination number of tensile bonds (negative forces). Note that z + includes only the contacts in which the repulsive elastic force is larger than the capillary force, and z - includes all the distant interacting pairs and also all the contacts in which the capillary force dominates. In wet granular materials, depending on the amount of the liquid in the system and on its resulting spatial distribution, dierent denitions for the liquid content are possible.

Several authors use the liquid volume fraction or V l /Ω, dened as the total volume of liquid, V l , divided by the total sample volume Ω [START_REF] Kohonen | On capillary bridges in wet granular materials[END_REF][START_REF] Herminghaus | Dynamics of wet granular matter[END_REF]. Denoting as V g the volume of all N solid grains in the system, one has

Ω = V g + V v , in which V v
Φ = V g /Ω = 1 -V v /Ω = πN d 3 /(6Ω). The liquid content W is then given by W = V l Ω = 3zΦV m π d 3 .
(3.27)

Other studies (e.g. in [START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF]) dene the liquid content in terms of masses. It is also possible to use degree of saturation S w instead of liquid content [START_REF] Mitarai | Wet granular materials[END_REF][START_REF] Scholtès | Micromechanics of granular materials with capillary eects[END_REF]. The present study, like a number of previous ones [START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF][START_REF] Radjaï | Bond anisotropy and cohesion of wet granular materials[END_REF][START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF], is restricted to the pendular state of low saturation, in which the wetting liquid is conned in bonds or menisci joining grains in contact or separated by a small distance. The degree of saturation is therefore dened as the ratio of the liquid volume, V l , to the void volume V v and is simply related to liquid content as

S w = V l V v = n i V m (1 -Φ)V = 3zV m Φ (1 -Φ)πd 3 .
V m = πa 2 h -2v i = 2πR 3 3
(2 cos 3 ϕ -3 cos 2 ϕ + 1).

(3.29)

The volume of the meniscus is therefore approximately 8.10 -3 d 3 when the lling angle reaches its maximum value ϕ = π/6. Note that, here, the so-called toroidal approximation is used: the curves limiting the meniscus in a plane containing its axis of symmetry are approximated as circular arcs. Formula (3.29) also assumes a contact angle equal to zero. Eq. (3.28) then predicts a maximum saturation between 0.05 and 0.1, similar to experimental observations [START_REF] Herminghaus | Dynamics of wet granular matter[END_REF][START_REF] Mitarai | Wet granular materials[END_REF]. With the selected value of V m /d 3 = 10 -3 , used, as a reference value, in most simulations of this study, Eq. (3.28) predicts the degree of saturation in this work to be of the order of 10 -3 .

Average normal force

Given the formula (3.24) for the stress tensor, it is easy to deduce a convenient expression for the average intergranular normal force [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF][START_REF] Agnolin | Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks[END_REF]. From

P = 1 3Ω i<j F ij • r ij , (3.30)
one gets, r ij = (d + h ij ) n ij , in which n ij denotes the normal unit vector pointing in contact i -j.

Eq. (3.30) can be rewritten as:

P = 1 3Ω i<j (R i + R j )F N ij + 1 3Ω i<j h ij F N ij (3.31)
The second term is usually neglected with grains in contact (as the elastic deection h ij < 0 is very small). For distant interactions through menisci (h ij > 0) it was also observed [START_REF] Khamseh | Flow of wet granular materials: a numerical study[END_REF] to be very small for the range of meniscus volumes corresponding to the pendular state. Replacing the sums with the average values and using the denitions of Φ, z, and z d gives us:

P = Φz πd 2 F N (3.32)
The relation (3.32) can be accurate in all simulated states with or without cohesion [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF].

However, in the presence of cohesion, it cannot be used to predict the typical contact force when the conning pressure is dominated by the attractive forces. Then, contact forces of order F 0 are quite common [START_REF] Richefeu | Stress transmission in wet granular materials[END_REF], while the average force is smaller.

In a moderately polydisperse system, (3.32) becomes

P zΦ d π d 3 F N , (3.33)
provided the correlation between the force value and the diameters is neglected.

/F 0 Interparticle distance, h/d Maugis, V m /d 3 =1.0x10 -3 Maugis, V m /d 3 =7.8x10 -6
Soulie, V m /d 3 =1.0x10 -3

Simulation procedures for the compression test 113

Our simulation procedure, to be detailed below, involves an aggregation stage to form low pressure initial states, in which grains are agitated with a characteristic velocity V 0 .

The microstructure of the resulting aggregates, as the initial kinetic energy is gradually dissipated, depends on whether the initial agitation is strong enough to cause breakage and restructuration of the contact networks, and thus on dimensionless ratio V 0 /V * , with V * dened in Eq. (3.8). This ratio will be observed to have a strong inuence on the initial coordination number.

The third dimensionless group introduced by capillary cohesion, an important control parameter, is the reduced pressure,

P * = d 2 P F 0 , (3.35)
comparing the applied pressure to the tensile strength (maximum attractive force) in contacts. P * 1 signals cohesion-dominated systems, for which attractive forces may stabilize loose structures. Conning forces dominate for P * 1, and attractive forces become negligible the bead assembly behaves as a dry granular material. The simulation parameters are chosen such that P * reaches large values before κ -1 decreases to 10 -3 , and thus the eect of pressure in this study is due to the competition between conning stress and capillary forces, not to compression induced contact deections.

Additional parameters are µ R /d = µ P /d in system with RPR, as well as the rotation angles for which rolling and pivoting friction thresholds are reached (small enough to be irrelevant in our case).

In the polydisperse system, a reduced pressure is redened as

P * po = d 2 P πΓ d . (3.36)
With the chosen diameter distribution (Sec. 3.2.1), one has

P * po = P πΓ d min d max d max -d min ln d max d min . (3.37)
3.4 Simulation procedures for the compression test

Specimen preparation

Before applying a controlled external pressure, we rst assemble disordered, loose congurations in which the grains form connected aggregates spanning the sample.

A preliminary step consists in placing the grains within the cubic, periodic simulation cell, forming a disordered conguration without contacts. We thus choose the initial solid fraction Φ 0 , which may be set to low values (down to 0.2) compared with the typical densities of dry granular systems (Φ hardly decreases below range 0.550.58, depending on the friction coecient, in conned cohesionless bead packs). Several techniques can be used to obtain disordered congurations of non-contacting beads. ).

A convenient one uses the melting of a regular lattice [START_REF] Agnolin | Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks[END_REF]. In this study we use samples of N = 4000 (often monodisperse) beads, initially placed on the sites of a Face-Centered Cubic (FCC) lattice within the cubic cell, as shown in Fig. 3.7. In order to obtain a disordered system, the beads are rst regarded as hard objects interacting with collisions that preserve kinetic energy. Once launched with random velocities, the initial order is quickly forgotten (provided Φ does not exceed the crystallization value 0.49 for monosized beads). This initial mixing stage is a purely numerical procedure in which the mechanical properties of the particles do not conform to the physical model.

Next comes an aggregation stage, for which the grains are endowed with all the properties of elasticity, friction (possibly RPR) and capillary cohesion as described in Sec. 3.2.

This numerical procedure is supposed to mimic, in some idealized version, the physical process by which low density equilibrium congurations of wet beads are stabilized. Once attributed random velocities, which are Gaussian distributed with a mean quadratic velocity V 0 , particles collide, and the capillary attraction causes them stick to one another.

A contact network gradually forms, and kinetic energy gets dissipated. The process is stopped when all particles are connected, forming one single cluster joined by liquid bonds, and when all particles reach equilibrium in the sense of Sec. In each compression step, external reduced pressure P * is multiplied by a constant factor 10 1/4 1.7783, and one waits until the new equilibrium conguration is reached, with the criteria stated in Sec. 3.3.2. The main objective of the present work is to study of the eect of a gradual compression, starting from cohesion-dominated loose states at small P * , and ending in connement-dominated denser states at large P * . The compression program is pursued up to P * max = 10 4 (or for some systems, to P * max = 10 3 ). Then, the reverse path of P * back to minimum value P * = 10 -3 from its highest value is also simulated. Based on the range of the reduced pressure P * , three regimes are identied from the compression curve (see the inset of Fig. 3.9(a)). A rst regime (Regime I) is observed for low reduced pressures, P * ≤ 10 -2 , in which the initial structure still sustains the increasing pressure without rearrangement, and void ratio e remains nearly constant (Fig. 3.9(a)). In a second stage (Regime II), roughly corresponding to interval 2×10 -2 ≤ P * ≤ 2×10 0 , the void ratio strongly decreases, and might be described as linearly varying with log P * . The loose structures formed at low P * are no longer able to support the increasing conning pressure, they collapse and restructure. Finally, in Regime III, e gradually approaches some minimum void ratio e min at the end of the loading process.

Φ 0 V 0 /V * V m /d 3 µ R /d = µ P /

Material behavior in isotropic compression

The slight decrease of e for P * ≥ 10 2 is similar to the behavior of cohesionless systems, and due to elastic deections in a stable contact network. The chosen parameters are such that κ remains large enough not to inuence the compression process taking place in Regime II.

The plastic compression behavior of the wet material is closely similar to the 2D results of Gilabert et al. [START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF], and to the experimental compression test results on metallic powders by Poquillon et al. [START_REF] Poquillon | Cold compaction of iron powders -relations between powder morphology and mechanical properties: Part i: Powder preparation and compaction[END_REF]. The void ratio curve in Regime II might be represented with a linear variation with log P * , assuming e ref is the void ratio for some reference Eq. 3.38 is usually used in soil mechanics for cohesive materials [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF].

Upon decompression, e increases slightly, remaining very close to e min : the compaction is essentially irreversible, unlike in the cohesionless granular assembly of Fig. 3.9(b), for which loading and unloading branches are not distinguishable (as observed in the numerical study of compression cycles of cohesionless spherical beads reported in Ref. [START_REF] Agnolin | Internal states of model isotropic granular packings. ii. compression and pressure cycles[END_REF]).

Note the use of a physical unit of pressure in Fig. 3.9(b), corresponding to glass beads of

#cycle

Path 1 A ⇒ E ⇒ E′ 2 A ⇒ B ⇒ B′ 3 A ⇒ C ⇒ C′ 4 A ⇒ D ⇒ D′ 5 C′ ⇒ C ⇒ E ⇒ E′ A B B′ C C′ D D′ E E′ Void ratio, e (-)
Reduced pressure, P * (-) Figure 3.10: Eect of dierent (isotropic) unloading and reloading histories on void ratio; the system does not rearrange along unloading paths BB', CC', DD', EE', which are reversible; Path 5 causes plastic response in Section CE, along which pressure increases beyond its past maximum; the primary curves (Path 1) is then retraced. diameter 118 µm, perfect wetting, and the surface tension of water. With such a system P * = 1 is obtained for P = 2 kPa.

Upon unloading and reloading at various pressure levels along the compression curve, a plastic response (irreversible structural rearrangement with density change) under isotropic pressure is only observed (as illustrated in Fig. 3.10) if the maximum pressure the system has been subjected to in the past (the overconsolidation pressure of soil mechanics) is exceeded. This maximum pressure value appears to fully characterize the history dependence of the system in isotropic compression [START_REF] Than | Basic mechanical properties of wet granular materials: A DEM study[END_REF].

Contact networks and force transmission

As dened in Sec. 3.3.3, an important microstructural characteristic, the coordination number z, dened as the average number of interactions per grain, is the sum of the contact coordination number z c and the coordination number of distant interactions, through menisci joining non-contacting grains, z d . and velocity indeterminacy or hyperstaticity and hypostaticity is not excluded, see Ref. [START_REF] Agnolin | Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks[END_REF]). z c and z d remain unchanged in Regime I (P * < 10 -2 ), as the initial structure is maintained. Both coordination numbers start to increase as the structure collapses and reorganizes in Regime II (10 -2 < P * ≤ 2 × 10 0 ). z exhibits little change in Regime III (P * > 2 × 10 0 ), as the increase of the number of contacts is then mainly due to the closing of narrow gaps between pairs of grains joined by a meniscus, as the structure is further compressed a moderate eect. Upon unloading, some contacts are opened and therefore z d rst increases while z c decreases a little. Then both coordination numbers tend to remain constant till to the nal state at low P * , reecting the stability of the dense structure formed at high P * .

The coordination number of contacts z c may be split into the contribution of compressive bonds, z + , and the one of tensile bonds, z -. Coordination numbers z c , z + , and 122 Chapter 3. Isotropic compression: DEM study z -are plotted versus reduced pressure P * in Fig. 3.11(b). At low P * , both compressive and tensile contact forces are present, with similar values of z + (slightly above 2) and z - (slightly below 2). Nearly constant in regime I, they vary in opposite directions under growing P * in regime II, with more and more compressive forces. At high pressure P * , z -gradually vanishes, as the system essentially behaves like a pack of compressed cohesionless beads. z + and z -vary in opposite ways upon reducing P * , without returning to their initial values. Line thicknesses for wet contacts (red lines), and for distant interaction (blue lines) are proportional to the normal force intensity |F N |.

Contact networks may be visualized in order to provide insight on the structural changes and the evolutions of z c and z d under growing P * . Fig. 3.12 shows the networks of wet contacts and of distant interactions through menisci at six typical stages of P * in the compression cycle. A slice is extracted from the sample with width 3d, its height and length being those of the whole cuboidal simulation cell. In Fig. 3.12, the red lines denote the wet contacts (the total normal force F N = 0), and the blue lines denote the interactions through menisci, without direct contact. Along the loading path, under growing of pressure P * , more distant interacting pairs are created, while the network density is increased (Figs. Similarly, compressive and tensile bonds are also depicted in Fig. 3.13. The contact networks of the z + and z -are selected at several values of P * in compression path (Fig. 3. 13(a, b, c, andd)) and decompression path (Fig. 3.13(e and d)). Initially, at the beginning of the compaction curve, the density of compressive bonds is higher than one of tensile bonds. There also exist several pairs having a null normal force (the magenta lines). This means the normal elastic force F e N is balanced with the capillary force F cap . Upon growing of P * , the compressive and tensile bonds are oppositely changed. While the compressive contacts are dramatically increased, the tensile ones are sharply decreased. We can obviously see the increase in the density of the read lines in Figs. 3.13(b andc). Until the maximal pressure (P * = 10 4 ), the compressive ones reach to the maximum (at z + = 5.75, approximately 99 % of the contacts) while the tensile ones reduce to minimum corresponding to z -= 0.04 (approximately 1 % of the contacts). For the unloading cycle (see Figs. force F N at larger reduced pressure P * . The normal force distributions are roughly symmetric about zero in the initial states and under low P * , (see e.g., the red and dark green curves). At the beginning of Regime I, P * = 10 -3 , all of the total normal forces F N are very small. This is due to the near lack of force indeterminacy: if the network were exactly isostatic, all forces would vanish under zero pressure. These P.D.Fs widen and take notable values near F N = -F 0 at the end of Regime I, signalling that more and more contacts are on the verge of opening. At pressure P * = 10 -1 , all P.D.Fs widen on the positive (repulsive) side, while the largest tensile force is xed at F 0 , thus becoming dissymmetric, as shown in Fig. 3.14(a). Given relation (3.32), the average normal force tends to scale with pressure P (up to moderate eects of the variations of contact density zΦ). In the last compression stage (regime III), these P.D.Fs, if normalized by the averge F N , become narrower, as apparent in Fig. 3.14(b). This is similar to the behavior of cohesionless systems, in which those PDF's tend to decrease exponentially, on a scale given by the average F N , the faster the larger the pressure (or the smaller κ) see Void ratio, e (-)

Reduced pressure, P * (-) the rearrangement of microstructure due to the decrease in liquid content under growing pressure P * . For P * = 10 -3 , close to the initial state, which is essentially V m independent (because of the absence of distant interactions), the connectivity is not aected. For larger pressures, the structure gets somewhat depleted (although this is not a strong eect)

F30V1Vm1 F30V1Vm2 F30V1Vm3 F30V1Vm4 F30V1Vm5 F30V1Vm6 F30V1Vm7 F30V1Vm8 F30V1Vm9 F30V1Vm10 (
The eect of V m on P.D.Fs of normal forces is also observed to be very small. Under low pressure P * ≤ 10 -2 , z(h) is almost similar (red and green curves) because hardly any structural change takes place in Regime I, and its increase with h is slow, with a slight upward curvature (convex curve). At P * = 10 -1 , z(h) slightly increases.

This corresponds to a slight increase of the z c in Fig. 3.11(a). For growing P * , function z(h) gradually retrieves the concave shape observed in cohesionless systems [START_REF] Agnolin | Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks[END_REF], with larger contact coordination numbers, and correlatively less incipient contacts at short distance. Beyond a few times d/100, z(h) steadily increases with density.

Within the interparticle the meniscus rupture distance D r , each grain has on average z(D r ) -z c non-contacting neighbors, among which z d neighbors are interacting by an attractive force. The ratio of z d to z(D r ) -z c is the proportion of the neighbors within the range D r that are bonded by a capillary bridge. This ratio for dierent values of P * in loading path of the reference case are given in Table 3 These values are similar to the ones of Khamseh et al. [START_REF] Khamseh | Flow of wet granular materials: a numerical study[END_REF] and to the proportion of about 50% reported by Kohonen et al. [START_REF] Kohonen | On capillary bridges in wet granular materials[END_REF] in static grain packs. In practice, one may act on a wet system by changing its saturation. The most drastic change should be obtained on entirely suppressing the capillary cohesion, either by drying, or by completely saturating the intergranular voids by the liquid. In numerical simulations, one may simply remove all capillary forces, leaving only the interactions present in a dry system. It is interesting to observe the eects of such an ideal transformation, carried out at various points along the irreversible compression curve. Fig. 3.20 illustrates the resulting void ratio evolution, if the system is suddenly deprived of capillary forces at dierent pressure levels. This ideal drying or saturation step produces a sudden collapse (a brutal compression step), unless all irreversible compression has already taken place, as for Points F and G in Fig. 3.20. In such cases, one may remove all capillary forces, as their mechanical role, at high P , is negligible.

F45V1Vm3, A ⇒ B ⇒ B′ C ⇒ C′ D ⇒ D′ E ⇒ E′ F ⇒ F′ G ⇒ G′ Non-cohesive case
A more gradual collapse due to progressive imbibition is reported in some experiments [START_REF] Bruchon | Full 3d investigation and characterisation of capillary collapse of a loose unsaturated sand using x-ray ct[END_REF].

Remarkably, the nal state after decompression keeps the same density, whether or not the system has been deprived of capillary cohesion [START_REF] Than | Basic mechanical properties of wet granular materials: A DEM study[END_REF].

Inuence of initial state

The assembling conditions notably aect the packing structure and microstructure under low pressure P * . It is therefore necessary to assess the role of initial congurations, which are produced by the initial assembling process in the presence of with capillary bonds.

The eects of various micromechanical parameters on the plastic compression curve are also discussed. The eects of both solid fraction Φ 0 and initial agitation intensity V 0 /V * are investigated.

Inuence of initial solid fraction

Among the features aected by the assembling process, the competition between compression and aggregation is the most important one. Fig. 3.21, obtained with standard values V 0 /V * = 0.2041, V m /d 3 = 10 -3 , compares the compression curves of specimens with dierent initial solid fractions Φ 0 . Denser systems are able to support larger pressures before rearranging, whence larger Regime I plateaus (see Fig. 3.21(a)) in dierent intervals of reduced pressure P * . For instance, from 10 -3 ≤ P * ≤ 4 × 10 -3 for the smallest Φ 0 value (F25V1Vm3) to 10 -3 ≤ P * ≤ 8 × 10 -2 for the largest one (F45V1Vm3).

However, in Regime II, these compression curves tend to converge and they very nearly coincide at the end of the collapse (for P * = 2 × 10 0 to P * = 5 × 10 0 ), as clearly apparent in Fig. 3.21(a), with the dierent slopes of the six curves corresponding to dierent Φ 0 .

In Regime III (P * > 5 × 10 0 ) and along the unloading path, the void ratio of the initially denser systems is only very slightly lower than the one of the looser systems.

Furthermore, the evolution of coordination numbers of z c and z d remain very similar during the loading and unloading paths whatever the initial solid fraction Φ 0 (see Fig. 3.21(b)). z c is only very slightly dierent in Regime II, depending on the moment when the collapse started. Meanwhile, in the second half of Regime II, the dierence in z d is also quite insignicant, with slightly less contact openings leading to distant interactions in the denser systems. Along the unloading path, the denser the systems, the lower z d too but this dierence is also quite small. One essential conclusion of these observations is that the coordination numbers depend much more on the system history (assembling process, loading path), than on its current density.

Looking at the coordination numbers of compressive and tensile bonds (Fig. one also observes quite insignicant dierences for dierent initial solid fractions. z + and z -at the beginning of loading curves (P * < 10 -1 ), despite the dierent conguration geometry. z + increases to extent, while z -decreases, for larger initial solid fraction Φ 0 (see clearly in the inset). Under growing pressure, z + and z -evolve in opposite directions, similarly for all compression curves. System states appear no longer to be inuenced at all by the initial density in Regime III and in unloading. Then, like other state variables, connectivities signal no dierence in Regime III and upon unloading.

Eects of polydispersity

Comparison between polydisperse and monodisperse systems

Previous results show that the introduction of polydispersity does not very notably change the observed behavior and microstructure evolution of initially loose systems subjected to a pressure cycle. We now turn to more quantitative comparisons. The monodisperse case with the initial solid fraction Φ 0 = 0.25 is selected to compare with the polydisperse one. These two cases are compared with the same values of initial agitation intensity and meniscus volume. One may note that both quantities should perhaps be adequately redened, accounting for the dierent grain scale: the meniscus volume, as mentioned here, is actually the value of V m corresponding to a pair of particles of equal diameter d min , while the meniscus volume is actually proportional to the surface area of the particle pair; and the value of V * should in principle be modied too, accounting for size dependences of the mass and of the meniscus rupture distance.

However, reference parameter values pertain to the limit of gentle aggregation, without cluster restructuration, and to a range of meniscus volumes with very little inuence on the behavior, and, consequently, we dispense with rescaling V 0 and V m . And the minimum void ratio e min is slightly smaller for the PDS the solid fraction of a polydisperse system is usually larger than for monosized grains.

The evolutions of coordination numbers z c and z d for both systems, plotted in Fig. 3.33(b), are also very close. At the end of the compression z d is somewhat smaller in the PDS this could be due to the dierence in ratios D r /d. In nal dense states after decompression, small dierences are maintained for z c . Generally, given the geometric dierences between monodisperse bead packs, and polydisperse ones with diameters varying by a factor of 2, it should be concluded that the observed behaviors do not exhibit very signicant size distribution eects, and thus the study of a mondisperse model 3.6. The case of wet beads with RPR 143 gives meaningful results which might be applicable in rather generic situations involving spherical beads with moderate polydispersity.

The case of wet beads with RPR

Rolling resistance was observed in two-dimensional [START_REF] Kadau | Pore stabilization in cohesive granular systems[END_REF][START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF][START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF] and three-dimensional [START_REF] Bartels | The eect of contact torques on porosity of cohesive powders[END_REF] DEM studies to have strong eects on loose cohesive granular materials for small P * .

We now report on the results of the isotropic compression cycle applied to the systems of Table [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF][START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF]. In all simulations with RPR, we set K R /(K N d 2 ) = 2.5 × 10 -2 , while K T /K N = 1 and K N /d = 400 GPa. The compression cycle is limited to P * ≤ 10 3 . Fig. 3.34 shows the resulting behavior, in terms of void ratio and coordination numbers.

Inuence of initial solid fraction

A rst observation is, in general, the absence of a plateau for Regime I, except in the densest case Φ 0 = 0.45. The irreversible density increase gradually takes place from the lowest investigated values of P * , (of order 10 -3 ). Under growing pressure P * , curves pertaining to dierent initial densities converge, until they nearly coincide in Regime III.

The slope of the linear part of the e versus log P * curve is no more an intrinsic material parameter as in the material without RPR, as it varies with Φ 0 . The nal (minimum) value of the void ratio, reached for large enough P * , and (approximately) maintained upon reducing the pressure, is also dependent on Φ 0 : one gets e min = 0.75 for Φ 0 = 0.40 and e min = 0.70 for Φ 0 = 0.45 (see Fig. 3.34(a)). The coordination numbers of contacts are still equal to the minimum value, which is 2 instead of 4 for systems with RPR: as noted in 2D systems [START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF], aggregates of particles with rolling resistance are rigid objects, and if they do not restructure due to the collisions while additional particles stick to the growing clusters in the formation stage, as in the case of a small initial agitation velocity V 0 , they should be generically devoid of loops, whence z c = 2. Contact loops entails force indeterminacy. The initial aggregate is an isostatic structure in the limit of small V 0 . Fig. 3.34(b) shows that this applies to V1 (small initial V 0 ) systems whatever the value of Φ 0 , and the variation of coordination numbers in the pressure cycle are remarkably independent on Φ 0 . Distant interaction coordination numbers, initially equal to zero, due to the instantaneous stability of aggregated structures, vary a little with Φ 0 in the later stages of the pressure cycle, with slightly less contact openings in denser systems. As noted before Regime I does not exist as a plateau in the compression curves for small initial agitation level (this applies to two rst cases red and green curves in Fig. 3.36). The initial tenuous structures (with coordination number approaching 2, i.e., virtually no loop) formed with RPR and low V 0 respond with some amount of plastic collapse at the rst pressure increment, as shown in Fig. 3.36(a). Regime I reappears on the pressure scale of Fig. 3.36(a) in systems assembled with larger initial agitation velocities, which are better coordinated (see Fig. 3.36(b)), and collapse more abruptly in plastic compression. This large change in the compression curve is due to strong variations of the initial values of coordination number z c , which nearly doubles, from approximately 2.15 to nearly 4.25, for increasing initial agitation. Interestingly, the system assembled with the lower coordination numbers rst react by plastic compression with no contact loss, as z d remains equal to zero (no contact opening creating a distant interaction). As in the 2D systems of Ref. [START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF], irrecoverable plastic restructuring may take the form of bending or folding deformation of single particle strands.

Distant interactions are present from the beginning in better coordinated, initially more strongly agitated systems, with z d = 0.25 for the largest V 0 . At intermediate pressures 1 ≤ P * ≤ 10, while z c slightly increases, z d increases faster, before slightly decreasing in Regime III. z c then reaches its maximum at the end of the loading path.

Upon unloading, many distant interactions replace interparticle contacts. Finally z c and z d remain constant when P * is smaller than 10. to vanish with the degree of force indeterminacy (hyperstaticity), for decreasing V 0 . The positive forces between grains dominate the maximum attractive forces at the maximum velocity (the cyan and brown curves). Like other cases of cohesive and non-cohesive grains, the stronger the pressure, the wider the P.D.Fs and the inuence of initial velocities to force distributions then fade out after the granular structure are compressed to higher pressures and densities.

Inuence of rolling and pivoting friction coecients

DIerent values of the rolling and pivoting friction coecients, µ R /d and µ P /d, are also tested in the same systems. The compression curves for four values of µ R = µ P (see Table 3.3) ranging from 0.004d (R1) to 0.5d (R4) are shown in Fig. 3.41.

As regards Regime I, for the evolution of void ratio, the eect of larger RPR coecients is somewhat similar to the one of a larger V 0 , creating more stable initial congurations. That does not exist in the systems of low µ R /d (R1 and R2). Coordination numbers, however, tend to decrease for increasing RPR, with initial states closer and closer to the isostatic, loopless limit z c = 2. The plastic compression of systems with larger RPR tends to be slower, both for void ratio and coordination numbers. As the irreversible compression proceeds, the dierences between small and large RPR persist Some further analysis of the plastic mechanisms by which network break or rearrange as the system is irreversibly compressed would be necessary to explain the relations between the evolutions of density, coordination numbers, position correlations and force distributions. It was suggested in Ref. [START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF] that the dominant plastic mechanism transits from tensile rupture to rolling friction mobilization and plastic bending of particle chains as the level of RPR increases. Firstly, the irreversible, plastic compression process, controlled by the reduced pressure, P * , goes through three regimes. In Regime I, which is sensitive to the assembling procedure, the microstructure can sustain small, nite pressure increments with strong enough without plastic collapse. Those increments are larger for denser initial states and for stronger initial structures, as obtained as agitation levels in the aggregation process entails cluster restructuration and reinforcement. Subsequent collapse as the conning pressure is increased and causes tensile contact rupture reduce the void volume in Regime II, for which the variation of void index with pressure has the same shape as reported in the soil mechanic literature. The preconsolidation pressure is a plastic threshold below which the quasistatic response of material is approximately elastic (similar to the behavior of a non-cohesive granular soils under isotropic load). The material eventually approaches a minimum void ratio (maximum density) in Regime III.

Secondly, our results emphasize the important dierences between cohesive and noncohesive granular assembles. The most important dierences, which are brought by cohesion in macroscopic behavior, are the existence of very loose stable structures (with solid fraction as low as 0.25 with identical spherical beads) and the plastic compression phenomenon. At the microscopic level, the attractive and repulsive contact forces are initially of the order of the maximum tensile strength F 0 , and tend to compensate each other under low pressure. The compaction path stabilizes self-balanced force networks with large compression forces.

Thirdly, our results show the remarkable, and somewhat unexpected, inuence of two microscopic features. Both strongly aect the shape of the compression curve. One is a material characteristic, the rolling and pivoting resistances at contacts between grains.

Even a small RPR has a strong microstructural role for small P * . The other one is the initial coordination number and internal state of connectivity, as determined by the assembling process. Coordination and connectivity are quite independent of density, and also crucially inuence the plasticity of loose materials.

Finally, the plastic behavior in compression is not very sensitive to the grain size distribution, as only a very moderate quantitative dierence is observed between singlesized grains and polydisperse ones with a diameter distribution ranging from 1 to 2.

CHAPTER 4

Oedometric compression: DEM study and experimental confrontation This chapter reports on numerical simulations of oedometric compression tests with the model material (with monodisperse beads and no RPR). Some specic analyses are introduced to investigate stress anisotropy and its microscopic origins. Such simulations enable direct comparisons with experiments, both for the plastic behavior (compression curve) and for some microstructural characteristics. The objective of this chapter is to investigate the behavior of the numerical material in oedometric compression, and then to quantitatively confront experiment and DEM simulation at macro-and microscopic levels, using the experimental results of Chapter 2.

The DEM simulations are carried out with the simplest model material, the monodisperse one without RPR. The force model and the initial assembling procedure are the same as the ones dened in Chapter 3. Several results of the isotropic compression tests are also reused for the comparison.

Simulation of oedometric compression tests

To the same initial samples of 4000 mono-sized spherical beads as for isotropic compression tests, stress and strain conditions are applied in order to simulate oedometric compression: while σ 1 (the axial stress) is increased, ε 2 = ε 3 , the lateral strains are maintained at zero. The axial strain rate is also constrained, during the compression, not to exceed a preset maximum value, in order to model a quasistatic test. This constraint is expressed through the inertial parameter I = ε m/dσ 1 where m is the particle mass. I is requested to remain below 10 -3 throughout the compression tests simulation.

Like in the isotropic case, we use dierent values of material parameters in order to explore the parameter dependence of the plastic compression behavior: four values of initial solid fraction (Φ 0 ), dierent values of initial agitation intensity V 0 /V * for the reference initial solid fraction Φ 0 = 0.30. As to the meniscus volume, one conclusion of Chapter 3 is its remarkably weak inuence on the material behavior, and therefore it is xed to V m /d 3 = 10 -3 in all oedometric test simulations. The list of numerical samples subjected to oedometric compression, with the corresponding parameter choice, is given in Table 4.1.

Φ 0 V 0 /V * V m /d
σ * 1 = d 2 σ 1 F 0 , (4.1)
The loading program, with controlled steps of externally applied axial stress, is similar to the isotropic compression cycles of Chapter 3. The coecient of lateral pressure, denoted as K 0 , is dened in oedometric compression as the ratio of lateral stress, to axial stress:

K 0 = σ 3 σ 1 . (4.2)
K 0 , in Fig. 4.3, is plotted on a logarithmic scale in order to observe clearly its variation along the compression and decompression curves. K 0 strongly depends on the initial state: it is about twice as large, under low σ * 1 , in the stronger systems (V4 to V6) than in the weaker ones (V1 to V3). This might be attributed to the initial isotropic organization of self-balanced forces of order F 0 , which cause the strong systems to transit less anisotropic stresses in Regime I. K 0 decreases slightly in Regime I, as the anisotropy of the loading process aects the initially isotropic microstructure more and more. Conversely, the collapse taking place in Regime II causes K 0 to gradually increase in all systems. K 0 values converge near the end of the collapse stage, at P * = 10 (4.

3)

The contribution of distant interactions to σ 22 is illustrated in Fig. Upon growing of vertical stress, the system is compressed and more distant interacting pairs are created. σ d 22 is therefore decreased while σ 22 is remarkably increased. This is why σ d 22 /σ 22 values signicantly decrease in the range of vertical stress 2 × 10 -2 ≤ σ * 1 ≤ 3 × 10 -1 for rst three cases (red, dark-green and blue curves) and in the shorter range of σ * 1 in the last three cases (magenta, cyan and dark-orange curves). When the conning forces dominate the attractive forces (σ * 1 > 10 0 ), this contribution is gradually increased until the end of loading path. Upon unloading, σ d 22 is steadily decreased until the end. This corresponds to the increase in the distant interactions in Fig. 4.1(b).

Another decomposition of total stress is given as follows:

σ αβ = σ cap αβ + σ N e αβ + σ T αβ , (4.4)
in which σ cap αβ is the contribution of capillary forces (either for the contacts or the dis- tant interacting pairs), σ N e αβ is the contribution of normal elastic forces, and σ T αβ is the contribution of tangential forces. Conversely, the σ d 22 /σ 22 is steadily decreased until the end of unloading path. When the vertical stress decreases to very low values, σ * 1 ≤ 10 -1 , the σ d 22 /σ 22 is remarkably dierent in four cases. Because in the looser system, the total stress σ 22 is smaller than in the denser system (as we already mentioned in the discussion on K 0 ). The σ d 22 /σ 22 of the F30OV1Vm3 (red curve) is therefore smaller than for the other cases. for F30OV1Vm3 (red curve), while the other case have higher ratios (dark-green, blue and magenta curves). Similar observations were discussed above for the smallest case of density. Additional information on the plateau of Regime I in both oedometric and isotropic models is provided in Fig. 4.12. In general, all P.D.Fs are roughly symmetric about 0 and diverge near 0 as a power law (with an exponent between 0.2 and 1.0) in the range of forces between 2 × 10 -4 F 0 and 10 -3 F 0 ). It is clearly apparent that at σ * 1 = 10 -2 , the P.D.F of the isotropic case (dark-green curve) is wider than the one of the oedometric case (magenta curve). The P.D.F of the isotropic case vanishes at -F 0 and F 0 while confrontation it vanishes at -0.5F 0 and 0.5F 0 for the oedometric test. That 8×10 0 kPa while its intervals in experimental curves are 7×10 -2 kPa ≤ σ 1 ≤ 2×10 1 kPa (OD01), and 7×10 -2 kPa ≤ σ 1 ≤ 3×10 1 kPa (CT01). However, for CT01, the void ratio at the end of the collapse is higher than the ones of OD01 as well as the F30OV1Vm3.

Oedometric versus isotropic compression

This reects the eect of the cell's size and the friction between grains and cell. This friction contributes to prevent the settlement of all grains. Therefore, the collapse of experimental specimen is slower than the simulated one. Under higher stress, at the end of loading path, the minimal void ratio of OD01 is lower than the simulated one. With the cell φ20, the void ratio is still higher than the ones of the cell φ50 and simulation.

This can be explained by the eect of the cell size, already mentioned in Chapter 2. In addition to possible friction on the cide walls, experimental test results are likely aected by initial large scale sample heterogenities evidence in the XRCT-scans of Chapter 2. As to the dierences in compression index and in the reduced stress corresponding to collapse, the results of the numerical parametric study indicate that it could be reduced with a choice of small RPR parameters and an appropriate initial connectivity, determined by agitation intensity V 0 . However, a more quantitative assessment of size and boundary eects in experiments is needed.

Geometrical structure change

In this section, only the case F30OV1Vm3 (see more details in Table 4.1) of the oedometric compression test is chosen in order to compare with the experimental results (from the local scan observations mentioned in Chapter 2). almost constant and insignicantly decreases at the starting of the collapse. Meanwhile, in Regime I of the experimental case, the initial structure is slightly collapsed with a very small rearrangement. Moreover, the rst local scan was performed at h 1 = 9.8 mm, not at the initial height (see Table 2.2 in Chapter 2). z exp is therefore higher than the z(h) at this stage. For the other scans (S06G, S09G, and S12G), z exp is higher than z(h = 0) (red curve) while it is globally lower than the other values of z(h). In experiments, a tolerance of 1 voxel is mainly used in all calculations. This value is equal to h/ d = 0.07 which is corresponding to z(h 4 ) at h/d = 0.075 (cyan curve). However, z exp at the S09G and S12G is signicantly lower than z(h 4 ). This is due to the heterogeneity of specimen. There still exists several large pores during the compaction test even at the highest level of compaction (S12G). This phenomenon is already described in Chapter 2 (see Fig. 2.26). Conversely, in simulation, the specimen is denitely homogeneous during the compression cycle, especially at the high pressure.

Nevertheless, when counting the lost contacts, z lc exp approaches closely to the z(h) 4 at h/d = 0.075, especially for the scan S06G. In other words, the total coordination number obtained in experiment is fairly underestimated when compared with the coordination numbers of neighbors in simulation.

Conclusion

In this Chapter, we have performed the simulation of oedometric compression tests for the monodisperse system without RPR. The results of this system are then compared in turn with the ones of the isotropic compression tests (in Chapter 3) and with the ones of the local scan observations (in Chapter 2). The results obtained from the simulation of oedometric compression tests, including the compaction behavior and the microstructural characterizations, are globally similar with the results of isotropic tests, especially if mean stress is regarded as the control parameter for both tests. The stress ratio, K 0 , in oedometric compression, as studied here for initially isotropic states, depends on initial connectivity and goes through a minimum in the course of oedometric compression, in Regime II (from 0.3 in initially well connected systems to about 0.5 in barely rigid ones. As the minimum void ratio is approached, it increases to a common value near 0.8. Strongly overconsolidated states with small σ * m may exhibit comparatively very large stresses in some directions, of order zΦF 0 /d 2 .

A comparison of compression behavior between the oedometric compression test and the experiments reveals several remarks. Regime I exists in the experimental tests with a very small rearrangement of the initial structure but the plateau of Regime I does not exist. In other words, it is really dicult to maintain the stable-state of initial structure under very low stress in experiment without rearrangement. Furthermore, the dierence of the collapse between simulation and experiment obviously conrms the signicant eect of the friction between the cell and grains during the test.

For the prediction of z, although the experimental result is underestimated in comparison to the DEM simulation result, we obtained the correct trends for the changes of the microstructure in both oedometric simulation and experimental tests.

General conclusions

The present work is a contribution to the understanding of the mechanical behavior of wet granular materials in very loose states. We have proposed both a numerical tool based on the DEM and experimental procedures to study the 3D plastic response of glass beads subjected to isotropic and oedometric compression tests. Microstructural changes of the specimen induced by increasing applied external force were also characterized in both experiments and DEM simulations. The main conclusions of this work are briey summarized as follows.

Firstly, we can produce in the laboratory specimens of wet granular soils presenting a very loose yet stable-state, even using a very small water content. This is one of the main originalities of the present work as compared with past studies which focused on denser states and higher initial water contents (e.g. Bruchon et al. [START_REF] Bruchon | Full 3d investigation and characterisation of capillary collapse of a loose unsaturated sand using x-ray ct[END_REF]). Similar solid fractions are obtained by a suitable numerical aggregation process prior to compression.

Secondly, the compression behavior of wet beads at loose states obtained from experiment and DEM simulations (isotropic and oedometric compression tests) has usually three regimes: (i) at stresses lower than an apparent precompression stress, the initial structure can sustain the load without any rearrangement (in the DEM simulations)

or with very small rearrangements (in the experiments); (ii) when the stresses exceed this apparent precompression stress, the loose structures collapse and restructure; and

(iii) at higher stresses, only elastic deformation of grains at contact points occurs. This behavior abides by the classical logarithmic law of soil mechanics in the intermediate stage. The compression index is not an intrinsic material property, but depends on its initial density and connectivity. This latter feature of initial loose states depends on the aggregation process, and aecting the subsequent behavior in the earlier stages of plastic compression. The properties of wet granular soils in plastic compression in the later stages of plastic collapse are not too much dependent on the assembling process.

Thirdly, we presented a numerical approach to precisely measure several fundamental microstructure properties of loose systems from 3D tomography images. Although this approach underestimates several microstructure properties in comparison with DEM 174

Conclusions and Perspectives simulations and other studies in the literature, it provides an appropriate characterization of the grain-scale behavior of wet granular soils at very loose state during compression.

Fourthly, we emphasized the important dierences between cohesive and non-cohesive granular assembles. The most important dierences, which are brought by cohesion in macroscopic behavior, are the existence of stable very loose structures and the response during compression. Moreover, the inuence of size polydispersity was investigated only under isotropic compression. Our results show that the size polydispersity negligibly aects the plastic behavior of wet beads.

Finally, our results also emphasize the remarkable inuence of the rolling and pivoting resistances at contacts between grains. Even a small RPR is enough to produce strong dierences in the plastic behavior, with the collapse under growing conning stress becoming notably more gradual.

Perspectives

The research program developed here can be extended in many potential research directions, as follows.

On the experimental side, in order to increase the accuracy of the microstructural characteristics in experiments, the entire specimen should be used for the detection process. However, this process requires a very large memory for the calculation. So the detection for the REV (mentioned in Sec. 2.4.4 of Chapter 2) is preferably proposed for the detection process. Moreover, the sub-voxel algorithm could be also applied in order to increase the resolution of the tomography images. Furthermore, the algorithm of detection should be improved to obtain the real value of the detected radii instead of the rounded value.

The pore-scale deformation of specimens during compression can be analyzed in details by using the Digital Volume Correlation (DVC or 3D DIC) method. This is one of the most interesting phenomenon of wet granular materials at very loose state. Besides, the pore size distribution can be directly analyzed from the 3D tomography images based on the preliminary analyses in Sec. 2.4.4 of Chapter 2.

On the numerical side, the quasistatic behavior of a 3D model, in pendular state, accounting for the possible eect of rolling and pivoting resistances at contacts under triaxial compression tests where the macroscopic deviatoric stress will play a major role, is a straightforward extension of the present work. Ref. [START_REF] Than | Basic mechanical properties of wet granular materials: A DEM study[END_REF] contains a preliminary comparison with the properties of the critical state in quasistatic simple shear ow. The wide range of stress and microstructural states between isotropically or oedometrically compressed states and steadily deformed materials at large strains remains to be explored.

In addition, the role played by solid bridges (e.g. cemented granular soils or methane hydrate soils) instead of capillary bonds should also be investigated.

  An example of XRCT scans of Hostun sand, with D is the diameter of cell and r is the resolution of image . . . . . . . . . . . . . . . . . . . . . 1.10 Examples of partitioning and subsequent ball-and-stick network construction, applied to an unconsolidated sand imaged at 2.8 µm . . . . . . . . . 1.11 Four dierent domains (fully saturated, funicular, pendular and hygroscopic) shown in slices from trinarised 3D images of Hostun sand . . . . . 1.12 Eect of soil compaction on interpores (a) and intrapores (b) volumes from soil aggregates with varying levels of soil moisture and compaction . 1.13 The cycle of the DEM calculation . . . . . . . . . . . . . . . . . . . . . . 1.14 The packing fraction of dierent sized glass beads as a function of the eective gravitational acceleration . . . . . . . . . . . . . . . . . . . . . . 1.15 Ballistic deposit before and after compaction by a constant external force on the piston . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.16 Compression and decompression curves: eect of initial agitation level in aggregation stage and inuence of rolling resistance parameter. . . . . . . 1.17 A map of tensile (green) and compressive (red) forces in a thin layer cut out in the packing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.18 Schematic view of liquid bridge between two smooth (a) polydisperse spheres and (b) monodisperse spheres. . . . . . . . . . . . . . . . . . . . xiv List of Figures 1.19 Capillary force and capillary pressure as a function of the gorge's radius for dierent values of contact angle θ. . . . . . . . . . . . . . . . . . . . . 1.20 Evolution of the capillary force F cap as a function of the interparticle distance h for a given capillary pressure ∆u. . . . . . . . . . . . . . . . . 1.21 Loose cohesive assembly of disk . . . . . . . . . . . . . . . . . . . . . . . 1.22 An example of compression curve of loose systems obtained with 2D and 3D models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.23 Illustration of the clumping of spheres method with no overlapping particles. 1.24 Pore size distributions of granules produced under high shear . . . . . . . 1.25 Comparison of coordination measurements on numerical sphere packings produced by DEM and by analysis of XRCT images. . . . . . . . . . . . 1.26 Comparison of the RDFs of (a) the monodisperse glass packing structure, and (b) polydisperse glass packing structure . . . . . . . . . . . . . . . . 1.27 (a) Cemented numerical sample and (b) laboratory sample. . . . . . . . .

2 . 7

 27 Applied force versus time during compression process. The inset shows the signal noise during the initial stage when the specimen approached and touched the piston . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Global view of in-situ compression test combined with XRCT at Laboratoire Navier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 (a) The PMMA cell with several gaskets and pistons (heights δ given in mm, see Fig. 2.10 for denition of δ). (b) Several wedges with dierent values of thickness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Illustration of the compressed heights of in-situ compaction test that were scanned for XRCT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 Void ratio versus vertical stress during oedometric compression tests for cells φ50 and φ20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 3D reconstructed images of the specimen at dierent stages of compaction 2.13 Illustration of grey level histograms for all global scans . . . . . . . . . . 2.14 Images at dierent steps of segmentation: (a) original image, (b) segmented image . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.15 Illustration of grey level as a function of the height of specimen . . . . . 2.16 Schematic view of several types of diagram for calculating the heterogeneity List of Figures xv 2.17 Void ratio distribution according to the compressed height of sample for dierent compaction steps (diagram 1) . . . . . . . . . . . . . . . . . . . 2.18 Void ratio distribution according to the average radii of the selected cylinders for dierent compaction steps (diagram 2) . . . . . . . . . . . . . . . 2.19 Void ratio distribution according to eight equal sectors for dierent compaction steps (diagram 3) . . . . . . . . . . . . . . . . . . . . . . . . . . 2.20 Void ratio distribution according to ve equal volumes for dierent compaction steps (diagram 4) . . . . . . . . . . . . . . . . . . . . . . . . . . 2.21 Void ratio distribution according to the x (a) and y (b) directions for dierent compaction steps (diagram 5) . . . . . . . . . . . . . . . . . . . 2.22 Total volume and total pore volume of volume V2 for dierent steps of compaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.23 An example of the local thickness for volume V2 (a) at step S01G and a slice (b) extracted from (a) . . . . . . . . . . . . . . . . . . . . . . . . . . 2.24 Distribution of pore volume of the volume V2 for dierent levels of compaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .2.25 Position of the local scan (b) from the global scan (a) and an example of investigated cube (c) . . . . . . . . . . . . . . . . . . . . . . . .

3 . 9 (

 39 (a) e versus log P * for reference system in compression cycle; (b) comparison of reference case and typical non-cohesive behavior . . . . . . . . 3.10 Eect of dierent (isotropic) unloading and reloading histories on void ratio120 3.11 a) Coordination numbers z, z c and z d ; and (b) coordination numbers z c , z + , and z -versus P * in compression and decompression cycles . . . . . . 3.12 Contact networks in loading (a, b, and c) and unloading (d, e, and f ) cycle 3.13 Networks of compressive contacts (red lines), tensile ones (green lines), and contacts with F N = 0 (magenta) . . . . . . . . . . . . . . . . . . . . 3.14 Distribution of normal forces for dierent values of P * from the compression path, (a) normalized by the maximum tensile force F 0 , and (b) normalized by the average normal force F N . . . . . . . . . . . . . . . . 3.15 Distribution of normal forces for dierent values of P * in non-cohesive case, normalized by the average normal force F N . . . . . . . . . . . . . 3.16 Void ratio (a), coordination numbers z c and z d (b) versus P * of the reference case for dierent values of meniscus V m /d 3 . . . . . . . . . . . . . . 3.17 Coordination numbers of z + and z -versus reduced pressure P * for dierent values of V m /d 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Distribution of the number of contacts per grain of the reference case for dierent values of V m /d 3 at three typical states of P * . . . . . . . . . . . 3.19 Evolution of coordination number of close neighbors z(h) versus dimensionless interparticle distance h/d for dierent values of P * in loading path 3.20 Void ratio versus pressure P , as cohesive forces are suppressed at the beginning of unloading . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 (a) Void ratio e, and (b) coordination numbers z c and z d versus P * for dierent values of initial solid fraction Φ 0 . . . . . . . . . . . . . . . . . . 3.22 Coordination numbers z + and z -versus P * for dierent values of Φ 0 . . List of Figures xvii 3.23 Distribution of the number of contacts per grain for dierent values of Φ 0 at three typical states of P * . . . . . . . . . . . . . . . . . . . . . . . . . 3.24 Distribution of normal forces for dierent values of Φ 0 under low pressure P * = 10 -2 , is normalized by the maximum tensile force F 0 . . . . . . . . 3.25 (a) Compression and decompression curves, and (b) coordination numbers z c and z d for dierent values of V 0 /V * . . . . . . . . . . . . . . . . . . . . 3.26 Coordination numbers z + and z -versus P * for dierent values of V 0 /V * 3.27 A 3D slice of specimen for dierent values of V 0 /V * at the beginning of loading path (P * = 10 -3 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.28 Number of contacts per grain for dierent values of V 0 /V * at three typical states of P * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.29 Distribution of normal forces for dierent values of V 0 /V * under P * = 10 -3 , normalized by F 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.30 (a) Void ratio e and (b) coordination numbers z c and z d versus P * for dierent values of V 0 /V * . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.31 Contacts (red lines) and distant interactions (blue lines) for dierent values of V 0 /V * at P * = 10 -3 . . . . . . . . . . . . . . . . . . . . . . . . . . 3.32 Grains contained in a slice (thickness 3d min ), at low P * , for dierent values of V 0 /V * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.33 (a) Void ratio e and (b) coordination numbers z c and z d versus P * for MDS and PDS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.34 (a) Void ratio e and (b) coordination numbers z c and z d versus P * in compression cycle for dierent values of Φ 0 . . . . . . . . . . . . . . . . . 3.35 Coordination numbers z + and z -versus reduced pressure P * for dierent values of Φ 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.36 (a) Void ratio e and (b) coordination numbers z c and z d versus P * in compression cycle for dierent values of V 0 /V * . . . . . . . . . . . . . . . 3.37 Coordination numbers z + and z -versus P * for dierent values of V 0 /V * . 3.38 3d-slices of specimens assembled with RPR for dierent values of V 0 /V * , under P * = 10 -3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.39 Number of contacts per grain for six values of V 0 /V * at three typical states P * = 10 -3 , P * = 10 -1 , and P * = 10 1 . . . . . . . . . . . . . . . . . . . . 3.40 Distribution of normal forces at low pressure P * = 10 -3 for dierent V 0 /V * . P.D.F normalized by F 0 . . . . . . . . . . . . . . . . . . . . . . . 3.41 (a) Void ratio e and (b) coordination numbers z c and z d versus P * for dierent values of µ R = µ P . . . . . . . . . . . . . . . . . . . . . . . . . . 3.42 Distribution of normal forces for dierent values of µ R /d at P * = 10 -2 . . 3.43 (a) Void ratio e and (b) coordination numbers z c and z d versus P * in compression cycle for systems both without RPR and with RPR . . .

( a ) 6 ( 8 K

 a68 Void ratio e, and (b) coordination numbers z c and z d versus σ * 1 for dierent values of V 0 /V * . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Coordination numbers z + (compressive contacts) and z -(tensile contacts) versus reduced axial stress σ * 1 for dierent values of V 0 /V * . . . . . . . . 4.3 Coecient of lateral pressure K 0 versus reduced vertical stress σ * 1 for dierent values of V 0 /V * . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Contribution of distant interactions to normal stress σ 22 . Red arrow denotes the increasing of values . . . . . . . . . . . . . . . . . . . . . . . . 4.5 a) Contributions of tangential (a) and capillary (b) interactions to normal stress σ 22 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.a) Void ratio e, and (b) coordination numbers z c and z d versus σ * 1 for dierent values of initial solid fraction Φ 0 . . . . . . . . . . . . . . . . . . 4.7 Coordination numbers of z + and z -versus reduced vertical stress σ * 1 for dierent values of Φ 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.0 versus reduced axial stress σ * 1 for dierent values of Φ 0 . . . . . . . . 4.9 Contribution of distant interactions to normal stress σ 22 for dierent values of Φ 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 (a) Contributions of tangential (a) and capillary (b) interactions to normal stress σ 22 for dierent values of Φ 0 . . . . . . . . . . . . . . . . . . . . . . 4.11 (a) Void ratio e, and (b) coordination numbers z c and z d versus reduced mean stress σ * m for the minimum and maximum V 0 /V * in both simulations of oedometric and isotropic tests . . . . . . . . . . . . . . . . . . . . . . . 4.12 Distribution of normal forces for both oedometric and isotropic models . 4.13 Oedometric compression curves in experiments and simulations (no RPR) for dierent values of Φ 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Coordination numbers of neighbors z(h) versus vertical stress σ 1 and the total coordination number z versus positions of scans . . . . . . . . . . .List of Tables

1. 1 . 9 Figure 1 . 1 :

 1911 Figure 1.1: Particles displacement (color shows the displacement magnitude) during triaxial compression. Results are shown as vertical slices. Subgures from left to right correspond to increasing vertical stresses (from Andò et al. [56]).

Figure 1 . 2 :

 12 Figure 1.2: Photoelastic image of a small assembly of disks (from Dantu [57]).

Figure 1 . 3 :

 13 Figure 1.3: Skeletal versus contact-level capillary and electrical forces. The upper part

Figure 1 . 4 :

 14 Figure 1.4: Schematic view of a laboratory XRCT scanner

Figure 1 . 5 :

 15 Figure 1.5: One-dimensional compression plots for carbonate and silica sands (from Golightly [136]).

Figure 1 . 6 :Figure 1 . 7 :

 1617 Figure 1.6: Capillary bridges in glass beads (from Fournier et al. [12])

Figure 1 . 8 :

 18 Figure 1.8: The dierences of water content for dierent regions 1, 2, 3, and 4 as a function of time (from Crestana et al. [103]).

Figure 1 . 9 :

 19 Figure 1.9: An example of XRCT scans of Hostun sand, with D is the diameter of cell and r is the resolution of image: (a) D=70 mm, r=50 µm, (b)D=13 mm, r=10 µm (from Bruchon [148])

Figure 1 . 10 :

 110 Figure 1.10: Examples of partitioning and subsequent ball-and-stick network construction, applied to an unconsolidated sand imaged at 2.8 µm. (a) Zoom into slice from original image, (b) segmented image, (c) partition of the pore space using watershed transform, (d) partition of the grain space, (e) equivalent pore-network of 3D region surrounding the 2D slice, and (f ) equivalent grain network (from Wildenschild and Sheppard [147]).

  1.11), complete continuity of the liquid phase is lost, with liquid bridges between the grains. At bottom right of Fig.1.11, the hygroscopic domain with only adsorbed water on the surface of the grains. The hygroscopic domain was also observed in some regions of the sample[START_REF] Riedel | Water retention behaviour explored by x-ray ct analysis[END_REF].

Figure 1 . 11 :

 111 Figure 1.11: Four dierent domains (fully saturated, funicular, pendular and hygroscopic) shown in slices from trinarised 3D images of Hostun sand; water in light blue, air in black and solid in light brown. The degrees of saturation noted are global (from Riedel et al. [146]).

Figure 1 . 12 :

 112 Figure 1.12: Eect of soil compaction on interpores (a) and intrapores (b) volumes from soil aggregates with varying levels of soil moisture and compaction. Note: W1 refers to moisture content of 9.3 %, W2 represents 18.3 %; BD1 and BD2 refer to bulk density increment of 0.28 and 0.71 g.cm -3 , respectively (from Menon et al. [149]).

Figure 1 . 13 :

 113 Figure 1.13: The cycle of the DEM calculation

Figure 1 . 14 :

 114 Figure 1.14: The packing fraction of dierent sized glass beads as a function of the eective gravitational acceleration. Points are the measured results: triangle, d = 500 µm;

Fig. 1 .

 1 Figure 1.15: Ballistic deposit before and after compaction by a constant external force on the piston. Contact properties include Coulomb friction, rolling friction and cohesion

Figure 1 . 34 Chapter 1 . Literature review 1 . 3 . 3

 1341133 Figure1.17: A map of tensile (green) and compressive (red) forces in a thin layer cut out in the packing. Line thickness is proportional to the magnitude of the force (from Richefeu et al.[START_REF] Richefeu | Stress transmission in wet granular materials[END_REF]).

Figure 1 . 18 :

 118 Figure 1.18: Schematic view of liquid bridge between two smooth (a) polydisperse spheres and (b) monodisperse spheres (not to scale).

Figure 1 . 19 :

 119 Figure 1.19: Capillary force and capillary pressure as a function of the gorge's radius for dierent values of contact angle θ (from Bruchon [148]).

Figure 1 . 20 :

 120 Figure 1.20: Evolution of the capillary force F cap as a function of the interparticle distance h for a given capillary pressure ∆u: a meniscus can form when h ≤ D c and breaks down when h ≥ D r (from Chareyre and Scholtès in [172]).

Figure 1 . 21 :

 121 Figure1.21: Loose cohesive assembly of disks (from Gilabert et al.[START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF]). Force magnitude is encoded as line thickness, with compressive forces in red and tensile ones in green.Blue lines correspond to null normal contact forces.

Figure 1 . 22 :

 122 Figure 1.22: An example of compression curve of loose systems obtained with 2D [14]and 3D[START_REF] Than | Basic mechanical properties of wet granular materials: A DEM study[END_REF] models.

Figure 1 .

 1 Figure 1.23: (a) Illustration of the clumping of spheres method with no overlapping particles. From left to right: the tomographic image, the actual clumping method, and the simplication of particle shape used by DEM. (b) Illustration of the clumping of shapes method with overlapping particles (from Moreno-Atanasio et al. [50]).

Figure 1 . 24 :

 124 Figure 1.24: Pore size distributions of granules produced under high shear (a) calculated from tomographic images; (b) measured by mercury porosimetry (from Farber et al. [144]).

  [START_REF] Marmottant | Coordination measurements in compacted NaCl irregular powders using X-ray microtomography[END_REF]. The authors observed the structure of compacted NaCl irregular powders, the coordination number is shown in Fig.1.25.

Figure 1 . 25 :

 125 Figure 1.25: Comparison of coordination measurements on numerical sphere packings produced by DEM and by analysis of XRCT images of (a) the mean coordination number z, and (b) the mean contact area A between touching particles normalized by the square of the particle radius R (from Marmottant et al. [210]).

Figure 1 . 26 :

 126 Figure 1.26: Comparison of the RDFs of (a) the monodisperse glass packing structure, and (b) polydisperse glass packing structure obtained from XMT (solid curve) and DEM simulation (dashed curve).

  particle displacements have not been used in conjunction with DEM, it oers a direct way of linking simulations and experiments. Here, Fig. 1.27 shows an example of the comparison between the DEM and the experiment under triaxial tests of cemented sand.

Fig. 1 .Figure 1 .

 11 Fig.1.27 shows the nal numerical specimen (a)[START_REF] De Bono | Investigating the micro mechanics of cemented sand using dem[END_REF], alongside the equivalent laboratory sample (b)[START_REF] Marri | Drained behaviour of cemented sand in high pressure triaxial compression tests[END_REF]. A conjugate shear plane can be observed clearly in the laboratory sample. A similar failure can be seen in the corresponding image of the numerical sample, which has fairly distinctive non-uniform deformations.

  beads are large enough to be satisfactorily resolved by the XRCT. The specic density of this material is equal to 2,460 kg.m -3 . The grain size distribution of the chosen glass beads is shown in Fig.2

Figure 2 . 1 :

 21 Figure 2.1: Cumulative grain size distribution of glass beads.

Figure 2 . 2 :

 22 Figure 2.2: Specimen preparation for cell φ50. (a) Metallic cell φ50, (b) specimen preparation, (c) specimen before testing, and (d) specimen after testing.

Figure 2 . 3 :

 23 Figure 2.3: Specimen preparation for cell φ20. (a) PMMA cell φ20, (b) specimen before testing, and (c) and (d) specimens after testing.

  2.2(b)). This action must be done carefully and gently in order to avoid any perturbation to the initial structure. Height of specimen is 20 mm in the case of cell φ50 and 10 mm in the case of cell φ20 (see Figs. 2.2(a) and 2.3(a)). The specimens before testing are shown in Figs. 2.2(c) and 2.3(b).

Figure 2 . 4 :

 24 Figure 2.4: Schematic view of the system for preparing the specimen.

Figure 2 . 5 :

 25 Figure 2.5: Initial void ratio e 0 and initial porosity n versus height of free fall h f .

Figure 2 . 6 :

 26 Figure 2.6: Schematic view of oedometric compression system.

  2.2(d)and 2.3(d)) to measure their mass and nal height. These measures were then used to validate the quantities measured after preparation and shown in Table2.1. These nal measurements allows verifying the nal state of the specimens after compression.

Figure 2 . 7 :Figure 2 . 8 :

 2728 Figure 2.7: Applied force versus time during compression process. The inset shows the signal noise during the initial stage when the specimen approached and touched the piston.

Chapter 2 .

 2 Experimental investigations also made with dierent values of thickness, as shown in Fig.2.9(b). These wedges are used to alter accurately the compressed height during the in-situ compression tests.

Figure 2 . 9 :

 29 Figure 2.9: (a) The PMMA cell with several gaskets and pistons (heights δ given in mm, see Fig. 2.10 for denition of δ). (b) Several wedges with dierent values of thickness.

Figure 2 . 10 :

 210 Figure 2.10: Illustration of the compressed heights of in-situ compaction test that were scanned for XRCT. The dashed-region corresponds to a 10-mm diameter cylinder and denotes the central part of specimen imaged during local scans.

Fig. 2 .

 2 [START_REF] Gröger | Modelling and measuring of cohesion in wet granular materials[END_REF] shows the compression curves for dierent values of initial void ratios. All the tests were duplicated to check repetability. The initial stress at the bottom of specimens including their weight varies equals 0.048 kPa and 0.024 kPa for cells φ50 and φ20, respectively.

Figure 2 . 11 :

 211 Figure 2.11: Void ratio versus vertical stress during oedometric compression tests for cells φ50 and φ20.

Figure 2 . 12 :

 212 Figure 2.12: 3D reconstructed images of the specimen at dierent stages of compaction.

  From (a) to (h), the images correspond to the global scans in

Figure 2 . 13 :

 213 Figure 2.13: Illustration of grey level histograms for all global scans.

  Fig. 2.14 shows an example (horizontal slice) of the segmentation process corresponding to two 66 Chapter 2. Experimental investigations consecutive stages: a) initial state (original image) and b) segmented image (or binarized image). After thresholding, the aggregates (black) and the air (white) are clearly distinguished. The void ratio is then recalculated by dividing the number of black voxels by the total number of voxels of the 3D segmented image.

Figure 2 .

 2 Figure 2.14: Images at dierent steps of segmentation: (a) original image, (b) segmented image (or binarized image).

Figure 2 . 15 :

 215 Figure 2.15: Illustration of grey level as a function of the height of specimen.

  order to observe the spatial distribution of void ratio of all global scans during the compaction test, ve types of diagram are proposed in Fig. 2.16. Firstly, the sample is divided in ten equal layers h i (Fig. 2.16(a)). Ten concentric hollow cylinders with average diameter R i are then distinguished (Fig. 2.16(b)). Thirdly, the sample is split into eight equal sectors (Fig. 2.16(c)). Fourthly, ve dierent volumes, with the size L×B×H = 10×10×H mm 3 = 770×770×H voxels 3 (H is specimen's height) are investigated (Fig. 2.16(d)). Twelve dierent volumes according to X and Y directions are nally identied (Fig. 2.16(e)). Their size is L×B×H 3.25×3.25×H mm 3 = 250×250×H voxels 3 . Hereafter, the ve above diagrams are referred to as diagram 1, diagram 2, diagram 3, diagram 4, and diagram 5, respectively.

Figure 2 . 16 :

 216 Figure 2.16: Schematic view of several types of diagram for calculating the heterogeneity. (a) Ten equal layers, (b) ten concentric cylinders, (c) eight equal sectors, (d) ve dierent volumes, and (e) dierent volumes according to X and Y directions.

Figure 2 . 17 :Figure 2 . 18 :

 217218 Figure 2.17: Void ratio distribution according to the compressed height of sample for dierent compaction steps (diagram 1).

For the diagram 4 (Figure 2 . 19 :Figure 2 . 20 : 4 . Global scan observations 71 Fig. 2 .Figure 2 . 21 :Figure 2 . 22 :Figure 2 . 23 :

 42192204712221222223 Figure 2.19: Void ratio distribution according to eight equal sectors for dierent compaction steps (diagram 3).

Fig. 2 .Figure 2 . 24 :

 2224 Fig. 2.23 shows an example of volume V2 at step S01G after applying the local thickness algorithm. In Fig. 2.23, the grains (only the black color) and the pores (the other colors) are clearly distinguished. Bright colors represent the larger pores. This can be easily observed via the bright volume in Fig. 2.23(a). Based on the histogram of the 3D images, the volumes of dierent pores are quickly calculated. Nevertheless, we can only compute the volume of the pores having the same size (or same color). Therefore, we cannot distinguish exactly the volume of every pore in the volume V2.

2. 5

 5 Local scan observations 2.5.1 General principles As described in Sec. 2.2.4.3, four local scans, with 10 mm in diameter, were captured from the central region of the global scans. An example of 3D reconstructed image from local scan is illustrated in Fig. 2.25. From the 3D tomography images of the local scans, cubical subregions (edge length = 1.95 mm) were extracted. One such subdomain is shown in Fig. 2.25(c). Ten such cubes (hereafter so-called Standard Volume, or SV) were randomly extracted at dierent positions within the local scan, to take into account the sample microstructure heterogeneity.

Figure 2 . 25 :Figure 2 . 26 :

 225226 Figure 2.25: Position of the local scan (b) from the global scan (a) and an example of investigated cube (c).

Figure 2 . 27 :

 227 Figure 2.27: The algorithm to detect the spherical structures.

  weight value is added to the voxels in the accumulation array upon the line segment which is dened by the vector ∇I(i, j, k). The length of the line segment is set to be the possible maximum diameter of the spheres, and the magnitude of ∇I(i, j, k) is used as the weight value. The accumulation array is built from the collecting of all nonzero gradient vectors from many directions. In Fig.2.28, the darkest voxel has the maximum weight.

Figure 2 . 28 : 78 Chapter 2 .

 228782 Figure 2.28: Construction of the accumulation array from the gradient eld (the added nonzero gradient vector to the accumulation array for the left, and the accumulation array after adding the nonzero gradient vector in dierent directions for the right) (not to scale).

Figure 2 . 29 :

 229 Figure 2.29: Denition of the signature curve (not to scale).

Fig. 2 .

 2 Fig.2.29 depicts the denition of the signature curve, in which (i c , j c , k c ) is the center position of a detected sphere. For a voxel (i, j, k) in the neighborhood of (i c , j c , k c ), let V (i, j, k) be the gradient vector of the voxel, (i, j, k) and q(i, j, k) be the vector joining

Fig. 2 .

 2 Fig. 2.30 shows an example of the signature curve. From this curve, the radius of the sphere corresponds to the distance to center position having the highest value of f (r). In the example shown, the radius equals seven times the voxel size. It is equivalent to 45.5 microns. Besides, other values of radius are also detected: 6 and 8 times the voxel size (equivalently 39 microns and 52 microns). Note that, in the present work, only rounded valued of radius were considered.

Fig. 2 .

 2 Fig. 2.31 illustrates an example for one slice in a 3D image. The blue disks are the detected spheres, and the red points represent their centers. Incomplete white disks at

Figure 2 . 30 :

 230 Figure 2.30: An example of a signature curve computed from the 3D image.

Figure 2 .

 2 Figure 2.31: A slice in a 3D detected image (not to scale).

  2.5.1) is created. It is called the Extended Volume (or EV) with the size b (see Fig. 2.32) which b = a+2d max , where a is the size of the SV, d max is the maximum diameter of the detected spheres. The EVs were then detected with the same input parameters and compared with the SVs. All grains inside the SVs and the

Figure 2 . 32 :

 232 Figure 2.32: Method to nd the lost contacts, inside cube is the standard volume; outside cube is the extended volume (not to scale).

Fig. 2 .

 2 Fig. 2.33 shows an example of a 3D reconstructed cube before and after detection. Each cube, from top to bottom, corresponds to each local scan. On the left-hand side, the cubes are directly reconstructed from 3D images (before detection); on the right-hand side, the cubes are reconstructed from the detected spherical structures. As mentioned in Sec. 2.5.2.3, the incomplete spheres, located on the border of the SVs, are not identied.

Figure 2 . 33 :Fig. 2 .Figure 2 . 34 :Fig. 2 .

 23322342 Figure 2.33: Example of 3D reconstructed cubes of four local scans. Left-hand side: before detection, right-hand side: after detection. Edge length = 1.95 mm.

Fig. 2 . 2 . 5 .

 225 Fig.2.35(a) shows the assembled spheres within the specimen, including the stacks (lighter grey levels) of spheres and the pores (darker grey levels). In this specimen, several large pores exist corresponding to a size up to 1 mm and a lot of stacks created from the spheres have dimension ranging from 0.3 to 0.5 mm. X-ray diraction (the opaque layers around the spheres) can be precisely shown in Fig.2.35(b), by zooming in several spheres in the specimen (a) and reducing the noise of original images. All small cubes (Figs. 2.35(b, d, f, and h)) show clearly the contact region at the boundary of particles.

Figure 2 . 35 :

 235 Figure 2.35: Example of 3D reconstructed specimen: (a) whole specimen in original images, (b) small cube extracted from (a), (c) solid part of specimen (binarized images), (d) small cube extracted from (c), (e) and (g) whole reconstructed specimens with tolerance 0 and 1 voxel, and (f ) and (h) small reconstructed cubes extracted from (e) and (g).

2. 5 .

 5 Local scan observations 85 2.5.6 HeterogeneityAs mentioned in Sec. 2.5.1, ten SVs were randomly extracted at ten dierent positions for each local scan. After the detection process, the number of particles for the ten SVs considering four steps (S03L, S06L, S09L, and S12L) are obtained, as shown in Fig.2.36.

Figure 2 . 36 :

 236 Figure 2.36: Number of particles of the ten SVs for four local scans after detection process.

Figure 2 . 37 :

 237 Figure 2.37: Macroscopic solid fraction (triangle) and the average solid fraction of 4 local

Fig. 2 .Figure 2 . 38 :

 2238 Fig.2.38 plots the average total coordination number of four steps for three values of tolerance (0, 1, and 2 voxels). For a tolerance of 0, z tends to decrease under growing

Fig. 2 .

 2 39 illustrates the average coordination number of close neighbors of four scans with the dimensionless interparticle distance h/ d , in which d is dened as follows: d = 2d max d min d max + d min .

Figure 2 . 39 :

 239 Figure 2.39: Coordination numbers of close neighbors for four local scans.

Figure 2 . 40 :

 240 Figure 2.40: Number of contacts per grain for three tolerances of 0, 1 and 2 voxels.

Figure 2 . 41 :

 241 Figure 2.41: An example of a sphere with one contact.

Figure 2 . 42 :

 242 Figure 2.42: Radial distribution functions for particle centers of four local scans.

Fig. 2 .

 2 Fig. 2.42 shows the radial distribution of particle centers in a spherical shell for four steps of the dimensionless distance r/ d . Each curve is represented by the average value of ten SVs. The probability distribution function is nearly similar at various states with their oscillations and peaks due to large size distribution. At very short r (r/ d < 1)

Figure 2 . 43 :Figure 2 .

 2432 Figure 2.43: Contact network of a typical SV under growing of compaction. Red lines denote the contacted grains with tolerance of 1 voxel, blue lines denote the neighbors of grains with their interstice of 1 voxel.

(3. 2 )

 2 This diameter distribution is shown in Fig.3.1, in the form of the cumulated probability function,

Figure 3 . 1 :

 31 Figure 3.1: Grain size distribution of the polydisperse system.

3. 2 . 2

 22 Interaction law: normal forces 3.2.2.1 Elastic normal force: Hertz law Consider two spherical beads i and j, with respective diameters d i and d j , made of a solid elastic material with Young modulus E and Poisson ratio ν. Using notations 3.2. Model material and interaction laws 101

(3. 5 )

 5 Although this choice does not entail any loss of generality, numerical results in physical units correspond here to glass beads, with Young modulus E = 70 GPa and Poisson ratio ν = 0.3. Formulae (3.4) and (3.5) are used with monosized beads (with

  contact in the past, or broken meniscus)0 h > D r . (3.7)Here, h < 0 corresponds to an elastic deection of the particles in contact. F 0 = πdΓ cos θ is the maximum tensile force, involving surface tension Γ of the liquid interface with the atmosphere (or the ambient continuous uid phase), and wetting angle θ. Numerical results below are obtained with Γ = 7.27 × 10 -2 N.m -1 (an appropriate value for pure water at 20 o C under atmospheric pressure), and θ = 0, assuming perfect wetting.

Fig. 3 .

 3 Fig. 3.2 shows how the static normal forces, with the repulsive, elastic (Hertzian) and the attractive, capillary contributions, vary with distance h.

Figure 3 . 2 :

 32 Figure 3.2: Static normal force, F N = F e N + F cap , as a function of the interparticle distance h. Inset: Geometry of a liquid bridge between two particles, with lling angle ϕ and contact angle θ. Short-dashed line: normal elastic repulsive force, long-dashed line: resultant normal force.

  DEM simulations consist in computing the collective motion of collections of N spherical beads (with N = 4000 in most cases in this study) by integration of the dierential equations resulting from Newton's laws and rigid-body mechanics. Grains are regarded as homogeneous solids, and thus a bead of diameter d has a mass m proportional to d 3 and a moment of inertia I = md 2 /10. The grains are enclosed in a cuboidal cell. The cell edges are parallel to the coordinate axes (x α ) α=1,3 , with respective lengths (L α ) α=1,3 . Periodic boundary conditions (PBC) are used, thereby avoiding wall eects. Two neighboring spheres, i and j, if brought into

Figure 3 . 4 :

 34 Figure 3.4: Motion of two spheres.

(3. 23 ) 108 Chapter 3 .

 231083 In Eqs.(3.23), only those spheres j interacting with i will contribute to the sums on the right-hand side. On dealing with periodic, deformable cells, those equations should be adapted: the linear velocities comprise two terms, one due to the variation of reduced coordinates xα = x α /L α , the other one due to the variation of L α . One possibility, used Isotropic compression: DEM study in Ref.[START_REF] Agnolin | Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks[END_REF], is to write down an equation for the vector of reduced coordinates exactly in the form of (3.23), with scaled down forces in the right-hand side. With appropriate choices of parameters (in particular, of the mass associated with the cell dimensions) the classical form (3.23) of the equations is only very slightly perturbed.

110 Chapter 3 .

 1103 Isotropic compression: DEM study 3.3.4 Liquid content and limit of pendular state

  is the volume of voids. With N b liquid bridges of volume V m , between N particles of diameter d, the coordination number of liquid bridges is z = 2N b /N and the solid fraction is

Figure 3 . 5 :

 35 Figure 3.5: Three spheres connected with liquid bridges.

Fig. 3 .

 3 Fig. 3.5 depicts an example of three spheres connected with the closest possible liquid bridges in the pendular state. For the liquid bridges not to merge, the lling angle ϕ should be smaller than π/6. The volume of the liquid is approximately the volume of a cylinder with radius a and height h, excluding two small volumes v i . For contacting

112 Chapter 3 . 34 )

 112334 Isotropic compression: DEM study3.3.6 Dimensionless control parametersIn order to achieve greater generality, dimensionless parameters are used in the presentation of all simulation results. Some are immediately apparent in the contact law: friction coecient µ, damping coecient ζ. From numerical observations[START_REF] Radjaï | Discrete-element Modeling of Granular Materials[END_REF], it is known that, while the inuence of µ on the macroscopic mechanical behavior of model granular materials is always quite important, ζ is rather irrelevant in the quasistatic regime. Under a pressure P , the typical contact deection |h|, relative to the grain diameter d, is related to a dimensionless stiness parameter, which is conveniently dened as It is such that typical values of |h|/d are of order κ -1[START_REF] Agnolin | Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks[END_REF]. If a linear contact elasticity is used instead of the Hertz law, then one may simply use κ = K N /(dP ). In all simulations carried out here, κ is kept large enough along the whole compression curve that the elastic deections are kept very small (κ ≥ 4000). Such κ values hardly aect coordination numbers (due to contact recruitment caused by shrinking intergranular distances)[START_REF] Agnolin | Internal states of model isotropic granular packings. ii. compression and pressure cycles[END_REF].To those parameters ruling the behavior of dry grains, one should add three new dimensionless groups which become relevant in the context of the present study, involving attractive forces.With capillary forces, an additional dimensionless control parameter associated with the force law is the ratio of meniscus volume V m to d 3 , determining the range of the corresponding capillary attraction, and consequently the number of liquid bonds. Fig. 3.6 compares function F cap (h) according to Maugis model and to the Soulié formula.

Figure 3 . 6 :

 36 Figure 3.6: Capillary force law F cap (h), for two dierent meniscus values, according to the Maugis approximation and the Soulié formula.

Figure 3 . 7 :

 37 Figure 3.7: Monodisperse system on FCC lattice (Φ 0 = 0.3): (a)whole specimen; (b) small cube extracted from (a); (c) a small cube of maximal density (Φ 0 = π 3 √ 2

3 . 3 . 2 .

 332 Typical resulting disordered congurations are shown in Fig. 3.8. The same aggregation procedure applies to polydisperse systems.

3. 5

 5 .1 A reference case 3.5.1.1 Macroscopic behavior The reference test is rst run with a low initial solid fraction, Φ 0 = 0.30. In the present works, the results are shown in the form of the conventional compression curve showing void index e versus log P * . The compression and decompression curves in the reference case is shown in Fig. 3.9.

Figure 3 . 9 :( 3 . 38 )

 39338 Figure 3.9: (a) e versus log P * for reference system in compression cycle (with schematic view of three regimes in the inset); (b) comparison of reference case and typical noncohesive behavior, with pressure in kPa for glass beads of diameter 0.118 mm wet by water.

120 Chapter 3 .

 1203 Isotropic compression: DEM study

Fig. 3 .Figure 3 .

 33 Fig. 3.11(a) plots z c , z d , and z versus P * in the compression cycle. Initially, one has z d = 0, as the velocities in the aggregation process, in the reference case, are low, and do not allow contact opening (let us recall that distant interaction through menisci without intergranular contact only appear when contacts open menisci do not spontaneously form between non-contacting grains). The obtaining of z c 4 at low P * is consistent with a minimally connected initial network, which does not reorganize or fold onto itself as aggregation proceeds. The value z c = 4 is the one corresponding to isostaticity in a pack of frictional beads (although some limited simultaneous amount of force indeterminacy

1 Figure 3 . 12 :

 1312 Figure 3.12: Contact networks in loading (a, b, and c) and unloading (d, e, and f ) cycle.

  3.12, a and b). Numbers of distant interacting pairs are slightly decreased till the end of loading path (Fig.3.12(c)). The density of distant 3.5. Material behavior in isotropic compression 123 interacting pairs reaches its maximum near P * = 10 1 , which corresponds to the end of Regime II. In unloading path, when the system is gradually opened, the numbers of distant interacting pairs are slightly increased. It remains nearly constant at the end of unloading path (beyond the stage P * = 10 -1 . We can see an insignicant change in the density of the blue lines in Figs.[START_REF] Atkinson | An Introduction to Mechanics of Soils and Foundations: Through Critical State Soil Mechanics[END_REF].12(d, e, and f ).

2 Figure 3 . 13 :

 2313 Figure 3.13: Networks of compressive contacts (red lines), tensile ones (green lines), and contacts with F N = 0 (magenta). Line thicknesses are proportional to normal force intensity |F N |. Loading path: (a), (b), (c) and (d); unloading path: (e) and (f ).

4 Figure 3 . 14 :Figure 3 . 15 :

 4314315 Figure 3.14: Distribution of normal forces for dierent values of P * from the compression path, (a) normalized by the maximum tensile force F 0 , and (b) normalized by the average normal force F N . Red arrow indicates increasing P * .

Fig. 3 . 15 . 126 Chapter 3 .

 3151263 Fig.3.15. These results are quite similar to the ones obtained in other numerical studies on cohesive granular materials[START_REF] Richefeu | Stress transmission in wet granular materials[END_REF][START_REF] Gilabert | Computer simulation of model cohesive powders: Inuence of assembling procedure and contact laws on low consolidation states[END_REF][START_REF] Gilabert | Computer simulation of model cohesive powders: Plastic consolidation, structural changes, and elasticity under isotropic states[END_REF]

Figure 3 . 16 :Figure 3 . 17 :Fig. 3 .

 3163173 Figure 3.16: Void ratio (a), coordination numbers z c and z d (b) versus P * of the reference case for dierent values of meniscus V m /d 3 . Red arrow shows direction of increasing V m values.

128 Chapter 3 . 12 P * =10 1 Figure 3 . 18 :

 1283121318 Figure 3.18: Distribution of the number of contacts per grain of the reference case for dierent values of V m /d 3 at three typical states of P * . Red arrow orders curves by increasing V m .

. 4 . 4 Figure 3 . 19 :

 44319 Figure 3.19: Evolution of coordination number of close neighbors z(h) versus dimensionless interparticle distance h/d for dierent values of P * in loading path.

Figure 3 . 20 :

 320 Figure 3.20: Void ratio versus pressure P , as cohesive forces are suppressed at the beginning of unloading, starting at dierent points on the primary compression curve (case with initial solid fraction Φ 0 = 0.45).

Figure 3 .

 3 Figure 3.21: (a) Void ratio e, and (b) coordination numbers z c and z d versus P * for dierent values of initial solid fraction Φ 0 . Red arrow denotes increasing Φ 0 .

1 Figure 3 . 22 :

 1322 Figure 3.22: Coordination numbers z + and z -versus P * for dierent values of Φ 0 .

Fig. 3 . 12 P * =10 1 Figure 3 . 23 :

 3121323 Fig.3.24 shows the distributions of normal forces for dierent values of initial solid fraction Φ 0 at P * = 10 -2 . At the lowest pressure P * = 10 -3 , the P.D.Fs are nearly symmetric about 0 and very similar to one another. At P * = 10 -2 , these P.D.Fs are also roughly symmetric about 0, except the P.D.F of case Φ 0 = 0.25, because it already transited to the Regime II, as indicated by the extension of the positive wing to 1.5 and the nite value in -1. The delayed approach to Regime II for larger densities causes narrower P.D.Fs (as indicated according to red arrow in this gure). For increasing pressure, these P.D.Fs widen and undergo similar evolutions, approaching the characteristic force PDF shape of cohesionless systems already described in Sec. 3.5.1.2.

Figure 3 . 24 :

 324 Figure 3.24: Distribution of normal forces for dierent values of Φ 0 under low pressure P * = 10 -2 , is normalized by the F 0 . Red arrow denotes increasing Φ 0 .

Figure 3 .

 3 Figure 3.25: (a) Compression and decompression curves, and (b) coordination numbers z c and z d for dierent values of V 0 /V * . Red arrow denotes increasing V 0 /V * .

Fig. 3 . 136 Chapter 3 . 1 Figure 3 . 26 :Fig. 3 .

 3136313263 Fig.3.[START_REF] Brewster | Plug ow and the breakdown of bagnold scaling in cohesive granular ows[END_REF] shows that this parameter mainly aects the beginning of the compression curve. Six dierent values of V 0 /V * are used, and the dierence between minimum and maximum values is two-hundredfold. Under low P * , higher initial agitation velocities entail wider Regime I plateaus. This corresponds to the intervals of reduced pressure from 10 -3 ≤ P * ≤ 10 -2 to 10 -3 ≤ P * ≤ 6 × 10 -2 . If initial velocity V 0 is of the order of V * , the initial packing structure is strongly perturbed. z c values in the initial state depart then from the minimum, one observes z c > 4, reaching z c = 4.75 for the largest V 0 . z d is also no longer equal to zero in the initial states. In other words, the larger the velocities, the stronger the initial structures (see Fig. 3.25(a)). More contacts are created, larger initial coordination numbers z c and z d are observed at low P * (see Fig. 3.25(b)). Conversely, with low agitation velocities (V 0 ≤ V * ), grains gently stick to one another and clusters of aggregated grains, once formed, are not disturbed. Lower

3. 5 . 137 (Figure 3 . 138 Chapter 3 . 12 P * =10 0 Figure 3 . 28 :

 513731383120328 Figure 3.27: A 3D slice of specimen for dierent values of V 0 /V * at the beginning of loading path (P * = 10 -3 ).

Figure 3 . 29 :

 329 Figure 3.29: Distribution of normal forces of reference case for dierent values of V 0 /V * under low pressure P * = 10 -3 , normalized by F 0 . Red arrow denotes increasing V 0 /V * .

3. 5 .Figure 3 . 140 Chapter 3 .Figure 3 . 31 :

 531403331 Figure 3.30: (a) Void ratio e and (b) coordination numbers z c and z d versus P * for dierent values of V 0 /V * . Red arrow denotes increasing V 0 /V * .

Figure 3 . 32 :

 332 Figure 3.32: Grains contained in a slice (thickness 3d min ), at low P * , for dierent values of V 0 /V * .

Fig. 3 .

 3 33 directly confronts the variations of void ratio and of coordination numbers in the compression cycle for those two dierent materials, the monodisperse (MDS) and the polydisperse (PDS) ones. This comparison reveals quantitatively close behaviors. Slight dierences can be Chapter 3. Isotropic compression: DEM study

Figure 3 .

 3 Figure 3.33: (a) Void ratio e and (b) coordination numbers z c and z d versus P * for MDS and PDS.

Figure 3 .

 3 Figure 3.34: (a) Void ratio e and (b) coordination numbers z c and z d versus P * in compression cycle for dierent values of Φ 0 . Red arrow denotes increasing Φ 0 .

1 Figure 3 . 35 :

 1335 Figure 3.35: Coordination numbers z + and z -versus reduced pressure P * for dierent values of Φ 0 . Red arrow denotes increasing Φ 0 .

Figure 3 .

 3 Figure 3.36: (a) Void ratio e and (b) coordination numbers z c and z d versus P * in compression cycle for dierent values of V 0 /V * . Red arrow indicates increasing V 0 /V * .

1 Figure 3 . 37 :

 1337 Figure 3.37: Coordination numbers z + and z -versus P * for dierent values of V 0 /V * (increasing according to red arrow).

Fig. 3 .Figure 3 . 38 : 3 12 P * =10 1 Figure 3 . 39 :Figure 3 . 40 :

 33383121339340 Fig.3.37 displays the variations of z + and z -in the compression cycle for six values of initial agitation level. z + and z -, just like z c = z + + z -, strongly depend on V 0 in the beginning of the compression cycle, and increase in parallel with V 0 , from z -= 0.82 and z + = 1.30 to z -= 1.95 and z + = 2.30. As Regime III is approached (P * > 10 0 ), the system becomes denser, the tensile contacts gradually vanish while z + approach its maximum. Part of this evolution is reversed along the decompression path, with,

150 Chapter 3 .Figure 3 .

 15033 Figure 3.41: (a) Void ratio e and (b) coordination numbers z c and z d versus P * for dierent values of µ R = µ P . Red arrow denotes increasing µ R /d.

3. 6 .Figure 3 . 42 :Fig. 3 . 152 Chapter 3 .

 634231523 Figure 3.42: Distribution of normal forces for dierent values of µ R /d at P * = 10 -2 . Red arrow denotes increasing µ R /d.

Figure 3 .

 3 Figure 3.43: (a) Void ratio e and (b) coordination numbers z c and z d versus P * in compression cycle for systems both without RPR and with RPR.

  The numerical study of isotropic compression of initially loose wet bead systems, forming continuous, random aggregates stabilized by capillary cohesion, has been carried out for dierent model materials: monodisperse or not, possibly with some resistance to rolling and pivoting in the contacts, and with varying values of initial density, coordination and connectivity (the latter features being determined by the agitation intensity in the aggregation stage), and meniscus volume. The plastic behavior of such a model wet granular materials is characterized at both macroscopic and microscopic levels. The main conclusions of this chapter are summarized as follows.

Contents 4 . 1 156 Chapter 4 .

 411564 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 4.2 Simulation of oedometric compression tests . . . . . . . . . . 156 4.2.1 Compression of dierent initial states: initial V 0 and coordination.157 4.2.2 Inuence of initial solid fraction . . . . . . . . . . . . . . . . . . 162 4.3 Oedometric versus isotropic compression . . . . . . . . . . . 166 4.4 Experiment versus simulation . . . . . . . . . . . . . . . . . . 168 4.4.1 Compaction behavior . . . . . . . . . . . . . . . . . . . . . . . . 168 4.4.2 Geometrical structure change . . . . . . . . . . . . . . . . . . . 170 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 Oedometric compression: DEM study and experimental confrontation 4.1 Introduction

4. 2 . 1

 21 Compression of dierent initial states: initial V 0 and coordination.

Figure 4 . 1 :

 41 Figure 4.1: (a) Void ratio e, and (b) coordination numbers z c and z d versus σ * 1 for

Figure 4 . 2 :

 42 Figure 4.2: Coordination numbers of z + (compressive contacts) and z -(tensile contacts) versus reduced vertical stress σ * 1 for dierent values of V 0 /V * , increasing as shown by read arrow.

Figure 4 . 3 :

 43 Figure 4.3: Coecient of lateral pressure K 0 versus reduced vertical stress σ * 1 for dierent

4

 4

. 4 .Figure 4 . 4 : 1 3 1 7 × 10 -2 and σ * 1 10 - 1 [

 444111011 Figure 4.4: Contribution of distant interactions to normal stress σ 22 . Red arrow denotes increasing V 0 /V * .

Fig. 4 . 161 (

 4161 Fig. 4.5 displays the contributions of tangential forces and capillary forces to the normal stress σ 22 . Remarkably, the contribution of tangential stress is almost negative

Figure 4 . 5 :

 45 Figure 4.5: (a) Contributions of tangential (a) and capillary (b) interactions to normal stress σ 22 . Red arrow denotes increasing V 0 /V * .

Figure 4 . 8 :

 48 Figure 4.8: Coecient of lateral pressure K 0 versus reduced axial stress σ * 1 for dierent

Figure 4 . 9 :

 49 Figure 4.9: Contribution of distant interactions to normal stress σ 22 for dierent values of Φ 0 . Red arrow denotes increasing Φ 0 .

Figure 4 . 166 Chapter 4 .

 41664 Figure 4.10: (a) Contributions of tangential (a) and capillary (b) interactions to normal stress σ 22 for dierent values of Φ 0 . Red arrow denotes increasing Φ 0 .

Figure 4 .

 4 Figure 4.11: (a) Void ratio e, and (b) coordination numbers z c and z d versus σ * m for the minimum and maximum V 0 /V * in both simulations of oedometric and isotropic tests.Red arrow denotes increasing V 0 /V * .

Figure 4 . 12 :

 412 Figure 4.12: Distribution of normal forces for both oedometric and isotropic models at low vertical stress, normalized by the maximal tensile force F 0 . Red arrow denotes the increasing of stress or pressure.

F 14 in Chapter 3 . 4 . 4

 14344 N reaches -F 0 signals the onset of in Regime II collapse (red curve in Fig. 4.11(a)). Meanwhile, a distribution bounded by |F N /F 0 | ≤ 0.5 in the oedometric case (magenta curve) implies that the initial network may support further stress increments without rearranging. This explains why the plateau of Regime I in the oedometric case is wider than the one in the isotropic case (Fig. 4.11(a)). Under high axial stress, the P.D.Fs of normal forces F N , normalized by the average normal force F N , for both models, are quite similar, as shown in Fig. 3.Experiment versus simulation This section aims to confront the results of the experiments of Chapter 2 and the ones of the simulated oedometric compression test in Sec. 4.2. The comparison is carried out for the macroscopic compression curve and for microstructural changes. All stresses are given here in kPa. 4.4.1 Compaction behavior A comparison of compression behavior between experiment and simulation of oedometric compression tests for four values of initial solid fraction Φ 0 (Φ 0 = 0.30, Φ 0 = 0.32, Φ 0 = 0.35, and Φ 0 = 0.40) is shown in Fig. 4.13. In general, both experimental and simulated curves have the same S-shape as already mentioned in Chapter 2. Three regimes can be equally identied in the compression curves, depending on the range of vertical stress σ 1 .

Fig. 4 . 1 ≤ 4 . 4 .

 4144 Fig.[START_REF] Tripodi | Elastoplastic nite element model development and validation for low pressure uniaxial compaction of dry cohesive powders[END_REF].13(a) shows the experimental and simulated compression curves for the initial solid fraction Φ 0 = 0.30 in which the dashed curve (CT01) was obtained with the cell φ20. The global features of those curves are in agreement. Some dierences are notable in the stress corresponding to the onset of plastic collapse, which is larger by a factor of 2 or 3 in the laboratory (whence a lateral shift in comparison to the numerical results); in the compression index (slope of the fastly decreasing part of the curves in Regime II), which is somewhat smaller in experimental results; and nal void ratio, which is larger in the numerical case.It is well recognized that Regime I in experimental curves exists in a very short range of vertical stress. The plateau of Regime I, if it exists, may be at the very low vertical stress (σ 1 < 10 -2 kPa). Under low stress, the initial structure of experiment is slightly collapsed while the one of simulation can sustain the increasing stress without rearrangement, with a void ratio e remaining almost constant. In Regime II, under growing of vertical stress, the collapse in simulation happens earlier in comparison with the experimental curves. This is conrmed by the slope of three curves in Fig.4.13(a).For instance, Regime II of simulated curve exists in interval 4 × 10 -2 kPa ≤ σ 1 ≤

Figure 4 . 13 :

 413 Figure 4.13: Oedometric compression curves in experiments and simulations (no RPR) for dierent values of Φ 0 . (a) Φ 0 = 0.30, (b) Φ 0 = 0.32, (c) Φ 0 = 0.35, (d) Φ 0 = 0.40.

Fig. 4 .

 4 Fig. 4.14 plots the coordination numbers of neighbours, z(h), versus axial stress σ 1 , and the total coordination number z in the local scan observations (see Sec. 2.5 of Chapter 2), versus positions of scans along the loading path. The total coordination number including lost contacts z lc exp is also plotted (brown triangle point). The results of the case F30OV1Vm3 is only displayed in the loading path with ve values of h,

denoted as h 1

 1 to h 5 , with h 1 /d = 0, h 2 /d = 0.025, h 3 /d = 0.05, h 4 /d = 0.075, and h 5 /d = 0.1. As mentioned in Chapter 3, when h/d = 0, z(h) is equal to z c , the contact coordination number (see Sec. 3.5.1 in Chapter 3). The positions of scans are directly interpolated from the compression curves of Fig.2.11 in Chapter 2. In general, the order of magnitude and the general trend of the z exp is well captured by the simulation. At σ 1 ≈ 4 × 10 -2 kPa, z exp = 4.75 is remarkably higher than all of the z(h). This could be explained by the stability of Regime I in the simulated case, where void ratio e is 4

Figure 4 . 14 :

 414 Figure 4.14: Coordination numbers of neighbors z(h) versus vertical stress σ 1 and the total coordination number z versus positions of scans.

Chapter 4 .

 4 Oedometric compression: DEM study and experimental confrontationSeveral conclusions can be drawn as follows.

  .34 Grain size distribution of the 3D reconstructed specimens . . . . . . . . . 2.35 Example of 3D reconstructed specimen . . . . . . . . . . . . . . . . . . . 2.36 Number of particles of the ten SVs for four local scans after detection process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.37 Macroscopic solid fraction (triangle) and the average solid fraction of 4 local scans before and after the detection process . . . . . . . . . . . . . .43 Contact network of a typical SV under growing of compaction . . . . . . 2.44 3d-slices of the typical SV under growing of compaction . . . . . . . . . .

. . . . . 2.26 3D reconstructed image of four local scans . . . . . . . . . . . . . . . . . 2.27 The algorithm to detect the spherical structures . . . . . . . . . . . . . . 2.28 Construction of the accumulation array from the gradient eld . . . . . . 2.29 Denition of the signature curve . . . . . . . . . . . . . . . . . . . . . . . 2.30 An example of a signature curve computed from the 3D image . . . . . . 2.31 A slice in a 3D detected image . . . . . . . . . . . . . . . . . . . . . . . . 2.32 Method to nd the lost contacts, inside cube is the standard volume; outside cube is the extended volume . . . . . . . . . . . . . . . . . . . . . 2.33 Example of 3D reconstructed cubes of four local scans. Left-hand side: before detection, right-hand side: after detection . . . . . . . . . . . . . . 22.38 Average total coordination numbers of 4 local scan with three values of tolerance (0, 1 and 2 voxels) . . . . . . . . . . . . . . . . . . . . . . . . . 2.39 Coordination numbers of close neighbors for four local scans . . . . . . . 2.40 Number of contacts per grain for three tolerances of 0, 1 and 2 voxels . . 2.41 An example of a sphere with one contact . . . . . . . . . . . . . . . . . . 2.42 Radial distribution functions for particle centers of four local scans . . . xvi List of Figures 2

  1.1 Schematic diagrams for classication of wet granular materials with various amount of liquid . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

			13
	1.2	Examples of imaging techniques to porous media systems . . . . . . . . .	17
	1.3	Combining experimental observations with computer simulations . . . . .	42
	2.1	Initial parameters of selected specimens after preparation. The sample

CT03 was prepared for the XRCT tests . . . . . . . . . . . . . . . . . . . 57 2.2 Details of the scanning program. G: global scan, L: local scan . . . . . . 61 3.1 Initial parameters of the monodisperse system without RPR . . . . . . . 116 3.2 Initial parameters of the polydisperse system without RPR . . . . . . . . 116 3.3 Initial parameters of the monodisperse system with RPR . . . . . . . . . 117 3.4 Probability of meniscus formation between close neighbors versus P * . . 130 4.1 Initial parameters of the oedometric simulations (without RPR) . . . . . 156

  Numerical investigations . . . . . . . . . . . . . . . . . . . . . 28 1.3.1 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.3.2 Typical observations . . . . . . . . . . . . . . . . . . . . . . . . 30 .5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

	1.1. Mechanical behavior of granular materials	7
	1.1 Mechanical behavior of granular materials
	1.1.1 General considerations

Contents 1.1 Mechanical behavior of granular materials . . . . . . . . . . 7 1.1.1 General considerations . . . . . . . . . . . . . . . . . . . . . . . 7 1.1.2 Bonded granular materials . . . . . . . . . . . . . . . . . . . . . 10 1.1.3 Wet granular materials . . . . . . . . . . . . . . . . . . . . . . . 12 1.2 Experimental investigations . . . . . . . . . . . . . . . . . . . 15 1.2.1 Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.2.2 Typical observations . . . . . . . . . . . . . . . . . . . . . . . . 22 1.3 1.3.3 Some characteristics of wet granular material in DEM simulation 34 1.4 Experimental investigation versus numerical investigation . 39 1.4.1 Existing studies combining experimental and numerical methods 39 1.4.2 Combination between the XRCT and the DEM . . . . . . . . . 43 6 Chapter 1. Literature review 1

Table 1 .

 1 non-destructive methods have been widely applied to observe eectively the properties of granular materials, from macro-to microscopic level. In order to have a general view of non-destructive techniques, we present here a sketchy history of some recent advanced techniques (type of imaging techniques) and its applications to granular materials research, as given in Table1.2 (following from Al-Raoush & Willson[START_REF] Al-Raoush | Extraction of physically realistic pore network properties from three-dimensional synchrotron X-ray microtomography images of unconsolidated porous media systems[END_REF]). These 2: Examples of imaging techniques to porous media systems

	typical studies were carried out over the last thirty years (from 1982 to 2016). The
	techniques are briey presented and their application will be clearly illustrated in the
	following sections.

Aside from the conventional destructive techniques (e.g. shear box, one-dimensional compression test, biaxial compression test, triaxial compression test, etc.), more and Chapter 1. Literature review more

Table 1 .

 1 3: Combining experimental observations with computer simulations

	ber and packing distribution	Bed structure, coordination num-	Application
		Glass beads	Material
	method	Optimization	Simulation
	isotropy	XRCT Homogeneity	Experiment EXP-SIM link
		and
	and Alsaleh	Al-Raoush,	Reference
		Al-Raoush

Table 2

 2

.1. In this table, the tests named using prex OD are performed with cell φ50 and the tests named using prex CT are carried out with cell φ20. The water content of all specimens

Table 2 .

 2 1: Initial parameters of selected specimens after preparation. The sample CT03 was prepared for the XRCT tests. w ap was measured after testing.

				Parameters
	Tests	h f (mm)	n (%)	w ap (%)	e 0 (-)	S (%)
	OD01	80	70.19	2.71	2.35	1.91
	OD02	80	70.12	2.70	2.35	2.12
	OD03	120	66.56	2.00	1.99	1.91
	OD04	120	66.86	2.30	2.02	2.49
	OD05	150	63.72	2.50	1.76	2.49
	OD06	150	63.95	2.20	1.77	2.48
	OD07	260	60.60	1.80	1.54	1.68
	OD08	260	60.02	1.60	1.50	2.19
	CT01	80	70.05	2.89	2.34	1.92
	CT02	80	69.93	2.85	2.33	1.42
	CT03	80	69.67	-	2.30	-

-05x 2 -0.0158x + 3.4126 R² = 0.9971

Table 2 .

 2 Mechanical behavior under increasing vertical stressAll the tests were carried out on specimens with very loose initial structures, from e 0 = 2.35 to e 0 = 1.50.

	Scan	h i (mm)	e i (-)	Notation
	01	h 0 10	2.2963	S01G
	02	h 1 9.8	2.2303	S02G
	03	h 1 9.8	-	S03L
	04	h 2 8.7	1.8678	S04G
	05	h 3 8.5	1.8678	S05G
	06	h 3 8.5	-	S06L
	07	h 4 7.5	1.4722	S07G
	08	h 5 7.3	1.4063	S08G
	09	h 5 7.3	-	S09L
	10	h 6 6.4	1.1096	S10G
	11	h 7 6.2	1.0437	S11G
	12	h 7 6.2	-	S12L

2: Details of the scanning program. G: global scan, L: local scan.

Table 2 .

 2 From the 3D image stacks of all global scans, the grey level histograms were successively extracted, as shown in Fig.2.13. For 8-bit images, grey levels range from 0 to 255.Depending on compaction steps, the two peaks more or less visible on each histogram are progressively changed. The rst one on the left-hand side strongly decreases while the second one on the right-hand side increases. At the highest level of compaction, the left peaks of the dark-green curves (continuous and dashed lines) almost no longer exist.Meanwhile the peaks of these curves on the right-hand side reach the highest level. This conrms that the specimen approached a very dense state.

	2.

  Conclusion 95 collapses negligibly because of the elastic deformation of grains at contact regions. This is completely demonstrated by the highest density of S12L and the slight increases of z and z(h) from the tolerance of 1 voxel (see Figs.2.38 and 2.39

2.12 and 2.26). Besides, at higher densities, the number of contacts is larger than the one at lower densities (see

Figs. 2.38 and 2.39).

In the stage 3 of the compaction curve, starting from the S11G and S12L, the structure 2.7.

  .3.4 Liquid content and limit of pendular state . . . . . . . . . . . . 3.3.5 Average normal force . . . . . . . . . . . . . . . . . . . . . . . . 3.3.6 Dimensionless control parameters . . . . . . . . . . . . . . . . . 3.4 Simulation procedures for the compression test . . . . . . . Compression cycle . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Material behavior in isotropic compression . . . . . . . . . . 3.5.1 A reference case . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Inuence of drying or of saturating . . . . . . . . . . . . . . . . 3.5.3 Inuence of initial state . . . . . . . . . . . . . . . . . . . . . .

	3.4.1 Specimen preparation . . . . . . . . . . . . . . . . . . . . . . .
	3.4.2 Simulation parameters . . . . . . . . . . . . . . . . . . . . . . .
	3.4.3
	. 99
	3.2.2 Interaction law: normal forces . . . . . . . . . . . . . . . . . . . 100
	3.2.3 Tangential forces . . . . . . . . . . . . . . . . . . . . . . . . . . 104
	3.2.4 Resistance to rolling and pivoting (RPR) . . . . . . . . . . . . 105

Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 3.2 Model material and interaction laws . . . . . . . . . . . . . . 99 3.2.1 Particle shape and size distribution . . . . . . . . . . . . . . . 3.3 Basic denitions and properties . . . . . . . . . . . . . . . . . 106 3.3.1 DEM equations, boundary conditions and stress control . . . . 106 3.3.2 Equilibrium condition . . . . . . . . . . . . . . . . . . . . . . . 108 3.3.3 Coordination numbers . . . . . . . . . . . . . . . . . . . . . . . 109 98 Chapter 3. Isotropic compression: DEM study 33.5.4 Eects of polydispersity . . . . . . . . . . . . . . . . . . . . . . 3.6 The case of wet beads with RPR . . . . . . . . . . . . . . . . 3.6.1 Inuence of initial solid fraction . . . . . . . . . . . . . . . . . . 3.6.2 Inuence of initial agitation intensity . . . . . . . . . . . . . . . 3.6.3 Inuence of rolling and pivoting friction coecients . . . . . . . 3.6.4 Comparison of bead assemblies with and without RPR . . . . . 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  The actual value of d is not relevant for our results, if they are expressed in dimensionless form. However, in order to provide physical orders of magnitude, results will also often be reported in dimensional form (e.g. pressures in kPa), in which case they pertain to glass beads of diameter d = 0.118 mm.

	100			Chapter 3. Isotropic compression: DEM study
	Moreover, we mainly consider equal-sized, identical spherical beads of diameter d.
	The inuence of a (moderate) grain polydispersity will be assessed on dealing with
	the following probability distribution function (P.D.F) for the distribution of diameter
	d, varying between a minimum value d min and a maximum value d max :
				p(d) =	2d 2 min d 2 max d 2 max -d 2 min	1 d 3 .	(3.1)
	This distribution is thus uniform by volume, and its rst moments are given by:
	d =	2d min d max d max -d min	;	d 2 =	2d 2 min d 2 max d 2 max -d 2 min	ln	d max d min	;	d 3 =	2d 2 min d 2

max d max + d min .

Table 3 .

 3 3: Initial parameters of the monodisperse system with RPR

  dAll assembled samples of the series of Tables 3.1, 3.2 and 3.3 are compressed under growing external isotropic pressure. A stepwise pressure-controlled loading path is applied.

	118			Chapter 3. Isotropic compression: DEM study
	3.4.3 Compression cycle
							Notation
	0.25	F25	0.2041	V1	1.00 × 10 -3 Vm3 0.02	R2	F25R2V1Vm3
	0.30	F30	0.2041	V1	1.00 × 10 -3 Vm3 0.02	R2	F30R2V1Vm3
			0.2041	V1	0.004	R1	F30R1V2Vm3
			0.2041	V1	0.1	R3	F30R3V2Vm3
			0.2041	V1	0.5	R4	F30R4V2Vm3
			0.4082	V2	0.02	R2	F30R2V2Vm3
			1.2247	V3	0.02		F30R2V3Vm3
			4.0825	V4	0.02		F30R2V4Vm3
			12.2474	V5	0.02		F30R2V5Vm3
			40.8248	V6	0.02		F30R2V6Vm3
	0.32	F32	0.2041	V1	1.00 × 10 -3 Vm3 0.02	R2	F32R2V1Vm3
	0.35	F35	0.2041	V1	1.00 × 10 -3 Vm3 0.02	R2	F35R2V1Vm3
	0.40	F40	0.2041	V1	1.00 × 10 -3 Vm3 0.02	R2	F40R2V1Vm3
	0.45	F45	0.2041	V1	1.00 × 10 -3 Vm3 0.02	R2	F45R2V1Vm3

Table 3 .

 3 4: Probability of meniscus formation between close neighbors versus P * . 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4

	P *		z d	F30V1Vm3 z d /[z(D r ) -z c ]	z d	F30V6Vm3 z d /[z(D r ) -z c ]
	10 -3	0.0000	0.0000	0.2510	0.3990
	10 -2	0.0005	0.0014	0.2510	0.3990
	10 -1	0.0515	0.1102	0.2990	0.4262
	10 0	0.9970	0.6115	1.1040	0.6404
	10 1	1.7720	0.7018	1.7990	0.6909
	10 2	1.6830	0.6826	1.8020	0.6875
	10 3	1.4540	0.6462	1.5590	0.6509
	10 4	1.0760	0.5565	1.1680	0.5713
	3.5.2 Inuence of drying or of saturating
	Void ratio, e (-)	0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 10 -3 A B′ C′	D′ E′ F′ G′	C Pressure, P (kPa) D E	F	B G

  The polydisperse system (without RPR) is studied with the diameter distribution, such that d max = 2d min , presented in Sec. 3.2.1, and the values of initial parameters given in Table3.2. The main objective of this study is to assess whether and to what extent

Table 4 .

 4 1: Initial parameters of the oedometric simulations (without RPR)

  The denition of a reduced stress used for pressure in Eq. (3.35) of Chapter 3, is applied, in the oedometric case, to the axial stress, compared to the tensile contact strength:

	4.2. Simulation of oedometric compression tests	157
					3	Notation
	0.30	F30O	0.2041	V1	1.00 × 10 -3 Vm3 F30OV1Vm3
			0.4082	V2		F30OV2Vm3
			1.2247	V3		F30OV3Vm3
			4.0825	V4		F30OV4Vm3
			12.2474	V5		F30OV5Vm3
			40.8248	V6		F30OV6Vm3
	0.32	F32O	0.0241	V1	1.00 × 10 -3 Vm3 F32OV1Vm3
	0.35	F35O				F35OV1Vm3
	0.40	F40O				F40OV1Vm3

  1 . At high 4.2. Simulation of oedometric compression tests 159 pressure, in Regime III, K 0 remains almost constant (K 0 ≈ 0.8) whatever the initial connectivity. In the rst part of the unloading path, K 0 slightly increases (from K 0 ≈ 0.7 to K 0 ≈ 0.8) with similar values in all cases (the inuence of the initial state seems now forgotten). K 0 starts to diverge in the latter half of unloading path, when the vertical stress σ * 1 < 10, reaching values of several tens or several hundreds. This behavior is simply due to the existence of residual lateral stresses of order F 0 /d 2 , while the axial stress is requested to decrease to much smaller values: σ * 1 F 0 /d 2 , by denition of σ *

				1 . The
	dierent residual lateral stresses observed for dierent initial coordination (or initial ag-
	itation, V1 to V6) correlates with the dierent density of tensile contacts (see Fig. 4.2)
	at the end of the unloading curve in the six specimens.
	10 -3 10 -1 10 0 10 1 10 2 Coefficent of lateral pressure, K 0 (-)	10 -2	10 -1 Reduced vertical stress, σ * 1 (-) 10 0 10 1 10 2 F30OV1Vm3 F30OV2Vm3 F30OV3Vm3 F30OV4Vm3 F30OV5Vm3 F30OV6Vm3	10 3

The scholarship was nanced by the Ministry of Education and Training of Vietnam (project 911). Ce travail a été réalisé au

Chapter 3. Isotropic compression: DEM study

Tangential forces

The tangential elastic force F T ij , as in [230], relates to the relative tangential displacement in the contact, δu T ij , involving a (deection-dependent) tangential stiness coecient K T assumed proportional to K N :

(3.15)

Tangential stiness K T has to be suitably adapted (rescaled) whenever the normal elastic force decreases, in order to avoid spurious elastic energy creation [START_REF] Elata | Contact force-displacement laws and the mechanical behavior of random packs of identical spheres[END_REF][START_REF] Agnolin | Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks[END_REF].

Let us recall that the relative displacement at the contact point combines the displacements of both grains, u i and u j , and their (small) rotations θ i and θ j , thus, with transparent notations for the radii, and n ij the unit normal vector, which for spherical objects points from the centre of i to the centre of j:

The normal component of relative displacement (3.16) is subtracted to dene the relative tangential displacement δu T ij appearing in the tangential contact law.

The Coulomb condition enforces inequality

Simulation parameters

The sample series and the parameter values for the monodisperse material without RPR are listed in Table 3.1. Sample F30V1Vm3 is chosen as the reference case, from which the dierent control parameters are varied one by one: initial solid fraction Φ 0 , initial agitation level V 0 /V * in the aggregation process, meniscus volume V m .

Polydisperse samples (without RPR) are prepared with six dierent initial agitation intensities, as listed in Table 3.2.

The sample series for the (monodisperse) material with RPR is given in Table 3.3.

The dimensionless rolling friction coecient, µ R /d, and its pivoting friction counterpart µ P /d are both set to 0.02 for most simulations.

Chapter 3. Isotropic compression: DEM study stress, the compressive and tensile contacts are approximately balanced. z -and z + uctuate around 2. Under growing of vertical stress, the population of compressive contacts gradually increases while the population of tensile ones is steadily depleted until the end of loading path. z -and z + vary in opposite directions upon unloading, the former decreasing and the latter increasing. At the end of unloading, both z -and z + remain almost constant, z + being larger than z -. The variations of z -and z + are similar with the ones obtained for the isotropic system (see Fig. Concerning the coecient of lateral pressure K 0 , its variation during the compression and decompression cycles for dierent values of Φ 0 is illustrated in Fig. 4.8. Initially, K 0 equals 0.4. In other words, the horizontal stress is approximately 40% of the vertical stress in all cases. K 0 is then slightly decreased in Regime I in the stress interval 10 -3 ≤ σ * 1 ≤ 10 -2 (to 3×10 -2 ). During the collapse of all specimens, K 0 is signicantly increased and it approaches to 1 at the end of the collapse (at σ * 1 = 10 1 ). Subsequently, in Regime III (σ * 1 ≥ 10 1 ), K 0 is negligibly decreased. At the end of loading path, the horizontal stress is approximately 80% of the vertical one (K 0 ≈ 0.8). Upon unloading, K 0 is slightly increased until the vertical stress decreases down to σ * 1 = 10 0 . Since then, K 0 strongly increases when the vertical stress decreases to very low values. In other words, the horizontal stress is higher than vertical one, from sevenfold to thirty-fold. This could be explained by the increase in the number of tensile contacts z -(see Fig. 4.7). Besides, the volumetric deformation of each specimen also contributes to the dierence of K 0 when the vertical stress is very low. In other words, the horizontal stress is directly proportional to the specimen's volume while the vertical one is inversely proportional to it. In our simulation, the denser the systems, the smaller the total volumes. That