Keywords: milieux granulaires, thermographie infrarouge, dynamique moléculaire, texture granulaire, thermoelastic stress analysis, contrainte hydrostatique granular materials, infrared thermography, molecular dynamics, granular texture, thermoelastic stress analysis, hydrostatic stress

L'objectif de la thèse est d'analyser le comportement mécanique de matériaux granulaires composites bidimensionels en terme de textures granulaires en utilisant deux approches : étude expérimentale par "thermoelastic stress analysis" et étude numérique par dynamique moléculaire. Les systèmes granulaires composites sont préparés à l'aide de cylindres en polyoxyméthylène (POM) et polyéthylène haute densité (PEHD), présentant un rapport de rigidité de 4 entre eux. Différents rapports de diamètres et de nombres de particules sont considérés. Les résultats expérimentaux et numériques sont en bon accord à l'échelle macroscopique. En particulier, le réseau fort (qui est ici caractérisé par des contraintes hydrostatiques supérieures à la valeur moyenne) contient moins de 50% des particules, et présente une distribution décroissance exponentielle quel que soit le type de particules considéré pour l'analyse (particules souples, particules rigides, toutes les particules). De plus, la distribution des contacts entre particules rigides (contacts POM-POM) est anisotrope et tend à s'organiser dans le sens de la direction du chargement extérieur appliqué, tandis que les autres types de contact agissent principalement pour maintenir le système en équilibre.

To

LIST OF TABLES

Page Table 1.1 Values of physical constants of interest in thermal radiation theory [START_REF] Maldague | Theory and practice of infrared technology for nondestructive testing[END_REF] Table 1.2 Technical specification of Cedip Jade III-MWIR infrared system [START_REF]DC001U-B Jade User Manual[END_REF] Table 2.1 Contact stiffness of POM for each configuration (in kN/mm) compression of photoelastic disks [START_REF] Majmudar | Contact force measurements and stress-induced anisotropy in granular materials[END_REF], b) experimental setup of the carbon paper method [START_REF] Mueth | Force distribution in a granular medium[END_REF], and c) experimental setup of the developed carbon paper method by using force sensors [START_REF] Løvoll | Force measurements on static granular materials[END_REF] 16

LIST OF FIGURES

Figure 1.6 Probability distribution of normalized contact forces in semilogarithmic scale [START_REF] Radjaï | Force distributions in dense two-dimensional granular systems[END_REF]. Four different 2D samples in term of particle size distributions and particle numbers are simulated under a uniaxial compression with size span s = 0.96 [START_REF] Voivret | Multiscale force networks in highly polydisperse granular media[END_REF]. The white particles are excluded from the force-bearing network. The amplitude of normal force is represented by the line thickness. b) Probability density functions of normal forces for various sizes span s of particle diameter [START_REF] Voivret | Multiscale force networks in highly polydisperse granular media[END_REF] Figure 1.9 a) Contact network of normal forces in a strongly cohesive granular system simulated by molecular dynamics method [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF].

The red lines represent compressive efforts, while the tensile efforts are presented the green lines. These lines are proportional to the amplitude of normal force. b) Semi-logarithmic plot of probability distribution function P of normalized normal forces for three different friction coefficients in the strongly cohesive system [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF] Figure 1.10 Schematic for definition of contact directions θ deposition under gravitational force [START_REF] Radjaï | Micro-statistical features of cohesionless granular media[END_REF], b) packing under isotropic compression [START_REF] Liu | Microscopic study on stress-strain relation of granular materials[END_REF], and c) packing under biaxial compression [START_REF] Liu | Microscopic study on stress-strain relation of granular materials[END_REF] Figure 1.12 a) Amplitude of mechanical anisotropy as a function of normalized contact forces [START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]. b) Angular distributions of contact forces for 2D granular packing obtained from numerical simulations [START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF] Figure 1.13 Sketch of molecular dynamics algorithm Figure 1.14 Geometry of contact between particle i with radius r i and particle j with radius r j regularized form [START_REF] Cambou | Micromechanics of Granular Materials[END_REF] Figure 1.17 Experimental study of Herschel [START_REF] Herschel | Experiments on the refrangibility of the invisible rays of the sun[END_REF]: A, B are the small stand, C, D are the prism at window, E is the spectrum of sunlight on the table, and 1, 2, 3 are the thermometers on the stand Figure 1.18 Electromagnetic spectrums with wavelength unit in micrometers [START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF] Figure 1.19 Spectral hemispherical emissive power of black body in vacuum for several absolute temperatures as a function of the wavelength [START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF] Figure 1.20 Incident radiations that impinge on a body [START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF] Figure 1.21 Transmission of a 6000 ft horizontal thick atmospheric layer at sea level containing precipitate water of 17 mm [START_REF] Rogalski | Infrared devices and techniques[END_REF] Figure 1.22 Ordinary components of infrared thermography systems

GENERAL INTRODUCTION

What are granular materials? This is a classic question that was often asked to the author. During an international congress in Thailand, the author was very surprised that no one knows precisely what granular materials are. Even though granular materials are plentiful around us, most of people do not even perceive them. The following examples present granular materials that we can encounter in our everyday lives: you will see coffee beans inside the coffee grinder when you go to the cafe shop, your feet will touch sand when you walk on the beach, you will find stone under the railroad when you travel by train, sugar is required when you would like to prepare your food to be sweet, etc. In this study, granular materials are defined as a collection of solid particles whose mechanical behavior at macroscopic level is governed by the interaction forces between the particles. They are composed of grains with distributions in sizes and shapes of particles, and can be made of different types of constitutive materials. Granular materials display complex behaviors which are distinguishable from those of common solids, liquids, and gases [START_REF] Jaeger | Physics of the granular state[END_REF][START_REF] Jaeger | Granular solids, liquids, and gases[END_REF]. For instance, a sand pile behaves like a solid when the angle of slope of the pile is lower than a limit angle which referred to as the angle of repose [START_REF] Jaeger | Physics of the granular state[END_REF][START_REF] Jaeger | Granular solids, liquids, and gases[END_REF][3]. When the slope exceeds the angle of repose, grains start to flow and a fluid-like behavior is observed.

It could be said that the industrial revolution was a beginning point for the development of sciences and technologies, particularly in the developed countries such as England, America, France, Germany, or even Japan. Such development required personnel with knowledge and skill. This was one of main reasons why the number of engineers and scientists increased in these countries. Moreover, it is undeniable that the development of science and technology mainly provides for human well-being. This ambition leads all the countries to construct and develop their infrastructure in the fields of transportation, energy supply, irrigation, telecommunication, and waste management.

The structural design and modeling of such big projects was based on the classical continuum mechanics, despite the fact that some employed materials was discrete by nature, e.g. the railway ballast which is an assembly of stones. It should be noted that the discrete materials can be considered as continuous materials at a macroscopic scale, i.e. when the particle assembly reaches a certain size, which characterizes a so-called representative elementary volume (REV). A few people turned their attention toward the study of the discrete materials, while the theory of continuum mechanics seems to be stable at that time. In 1944, Casagrande and Carillo [4] seemed to be the pioneers in the study of discrete materials. They modeled strength anisotropy in soils, which can be distinguished into the inherent and the induced anisotropy. In the following decade, Dantu [5] discovered the contact network in granular materials by performing photoelastic experiments on a pile of translucent disks. Such network was disorder, inhomogeneous, and anisotropic. The well-known stress-dilatancy relationship in soil mechanics, which was assumed to depend on the internal geometry constraint under the basis of the minimum energy principle, was derived by Rowe [6] in 1962. The mechanism of deformation implied in Rowe's stress-dilatancy was then proved by Horne [START_REF] Horne | The behaviour of an assembly of rotound, rigid, cohesionless particles. I[END_REF][START_REF] Horne | The behaviour of an assembly of rotound, rigid, cohesionless particles. II[END_REF][START_REF] Horne | The behaviour of an assembly of rotound, rigid, cohesionless particles. III[END_REF]. He provided a fact that the evolution of fabric anisotropy is dependent of the volume, i.e. there is a limiting state of anisotropy which is approached after large deformations have occurred. Weber [START_REF] Weber | Recherches concernant les contraintes intergranulaires dans les milieux pulvérulents[END_REF] introduced the concept of stress tensor in a granular media considered as a continuum material, namely an average stress tensor in a volume V, which is related to the average of the local contact forces in this volume.

Nonetheless, Kakalios [3] noted that, dating back to the 17 th century, the basic properties of granular media were first investigated by the work of Coulomb on slope stability. Coulomb proposed the idea of static friction. In fact, the concept of dilatancy phenomenon was introduced by Reynolds in the late 18 th century. This phenomenon can describe why the wet sand around your feet becomes dry when you are walking on the beach. The deformation induced by the pressure of the foot causes the sand to expand and permit the water to drain out. Both works of these two pioneers indicate a longstanding history of the science of granular materials started with outstanding researchers, who dedicated their time to understand and find the way to deal with the mechanical responses of these materials.

The understanding of granular materials at that moment was mostly obtained from experiments conducting on both analogue and real materials. Some theoretical aspects relied on idealized materials were also established. However, the lack of micromechanical data was encountered at the end of 1960s. Such problem can be partially solved by using computer simulations, which provide a relationship between the macroscopic mechanical behavior and the local variables based on the statistical mechanics. In 1971, Cundall [START_REF] Cundall | A computer model for simulating progressive large scale movements in blocky rock systems[END_REF] was the first pioneer who developed a new approach called distinct element method (DEM) for rock mechanics. This approach is actually suited for any discrete materials, contrary to the finite element method which is limited to continuous materials. Note that the term of "distinct" was replaced later by "discrete" in order to emphasize that the DEM is convenient for the nature of discrete system [12].

The DEM considers all grains of a granular medium as individual rigid bodies interacting between them through contact laws. The grain displacements and rotations are governed by the Newton's equation of motion [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. There are two widely used approaches applied to granular materials, which belong to the DEM: the smooth DEM so-called molecular dynamics (MD) method and the non-smooth DEM referred to as contact dynamics (CD) method. The former method is based on an explicit algorithm which considers the contact interactions as a function of the relative displacements and the relative velocities. In addition, the regularized form of the Coulomb's law of friction is implemented. The latter method was proposed by Moreau and Jean [START_REF] Radjaï | Contact dynamics as a nonsmooth discrete element method[END_REF]. Noncontinuous functions, i.e. Unilaterallity condition and Coulomb's law of friction, can be described within the framework of convex analysis. Hence, an implicit algorithm is required to tackle with the dynamical equation and the non-smooth frictional contact laws. The fundamental difference between both methods [START_REF] Radjaï | Contact dynamics as a nonsmooth discrete element method[END_REF][START_REF] Radjaï | Multicontact dynamics of granular systems[END_REF][START_REF] Radjaï | Contact forces in a granular packing[END_REF] is that MD method requires a fine resolution of the time for the numerical integration of the equations of motion, that is why the MD method is slower than the CD method in some cases. The MD method also seems to be less precise than the CD method in regard to a friction handling. In contrast, the CD method cannot be applied in some cases involving the deformation of particles since this method does not consider the elasticity of contacts.

The numerical simulations have been a powerful and reliable research tool in the domain of granular materials since then until the present.

Due to the dramatic advances in camera technology, microcomputer, and image processing in terms of algorithms and softwares, non-contact measurement techniques so-called full-field measurements become increasingly popular and attractive for engineers and scientists in the experimental mechanics field [START_REF] Grédiac | The use of full-field measurement methods in composite material characterization: interest and limitations[END_REF][START_REF] Viggiani | Full-field measurements, a new tool for laboratory experimental geomechanics[END_REF]. Under a given mechanical loading, such methods can measure the spatial distributions of physical variables at the surface of specimens or even in their bulk in some cases, such as temperatures (scalars), displacements (vectors in two or three dimensions), and strains (second-order tensors) [START_REF] Viggiani | Full-field measurements, a new tool for laboratory experimental geomechanics[END_REF][START_REF] Grédiac | Full-field measurements and identification in solid mechanics[END_REF]. [START_REF] Misra | Measured kinematic fields in the biaxial shear of granular materials[END_REF][START_REF] Hall | Localised deformation patterning in 2D granular materials revealed by digital image correlation[END_REF][START_REF] Rechenmacher | Characterization of mesoscale instabilities in localized granular shear using digital image correlation[END_REF][START_REF] Richefeu | An experimental assessment of displacement fluctuations in a 2D granular material subjected to shear[END_REF], particle image velocimetry [START_REF] Tischer | Kinematics of sand avalanches using particle-image velocimetry[END_REF][START_REF] White | Soil deformation measurement using particle image velocimetry (PIV) and photogrammetry[END_REF][START_REF] Pudasaini | Velocity measurements in dry granular avalanches using particle image velocimetry technique and comparison with theoretical predictions[END_REF][START_REF] Słomiński | Deformation measurements in granular bodies using a particle image velocimetry technique[END_REF][START_REF] Słomiński | Application of particle image velocimetry (PIV) for deformation measurement during granular silo flow[END_REF], digital speckle radiography [START_REF] Grantham | Measurement of granular flow in a silo using digital speckle radiography[END_REF], and stereo-photogrammetry [START_REF] Finno | Strain localization and undrained steady state of sand[END_REF][START_REF] Mokni | Strain localization measurements in undrained plane-strain biaxial tests on Hostun RF sand[END_REF][START_REF] Desrues | Strain localization in sand: an overview of the experimental results obtained in Grenoble using stereophotogrammetry[END_REF]; or in volume measurement techniques such as x-ray tomography [START_REF] Bransby | Deformations near rupture surfaces in flowing sand[END_REF][START_REF] Drescher | Kinematics of the mass flow of granular material through a plane hopper[END_REF][START_REF] Aste | Geometrical structure of disordered sphere packings[END_REF][START_REF] Jerkins | Onset of mechanical stability in random packings of frictional spheres[END_REF][START_REF] Zou | The packing of granular polymer chains[END_REF], electrical capacitance tomography [START_REF] Zhu | Electrical capacitance tomography measurements on vertical and inclined pneumatic conveying of granular solids[END_REF][START_REF] Chaniecki | Application of electrical capacitance tomography for bulk solids flow analysis in silos[END_REF][START_REF] Grudzien | Description of the silo flow and bulk solid pulsation detection using ECT[END_REF][START_REF] Grudzien | ECT image analysis methods for shear zone measurements during silo discharging process[END_REF], magnetic resonance imaging [START_REF] Nakagawa | Non-invasive measurements of granular flows by magnetic resonance imaging[END_REF][START_REF] Hill | Bulk segregation in rotated granular material measured by magnetic resonance imaging[END_REF][START_REF] Nguyễn | Segregation in horizontal rotating cylinders using magnetic resonance imaging[END_REF], positron emission particle tracking [START_REF] Parker | Positron emission particle tracking studies of spherical particle motion in rotating drums[END_REF][START_REF] Stewart | Simulated and measured flow of granules in a bladed mixer -a detailed comparison[END_REF][START_REF] Ingram | Axial and radial dispersion in rolling mode rotating drums[END_REF][START_REF] Laurent | Comparative study by PEPT and DEM for flow and mixing in a ploughshare mixer[END_REF][START_REF] Windows-Yule | Effects of packing density on the segregative behaviors of granular systems[END_REF],

and radar-based sensing [START_REF] Hill | Granular segregration studies for the development of a radar-based three-dimensional sensing system[END_REF]. Sometimes, these methods are mutually applied to extract the kinematic data in the granular material under study [START_REF] Michalowski | Flow of granular material through a plane hopper[END_REF][START_REF] Oda | Microstructure of shear bands and its relation to the mechanisms of dilatancy and failure of dense granular soils[END_REF][START_REF] Mueth | Signatures of granular microstructure in dense shear flows[END_REF][START_REF] Wolf | Experimental investigation of shear band patterns in granular material[END_REF][START_REF] Mcdonald | Particle movement during the deep penetration of a granular material studied by x-ray microtomography[END_REF][START_REF] Hall | Discrete and continuum analysis of localised deformation in sand using xray µCT and volumetric digital image correlation[END_REF][START_REF] Andò | Grain-scale experimental investigation of localised deformation in sand: a discrete particle tracking approach[END_REF]. It is worth noting that contact forces are not easily obtained from the experiments, whereas they can directly evaluate from simulations through various mathematical and modeling hypotheses.

There are two available full-field techniques which can be used to measure stresses inside the particles. First, photoelasticimetry is a powerful method to provide the difference between the principal stresses in each particle made from birefringent materials. This difference is proportional to the shear stress. The method has been employed for a long time [5,[START_REF] Drescher | Photoelastic verification of a mechanical model for the flow of a granular material[END_REF][START_REF] Oda | Experimental micromechanical evaluation of strength of granular materials: effects of particle rolling[END_REF][START_REF] Sadd | Contact law effects on wave propagation in particulate materials using distinct element modeling[END_REF][START_REF] Veje | Kinematics of a twodimensional granular couette experiment at the transition to shearing[END_REF][START_REF] Howell | Stress fluctuations in a 2D granular couette experiment: a continuous transition[END_REF]. It is now common and broad in the literature on granular materials in the past fifteen years [START_REF] Geng | Footprints in sand: the response of a granular material to local perturbations[END_REF][START_REF] Geng | Memory in twodimensional heap experiment[END_REF][START_REF] Hartley | Logarithmic rate dependence of force networks in sheared granular materials[END_REF][START_REF] Geng | Green's function measurements of force transmission in 2D granular materials[END_REF][START_REF] Majmudar | Contact force measurements and stress-induced anisotropy in granular materials[END_REF][START_REF] Majmudar | Jamming transition in granular systems[END_REF][START_REF] Zhang | Statistical properties of a 2D granular material subjected to cyclic shear[END_REF][START_REF] Tordesillas | Structural stability and jamming of self-organized cluster conformations in dense granular materials[END_REF][START_REF] Clark | Particle scale dynamics in granular impact[END_REF][START_REF] Puckett | Equilibrating temperaturelike variables in jammed granular subsystems[END_REF][START_REF] Zhang | A comparison between bridges and force-chains in photoelastic disk packing[END_REF]. The second method is the so-called thermoelastic stress analysis (TSA) [START_REF] Grédiac | Full-field measurements and identification in solid mechanics[END_REF][START_REF] Dulieu-Barton | Development and applications of thermoelastic stress analysis[END_REF][START_REF] Pitarresi | A review of general theory of thermoelastic stress analysis[END_REF][START_REF] Greene | Thermoelastic stress analysis[END_REF]. It is based on infrared (IR) thermography, which captures the temperature fields on the surface of specimens. TSA method provides the sum of the principal stresses, which is proportional to the hydrostatic stress, at the surface of a body. Let us note that there are a few studies based on IR thermography on granular medium: the measurement of mechanical dissipation in sands and soils [START_REF] Luong | Charateristic threshold and infrared vibrothermography of sand[END_REF][START_REF] Luong | Introducing infrared thermography in soil dynamics[END_REF], the characterization of heat transfer in moving spheres [START_REF] Sepúlveda | Infrared tracking from morphological image processing tools -application to heat transfer characterization in granular media[END_REF],

and estimation of soil surface microrelief and rill morphology [START_REF] De Lima | Can infrared thermography be used to estimate soil surface microrelief and rill morphology?[END_REF]. It can be seen that the IR thermography has never been used for reveal the stress network in granular assembly. This inspired the author to apply the TSA method as a new experimental approach to study granular materials. In the author's opinion, the measurement of hydrostatic stress networks using the TSA technique would be a novel approach for the community of granular materials, which is complementary to other full-field measurement techniques and also numerical simulations.

The main motivation of the present work is therefore to reveal the hydrostatic stress network in granular materials by means of TSA method. For that purpose, IR thermographic measurements are performed during mechanical tests at ambient temperature. MD simulation is also used as a complementary tool. In order to investigate the mechanical behavior of granular medium, the experimental and numerical results will be then analyzed through the granular texture of the tested specimens.

Most of former studies on granular materials were focused on non-composite media, i.e. made of only one type of constitutive materials. In this work, two different polymeric materials are utilized to prepare an assembly of two-dimensional analogue cohesionless composite granular materials. In fact, an interesting idea for the future is to design and build composite granular materials for new applications by tuning the proportions of particles made of different constitutive materials. An advantage would be cheap and light-weight materials with other materials. In the present study, twelve configurations differing in terms of diameter size and number of particle are studied both with TSA experiments and MD simulations (chapter 3). Due to some restrictions in the experiments, only simulation is employed for two additional studies with the objective to investigate the influence of particle number ratio and contact stiffness ratio on the granular texture (chapters 4 and 5). Finally, a short study is started on cohesive granular media (chapter 6).

The outline of this dissertation is organized into six chapters followed by conclusion and perspectives.

Chapter 1 is devoted to review the useful theoretical aspects for the present study: granular materials, MD simulation, and TSA technique. The first topic concerns the static properties of granular materials. Many interesting mechanical behaviors are presented through the granular texture. A definition of the stress tensor in granular materials also defined both at continuum scale and grain scale. The second topic concerns the numerical integration scheme so-called "Gear predictor-corrector method" that exists in the literature for solving the Newton's equations of motion. The integration scheme and the contact forces model used in this study are then explained.

The workflow of MD method is also introduced. The final topic concerns the TSA technique. The basis of IR radiation and infrared systems is then described. The IR camera utilized in this study is presented in terms of specifications and calibration methods. The processing related to TSA is also described.

Before performing experiments and simulations, some material parameters need to be identified. Chapter 2 provides three elementary experiments to evaluate the values of contact stiffness, the thermoelastic calibration factor, and the surface energy.

Contact stiffnesses will be used later to perform MD simulations. The thermoelastic calibration factors will be used for converting temperature variations to hydrostatic stresses. The surface energy will be used for the final experiment on cohesive granular meida (in chapter 6).

Chapter 3 concerns the main study in this thesis. It aims at analyzing the mechanical behavior of composite granular materials through the granular texture. TSA experiments and MD simulations are performed for this purpose. The experimental and numerical results are presented in term of normalized hydrostatic stress fields. Note that, the normalized hydrostatic stress in each particle is obtained by the ratio of the hydrostatic stress in that particle to the average value of the hydrostatic stress over all particles in the granular system. A statistical processing of the results is done in term of angular distribution of contact directions and probability distribution of hydrostatic stresses.

In chapters 4 and 5, MD simulations are employed to go further in the analysis.

The fourth chapter handles with the influence of the ratio of particle numbers on the granular texture, while the effect of the ratio of contact stiffnesses on the granular texture is described in the fifth chapter. It should be noticed that the contact forces are here used in the analysis instead of the hydrostatic stresses. The textural parameters used in both chapters consist of coordination number, probability distribution of contact forces, as well as angular distribution of contact directions, except for the fifth chapter which considers one more parameter (the angular distribution of the contact forces). Let us finally underline that the choice of constitutive materials employed in this study (polymeric materials) was done for convenient reasons in the experiments. The application to particles made of rocks is a perspective of this work.

The author expects this thesis to be more or less useful for the reader. However, if this thesis has any error particularly in regard to English language, the author would like to apologize for any inconvenience caused by this error.

CHAPTER 1 THEORETICAL BACKGROUND

Introduction

The purpose of this chapter is to review and to summarize the relevant theories used in this dissertation. This bibliographical chapter is organized as follows:

 Section 1.
2 is a general presentation of the granular materials. Due to a discrete nature, these materials exhibit complex mechanical responses, i.e.

solid-like or fluid-like behaviors, which differ from ordinary substances.

Nonetheless, this section will consider only the solid-like case at the static equilibrium, which is relevant in the framework of this thesis. The interesting static properties of granular material will be discussed using the granular texture, in terms of compactness, coordination number, geometrical fabric, and mechanical fabric. The distributions of contact forces and the definition of stress tensor in granular materials are also discussed in this section;

 Section 1.3 presents the framework of the molecular dynamics (MD) simulations. Concerning the solving of the Newton's equations of motion, Gear predictor-corrector integration scheme based on the finite difference methods is described. In addition, the workflow of MD simulation is given in this section. The contact force models that will be used for both non- and IR thermographic systems are also given in this section. Furthermore, the IR camera used in this study is presented in terms of specifications and calibration procedures. In order to achieve an accurate stress data, the procedure to determine the so-called thermoelastic calibration factor is also described.

Granular Materials

The definition of granular materials is here emphasized again that these material are a collection of solid particles whose macroscopic mechanical behavior is governed by the interaction forces between particles. They are constituted from grains with a variety of shapes, sizes, and material types. An example of granular materials is demonstrated in Figure 1.1. Despite a great diversity in grains, all these media share elementary characteristics: they are chaotic at the grain scale but solid-like or fluid-like behaviors are displayed at the continuum scale. However, the nature of granular materials is generally complex with a broad variety of phenomena that are still partly understood.

Figure 1.1 examples of granular materials which can be found in our everyday life Several difficulties dealing with granular materials are encountered apart from common substances like solids, liquids, and gases. These difficulties can be summarized as follows [START_REF] Jaeger | Physics of the granular state[END_REF][START_REF] Jaeger | Granular solids, liquids, and gases[END_REF][3][START_REF] De Gennes | Granular matter: a tentative view[END_REF][START_REF] Jaeger | Does the granular matter?[END_REF][START_REF] Herrmann | Granular matter[END_REF][START_REF] Jaeger | Sand, jams, and jets[END_REF][START_REF] Andreotti | Granular media: between fluid and solid[END_REF]. First, granular media cannot be easily categorized into solids, liquids, or gases as same as ordinary matters: a sand pile at rest which its angle less than angle of repose exhibits solid-like behavior, while fluid-like behavior is observed when grains at the surface starts to flow (pile angle higher than angle of repose). Second, the contact interactions between individual grains are non-trivial, highly non-linear, and even discontinuous owing to friction and inelastic collision. Energy is also simply dissipated by friction and collision between grains. Third, the thermal agitation in granular medium is negligible compared to external forces such as the gravity, when the particle size greater than 1 m. Brownian motion cannot thus be seen. Next, granular medium are composed of a large number of particles, e.g. sand in a bucket which a number of sand grains in the bucket is approximately more than a million grains. In order to model a huge number of particles in granular media as a continuum, a determination of an average quantity would be a great strategy. However, the use of continuum mechanics to describe the mechanical behavior of granular media seems questionable. Finally, a lack of clear scale separation between microscopic (e.g. particle sizes) and macroscopic scale (e.g. size of flow) is occurred from the difficulty of the continuum description: sand pile generally flows down with 10-20 particle diameters of the flow thickness; likewise the thickness of shear band localized in a granular soil is about a few tens of particle width.

This section is not an overall review of all aspects of granular materials, only granular solids at the static equilibrium condition are considered. Their interestingly major properties are described hereinafter. It is important to note that the granular texture, which means both force and contact networks for multi-contact systems, is also introduced in this section. The granular texture can be described by the scalar parameter such as solid fraction and coordination number, including described the anisotropy of the texture by the second-order tensors i.e. fabric tensors. In this dissertation, the words "textural properties" or "micromechanical analysis" are sometimes mentioned, which both are consistent with the granular texture. Additionally, the characteristics of force transmission in granular materials are considered as one of parameters in the granular texture.

Granular Packing

A heterogeneous structure is generally formed by grains, which a packing of grains becomes an important property to characterize the rigidity and flow of granular materials [START_REF] Jaeger | Physics of the granular state[END_REF]. That is a reason why the packing of granular materials have been studied for a long time. These packing issues are also primary aspects to many industries. The studies in the field of packing can be characterized by some variables, such as volume fraction, porosity, and void fraction. The volume fraction η or sometimes called "solid fraction" is defined as the ratio between the volume filled by the grains inside the representative elementary volume (REV) and the total volume of the REV:

  grains REV V V (1.1)
Note that the volume fraction cannot exceed 1. Considering the porosity ϵ, it is equal to the volume of voids divided by the whole volume, which corresponds to the volume fraction by ϵ = 1 -η. This quantity is often used in the field of porous media. Another is the void fraction, which is usually found in the soil mechanics. This variable is represented by the ratio of the volume of voids to the volume of grains. However, all these variables are equivalent. In this dissertation, only the concept of volume fraction η will be employed.

To get a rid of complexity of non-spherical grains, most of studies in the beginning focus on spheres with the same size. In theory, the famous Kepler conjecture speculates that the densest possible packing of spheres in three-dimensional Euclidean space is the face-centered cubic and the closepacked hexagonal structure, which has the volume fraction η = 0.74 [START_REF] Jaeger | Physics of the granular state[END_REF][START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Berryman | Random close packing of hard spheres and disks[END_REF][START_REF] Donev | Unusually dense crystal packings of ellipsoids[END_REF]. This conjecture has already been mathematically proven by Hales [START_REF] Hales | A proof of the Kepler conjecture[END_REF].

Nevertheless, it is difficult to achieve the theoretical densest packing in practice. When a collection of spherical grains is compacted into a packing, the grains will be typically formed to a stable random arrangement. The final volume fraction for this stable packing cannot be outside the range between two limits so-called "random close packing" (RCP) limit and "random loose packing" (RLP) limit. Generally, the random close packing is referred to as the densest state achieved by randomly packed [START_REF] Dong | Role of interparticle forces in the formation of random loose packing[END_REF]. For uniform spheres, the well-known volume fraction in the RCP limit is about η RCP = 0.64 [START_REF] Berryman | Random close packing of hard spheres and disks[END_REF][START_REF] Zhang | A simulation study of the effects of dynamic variables on the packing of spheres[END_REF][START_REF] Song | A phase diagram for jammed matter[END_REF][START_REF] Baranau | Random-close packing limits for monodisperse and polydisperse hard spheres[END_REF]. Contrary to the RCP, the random loose packing is defined as the loosest possible random packing that is mechanically stable [START_REF] Jaeger | Physics of the granular state[END_REF][START_REF] Dong | Role of interparticle forces in the formation of random loose packing[END_REF][START_REF] Onoda | Random loose packing of uniform spheres and the dilatancy onset[END_REF]. It corresponds to the volume fraction of about η RLP = 0.56 for uniform spheres [START_REF] Dong | Role of interparticle forces in the formation of random loose packing[END_REF][START_REF] Onoda | Random loose packing of uniform spheres and the dilatancy onset[END_REF]. These results are relevant to the dilatancy phenomena, which occurs from spreading of a dense packing to make a space for grains motion. The onset of dilatancy corresponds to the RLP limit without the gravitational effect g  0 [START_REF] Onoda | Random loose packing of uniform spheres and the dilatancy onset[END_REF].

In fact, the packing of particles depends on several factors [START_REF] Zhang | A simulation study of the effects of dynamic variables on the packing of spheres[END_REF][START_REF] Yang | Computer simulation of the packing of fine particles[END_REF] such as preparation method (dropping rate, and dropping height), material properties (friction coefficient and restitution coefficient), particle size distribution (monodisperse to polydisperse), particle shape (spherical to non-spherical grains), and absolute particle size (coarse to fine grains). The interaction forces are also significant for the packing of particles [START_REF] Yu | On the relationship between porosity and interparticle forces[END_REF]. As a consequence, there are many attempts to investigate the influence of such factors on the packing of particles in terms of RCP and RLP [START_REF] Berryman | Random close packing of hard spheres and disks[END_REF][START_REF] Donev | Unusually dense crystal packings of ellipsoids[END_REF][START_REF] Dong | Role of interparticle forces in the formation of random loose packing[END_REF][START_REF] Yang | Computer simulation of the packing of fine particles[END_REF][START_REF] Yu | On the relationship between porosity and interparticle forces[END_REF][START_REF] He | Computer simulation of random packing of unequal particles[END_REF][START_REF] Kansal | Computer generation of dense polydisperse sphere packings[END_REF][START_REF] Williams | Random packings of spheres and spherocylinders simulated by mechanical contraction[END_REF][START_REF] Okubo | Random packing of binary hard discs[END_REF][START_REF] Donev | Improving the density of jammed disordered packings using ellipsoids[END_REF][START_REF] Xu | Packing structure of cohesive spheres[END_REF][START_REF] Man | Experiments on random packings of ellipsoids[END_REF][START_REF] Valverde | Random loose packing of cohesive granular materials[END_REF][START_REF] Lochmann | Statistical analysis of random sphere packings with variable radius distribution[END_REF][START_REF] Al-Raoush | Simulation of random packing of polydisperse particles[END_REF][START_REF] Delaney | Defining random loose packing for nonspherical grains[END_REF][START_REF] Jin | Statistical theory of correlations in random packings of hard particles[END_REF]. Figure 1.2 presents the relationship between normalized shear stress q/p and cumulative shear strain ɛ for two different initial solid fractions [START_REF] Radjaï | Shear strength of granular materials[END_REF]. It is observed that the normalized shear stress increases with the cumulative shear strain. These increments are clearly different for dense system and loose system. Nonetheless, the normalized shear stress for both systems eventually approaches to the characteristic value of shear stress.

Base on the fact that all particles in the granular packing interact at the contact points, another parameter regarding the compactness [START_REF] Radjaï | Contact dynamics study of 2D granular media: Critical states and relevant internal variables[END_REF] could be thinking of. This parameter is a coordination number z [START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF], which is equivalent to an average number of contacts per a particle. z can be calculated by the ratio of the total number of contacts N c to the total number of particles N p as defined in Eq. (1.2). The factor of two in this equation means that there are in fact two particles in one contact. In any granular packing, the coordination number is an important quantity because contacts between particles in the packing need a sufficient number of the mechanical constraints to ensure the stability of packing [START_REF] Silbert | Geometry of frictionless and frictional sphere packings[END_REF][START_REF] Blouwolff | The coordination number of granular cylinders[END_REF]. Considering the d-dimensional static packing consisting of N frictionless particles, there are equilibrium equations of dN for determination of N c contact forces between each particle. In this static packing that are stable against external disturbances, the coordination number is equal to 4 for two-dimensional and 6 for three-dimensional [START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Silbert | Geometry of frictionless and frictional sphere packings[END_REF]. Such packings are called "isostatic" [START_REF] Krapf | Force propagation in isostatic granular packs[END_REF][START_REF] Roux | Geometric origin of mechanical properties of granular materials[END_REF], which implies the statistically determinate systems. In other words, this isostatic packing mean the total number of contacts is equal to the equilibrium equations of dN. It is significant to note that the isostatic packing is highly sensitive to disturbances, because the stability of the packing can be disrupted by a breakage of only one contact [START_REF] Andreotti | Granular media: between fluid and solid[END_REF]. When friction is taken into account, there are 2N c contact forces in 2D and 3N c contact forces in 3D, while the static packing must satisfy 3N equilibrium equations in 2D and 6N equilibrium equations in 3D. Hence, the coordination number of the dry granular packing with frictional particles is in a range of 3 ≤ z ≤ 4 for 2D and 4 ≤ z ≤ 6 for 3D [START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Silbert | Geometry of frictionless and frictional sphere packings[END_REF]. The static packing with frictional particles might be "hyperstatic", which the contact forces are undetermined [START_REF] Andreotti | Granular media: between fluid and solid[END_REF].

Even though frictionless particles are constructed to the isostatic packing regardless of construction history and coefficient of restitution, the packing with frictional particles is dependent of construction history and system parameters [START_REF] Silbert | Geometry of frictionless and frictional sphere packings[END_REF]. Moreover, it could be remarked that the distribution of particle shapes (circle, polygon, polyhedron, etc.), different types of contact (edge-to-edge, face-to-face, etc.), coefficient of friction, and particle stiffness are several parameters that affect to the coordination number of the packing [START_REF] Silbert | Geometry of frictionless and frictional sphere packings[END_REF][START_REF] Roux | Geometric origin of mechanical properties of granular materials[END_REF][START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF][START_REF] Cambou | Micromechanics of Granular Materials[END_REF]. Figure 1.3 displays the influence of particle shape on the coordination number during shear [START_REF] Azéma | Force transmission in a packing of pentagonal particles[END_REF]. The effect of different types of contact as mentioned above is also taken into account in this figure. In fact, the solid fraction and the coordination number are simple zero-order tensors which provide information in term of compactness of the system. These parameters cannot describe the complexity of particles arrangement in the granular packing, but the stress-strain relationship under a monotonic shear is governed by the compactness of the granular system as illustrated in Figure 1.2 [START_REF] Radjaï | Contact dynamics study of 2D granular media: Critical states and relevant internal variables[END_REF][START_REF] Cambou | Micromechanics of Granular Materials[END_REF]. 1.2.2 Forces in the Granular Packing Figure 1.4 shows an inhomogeneous force distribution in a static granular assembly. In this figure, the "chain-like" structures can be observed in some grains. These force chains seem to carry the most of external applied loading [START_REF] Radjaï | Force distributions in dense two-dimensional granular systems[END_REF]. In fact, it is well-known that the contact forces in granular materials are transmitted through the contacts from one particle to another. The contact networks play as the skeleton of static granular matter, which are an important feature to determine its mechanical properties [START_REF] Ostojic | Scale invariance and universality of force networks in static granular matter[END_REF].

Absolutely, the fundamental objective of mechanics of granular materials is to better understand the contact forces and their spatial correlations, especially in the effect of external applied loading at the system boundaries.

In this section, the distributions and correlations of both force and contact networks are thus discussed.

Contact Forces Measurement

In order to measure contact forces in granular materials, there are various already used methods. The first is photoelasticity, which was initially applied in granular materials by Dantu [5] in 1957. This method uses particles made of a photoelastic material, e.g. pyrex, plexiglass, or polymer, that the property of temporary birefringence is exhibited when subjected to the strain [START_REF] Patterson | Optical methods for solid mechanics: A full-field approach[END_REF]. Figure 1.5a displays the stress pattern in twodimensional granular packing provided by the photoelasticity. Nevertheless, it is very complicate task to extract the contact forces at each contact from the stress pattern of the granular packing. Recently, such problem was successfully solved by Majmudar and Behringer [START_REF] Majmudar | Contact force measurements and stress-induced anisotropy in granular materials[END_REF].

The second method is carbon paper method, which can measure the contact forces exerted by the grains at the boundary of packing. In these measurements, a carbon paper is placed between the system boundary and the materials to make a size of spots proportional to the normal forces as shown in Figure 1.5b [START_REF] Mueth | Force distribution in a granular medium[END_REF][START_REF] Liu | Force fluctuations in bead packs[END_REF]. The first experiments were performed by Liu et al. [START_REF] Liu | Force fluctuations in bead packs[END_REF]. It is possible to practically measure the contact forces in 2D and 3D by the carbon paper method, but it can just measure at the system boundaries. Furthermore, the carbon paper method is required a large sufficient external loading to produce the spots on the paper. This measurement technique was later developed by means of force sensors instead of the carbon paper [START_REF] Løvoll | Force measurements on static granular materials[END_REF] as demonstrated in Figure 1.5c. experimental setup of the developed carbon paper method by using force sensors [START_REF] Løvoll | Force measurements on static granular materials[END_REF] The last approach is studied by numerical simulations based on the discrete element method. The numerical simulations are able to measure the contact forces in the bulk of 3D granular packing, and to eliminate the limitation of grain parameters in the experiments. Moreover, the simulations can provide higher measurement resolution for the small contact forces than the experiments. However, it is a great challenge to model the realistic granular materials by the simulations. An overview of the molecular dynamics (MD) simulation used in this dissertation, which belongs to the discrete element method, is explained next in section 1.3.

Probability Distribution of Contact Forces

In order to statistically analyze the amplitude of contact forces that obtained by various methods as mentioned in the previous section, the probability distribution function P(F) of contact forces between grains is used. The simulation results [START_REF] Radjaï | Force distributions in dense two-dimensional granular systems[END_REF] shown that the contact forces can be distinguished into two complementary networks by the mean force: strong and weak networks. The "strong network" is the contact forces greater than the mean force, while the "weak network" is the contact forces lower than the mean force. The shape of distribution for these two complementary networks is concerned. The distribution in the strong network can be wellestimated by an exponential decay, while a power law can be approximately fitted for the distribution in the weak network. These distributions can roughly written [START_REF] Cambou | Micromechanics of Granular Materials[END_REF][START_REF] Radjaï | Force distributions in dense two-dimensional granular systems[END_REF] by

      1 1 1             F F k F F F F PF F F ke (1.3)
where k is the normalization factor given by

1 1 1 1      k (1.4)
In addition, the relationship between the exponents of the two distributions is obtained in the following:

   2 1 2       (1.5)
Other researchers were also discovered the similar results in terms of the contact force distribution as described above by using experiments [START_REF] Majmudar | Contact force measurements and stress-induced anisotropy in granular materials[END_REF][START_REF] Mueth | Force distribution in a granular medium[END_REF][START_REF] Liu | Force fluctuations in bead packs[END_REF][START_REF] Løvoll | Force measurements on static granular materials[END_REF][START_REF] Tsoungui | Use of contact area trace to study the force distributions inside 2D granular systems[END_REF][START_REF] Brockbank | Contact force distribution beneath a three-dimensional granular pile[END_REF][START_REF] Erikson | Force distributions in three-dimensional compressible granular packs[END_REF], numerical simulations [START_REF] Radjaï | Contact forces in a granular packing[END_REF][START_REF] Antony | Evolution of force distribution in three-dimensional granular media[END_REF][START_REF] Luding | Stress distribution in static two-dimensional granular model media in the absence of friction[END_REF][START_REF] Silbert | Statistics of the contact network in frictional and frictionless granular packings[END_REF][START_REF] Bagi | Statistical analysis of contact force components in random granular assemblies[END_REF], and theoretical models [START_REF] Bagi | Statistical analysis of contact force components in random granular assemblies[END_REF][START_REF] Coppersmith | Model for force fluctuations in bead packs[END_REF][START_REF] Socolar | Average stress and force fluctuations in noncohesive granular materials[END_REF]. Some of these studies indicated that the exponential decreasing of the strong forces is a vigorous characteristic of the force distribution in both 2D and 3D granular materials. On the other hand, the force distribution in the range of weak network is sensitive to the preparation method or the internal state of the packing [START_REF] Radjaï | Contact forces in a granular packing[END_REF][START_REF] Majmudar | Contact force measurements and stress-induced anisotropy in granular materials[END_REF][START_REF] Azéma | Force transmission in a packing of pentagonal particles[END_REF][START_REF] Erikson | Force distributions in three-dimensional compressible granular packs[END_REF][START_REF] Antony | Evolution of force distribution in three-dimensional granular media[END_REF][START_REF] Silbert | Statistics of the contact network in frictional and frictionless granular packings[END_REF]. Furthermore, it is interesting to note that the weak network is nearly 60% of the whole set of contacts [START_REF] Radjaï | Contact forces in a granular packing[END_REF][START_REF] Radjaï | Force distributions in dense two-dimensional granular systems[END_REF].

The influence of the numerical methods, anisotropy, and particle shape on the contact forces distribution was investigated by Richefeu et al. [START_REF] Richefeu | Force transmission in dry and wet granular media[END_REF]. The CD simulation was applied to examine the effect of an anisotropic state as presented in Figure 1.7b. The spherical samples were prepared by an isotropic compression (isotropic state) and a triaxial compression (anisotropic state). The effect of the anisotropy can be observed that the exponent of the exponential reduction is almost unchanged, while the peak close to the mean force vanishes.

Figure 1.7c shows the normal forces distribution of a polyhedral packing under an isotropic compaction simulated by the CD method. The same shape of force distribution for the strong and weak network is also observed: the exponential decreasing and the power-law distribution.

Moreover, it seems that the number of very weak forces greatly increases with the particle shape by enhancing the arching effect. As a consequence, the value of exponent  decrease to 0.97 compared to 1.4 for spheres. In this manner, it can be said that the contact forces in the strong network are stronger but less in number.

Figure 1.8 a) Transmission of normal forces in a highly polydisperse system with size span s = 0.96 [START_REF] Voivret | Multiscale force networks in highly polydisperse granular media[END_REF]. The white particles are excluded from the force-bearing network.

The amplitude of normal force is represented by the line thickness. b) Probability density functions of normal forces for various sizes span s of particle diameter [START_REF] Voivret | Multiscale force networks in highly polydisperse granular media[END_REF] The force transmission in highly polydisperse granular media subjected to a simple shear was investigated by the CD simulation [START_REF] Voivret | Multiscale force networks in highly polydisperse granular media[END_REF].

The particle size span, which is defined by s = (d max -d min )/(d max +d min ), was varied from 0.2 to 0.96 in this study. Figure 1.8a demonstrates a snapshot of the force chains for s = 0.96. It can be observed that the larger particles mostly transmit the strong force chains, while an enormous number of small particles do not carry the force chains. As a result, more inhomogeneous force distribution is exhibited when the size span becomes broader. This can be illustrated by probability density functions of normalized contact forces in Figure 1.8b for various values of s. The force distribution is wider when the size span increases. For the weak forces, the exponent α of their distribution increases as a function of s.

Concerning the distribution of contact forces in cohesive granular media, the existence of cohesive interaction between particles does not change the heterogeneous force network as shown in figure 1.9a. In this figure, the strong force chains both in compressive and tensile network are vigorously correlated. Indeed, the tensile force chains are distributed in parallel to the compressive force chains. Moreover, some particles are in the equilibrium by only tensile forces, referred to as "self-stress" structure [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF].

Figure 1.9 a) Contact network of normal forces in a strongly cohesive granular system simulated by molecular dynamics method [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF]. The red lines represent compressive efforts, while the tensile efforts are presented the green lines. These lines are proportional to the amplitude of normal force. b) Semi-logarithmic plot of probability distribution function P of normalized normal forces for three different friction coefficients in the strongly cohesive system [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF] Figure 1.9b presents the histogram in semi-logarithmic scale of probability of normal forces normalized by the mean force in the granular system with the high values of adhesion. It is observed that the force distribution in the strong network for both compressive and tensile forces fall off exponentially [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF][START_REF] Preechawuttipong | Modélisation du comportment mécanique des milieu granulaires cohésifs[END_REF], as found in non-cohesive case [START_REF] Radjaï | Force distributions in dense two-dimensional granular systems[END_REF]. This property was also discovered in the presence of capillary forces in a wet granular media by Richefeu et al. [START_REF] Richefeu | Force transmission in dry and wet granular media[END_REF].

Angular Distribution of Contact Directions and Contact Forces

From the previous section, the force transmission through the interparticle contacts leads to a highly inhomogeneity of the forces in granular materials. In order to clearly characterize the force network and to describe a chance of contact occurrence in each direction of the space, two useful statistical quantities [START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Cambou | Anisotropy in granular materials at different scales[END_REF], which are relevant to the "fabric tensor ϕ", are considered. The first quantity is the "geometrical fabric", which describes the angular distribution of contact directions [START_REF] Andreotti | Granular media: between fluid and solid[END_REF]. The second quantity is the "mechanical fabric", which describes the angular distribution of contact forces [START_REF] Andreotti | Granular media: between fluid and solid[END_REF]. A difference between the geometrical and mechanical fabric is that the magnitude of the contact forces is considered for the mechanical fabric. It is significant to notice that the fabric tensors consider the geometrical structure of the medium, which contains all the directional information in the space of contact orientations [START_REF] Cambou | Micromechanics of Granular Materials[END_REF][START_REF] Peña | Influence of particle shape on sheared dense granular media[END_REF][START_REF] Kruyt | Micromechanical study of fabric evolution in quasi-static deformation of granular materials[END_REF]. This fabric tensor is a symmetrical second-order tensor that is defined by

1      c c c ij i j i j c c N c n n n n N (1.6)
where N c is the total number of contacts in the medium. The i and j components of the unit vector along the normal direction at contact point c are represented by n i c and n j c , respectively. Note that the symbol … c denotes the average value for the whole set of contacts inside the sample.

Kruyt [START_REF] Kruyt | Micromechanical study of fabric evolution in quasi-static deformation of granular materials[END_REF] mentioned that the fabric tensor can alter owing to the three micromechanical mechanism: contact disruption, contact creation, and contact reorientation. The contact disruption is occurred when contacting particles at an exist contact move away from one another, while the contact creation is a new contact created by non-contacting particles move to touch one another. For the contact reorientation, it means that its orientation changes due to the particle movement that remain in contact.

In case of two-dimensional data, the angular distribution of contact directions θ in the normal direction n can be simply given by the function P(θ) defined that P(θ)dθ provides the number of contacts in normal direction within an angular interval between θ -dθ/2 and θ + dθ/2 [START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular materials[END_REF].

The integration of P(θ) over the whole set of contact orientations is equal to the coordination number z [START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Kruyt | Micromechanical study of fabric evolution in quasi-static deformation of granular materials[END_REF]. Figure 1.10 demonstrates how to determine the angular direction of contacts in 2D granular media. When two particles are in contact, the angular direction of contact between them can be measure by using the angle of normal vector connected between their centroids with respect to the horizontal. The first sample is a granular bed prepared by particle free falls from a line source, i.e. a uniform "rain" of grains [START_REF] Radjaï | Micro-statistical features of cohesionless granular media[END_REF].

The effect of gravity leads these particles go down and reach a compaction.

This sample clearly exhibits two peaks along an approximate direction of 45 with respect to the vertical. This feature is the inherent anisotropy of granular medium, which was also reported in experiments [START_REF] Oda | Some experimentally based fundamental results on the mechanical behavior of granular materials[END_REF][START_REF] Calvetti | Experimental micromechanical analysis of a 2D granular material: relation between structure evolution and loading path[END_REF]. The second sample is created by an isotropic compression under the condition of the stress ratio ( 1 / 3 ) equals to 1 [START_REF] Liu | Microscopic study on stress-strain relation of granular materials[END_REF]. An isotropic angular distribution of contacts is observed in this sample. The last sample is compacted by a biaxial compression under the condition of a constant pressure, i.e.  1 increases with the same value as  3 decreases [START_REF] Liu | Microscopic study on stress-strain relation of granular materials[END_REF]. In this case, the contacts are preferentially established in the direction of  1 which is the highest applied loading in the system. It is evident that the preparation method and the external applied loading powerfully affect the geometrical fabric of granular materials [START_REF] Cambou | Anisotropy in granular materials at different scales[END_REF][START_REF] Calvetti | Experimental micromechanical analysis of a 2D granular material: relation between structure evolution and loading path[END_REF][START_REF] Nouguier-Lehon | Influence of particle shape and angularity on the behaviour of granular materials: a numerical analysis[END_REF][START_REF] Li | Influence of loading direction on the behavior of anisotropic granular materials[END_REF][START_REF] Li | Numerical investigation of granular material behaviour under rotational shear[END_REF][START_REF] Claudin | Static properties of granular materials[END_REF][START_REF] Kanatani | Distribution of directional data and fabric tensors[END_REF]. In general, the angular distribution P(θ) of contact directions can be estimated by using the second-order Fourier expansion [START_REF] Cambou | Micromechanics of Granular Materials[END_REF][START_REF] Cambou | Anisotropy in granular materials at different scales[END_REF][START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular materials[END_REF][START_REF] Calvetti | Experimental micromechanical analysis of a 2D granular material: relation between structure evolution and loading path[END_REF][START_REF] Claudin | Static properties of granular materials[END_REF][START_REF] Kanatani | Distribution of directional data and fabric tensors[END_REF][START_REF] Thornton | Computer simulated deformation of compact granular assemblies[END_REF][START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]:

    1 1 cos 2 2            c c P a (1.7)
where a c is the second-order geometrical anisotropy in the contact orientations θ and θ c is the orientation of the principal direction. These parameters can be calculated by using the definition of the fabric tensor as defined in Eq. (1.6). By this definition of the fabric tensor, tr(ϕ) = 1. The geometrical anisotropy of the contact network is given by the difference between the principal value ϕ 1 and ϕ 2 of fabric tensor [START_REF] Cambou | Micromechanics of Granular Materials[END_REF][START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF], which can be written by

  1 2 2     c a (1.8)
It is interesting to distinguish between the strong and weak networks as analyzed regarding the force distribution P(F) in the previous section. In this manner, the angular distribution of contact forces greater than the mean force and the angular distribution of contact forces less than the mean force are separately plotted in the same polar diagrams. Figure 1.12 presents the amplitude of mechanical isotropy and the angular distribution of contact forces for a 2D system of circular particles subjected to a biaxial compression [START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF].

Figure 1.12 a) Amplitude of mechanical anisotropy as a function of normalized contact forces [START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]. b) Angular distributions of contact forces for 2D granular packing obtained from numerical simulations [START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF] It can be observed from Figure 1.12a that the geometrical anisotropy for the two complementary networks is dissimilar: the contact network corresponding to strong forces is much more anisotropic than the contact network corresponding to weak forces. In addition, the geometrical anisotropy for the weak network is distributed in the opposite direction to the axis of compression. Note that a positive value of a c relates to the major principal stress direction, whereas a negative value corresponds to the orthogonal direction [START_REF] Cambou | Micromechanics of Granular Materials[END_REF]. As demonstrated in Figure 1.12b, the contact orientations of strong forces are preferentially oriented along the axis of compression. It should be note that the contacts in the strong network are less than 40% of all contacts in the granular system [START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF].

Although a stress tensor in granular materials will be defined in the next section, the contribution of the strong and weak networks to the stress tensor is discussed as follows. The strong network totally carries the whole deviatoric stress that is shear stress in the medium, while the hydrostatic stress is supported by the weak network [START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Cambou | Micromechanics of Granular Materials[END_REF][START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]. This means the granular materials withstand shear and exhibit solid-like behavior due to the strong forces network. However, the strong network could not be stable in the absence of the weak lateral forces, which are required to the buckling of the strong force chains.

from Forces to Stresses

Granular materials are in fact discrete materials, which the Cauchy stress tensor in the meaning of continuous materials cannot be defined.

Moreover, the study of forces distribution as described in the previous section was demonstrated that the forces between contacting particles are greatly heterogeneous. As a result, some questions regarding the definition of stress tensor in granular materials are arisen. In this section, the stress tensor at continuum and grain scales will be defined from the inter-particles forces in granular medium.

Stress Tensor at Continuum Scale

Many researchers attempt to develop an averaging method with different theoretical considerations to relate the contact forces in discrete materials to the stress tensor in continuum materials [START_REF] Drescher | Photoelastic verification of a mechanical model for the flow of a granular material[END_REF][START_REF] Christoffersen | A micromechanical description of granular material behavior[END_REF][START_REF] Rothenburg | A micromechanical definition of the Cauchy stress tensor for particulate media[END_REF][START_REF] Kanatani | A theory of contact force distribution in granular materials[END_REF][START_REF] Chang | Constitutive relation for a particulate medium with the effect of particle rotation[END_REF][START_REF] Bagi | Stress and strain in granular assemblies[END_REF].

However, the stress tensors proposed by these researchers are in fact equivalent to each other under the equilibrium conditions. In this case, the stress tensor is defined as the average stress of an assembly of particles with arbitrary shape. Here, the definition of the stress tensor at the continuum scale is in brief described by following a grain-by-grain approach [START_REF] Cambou | Micromechanics of Granular Materials[END_REF]. A concept of this approach lies in the determination of the internal moment, i.e. the tensor moment of the internal efforts of a bounded mechanical system.

For the calculation of the stress tensor, the tensorial internal moment M i of each particle p is taken into account [START_REF] Cambou | Micromechanics of Granular Materials[END_REF]. It can be defined by

   p c c ij i j c p M f r (1.9)
where c i f is the i th component of the force acting on the particle at the contact point c, c j r is the j th component of the radius vector of the particle at the same contact c, and the summation in this equation considers over all contact c of connecting particles with the particle i. The component (i, j) of the stress tensor for a granular packing of volume V can be simply written as [START_REF] Cambou | Micromechanics of Granular Materials[END_REF]:

1 1        p c c ij ij i j p V c V M f l V V (1.10)
where c j l is the j th component of the branch vector connecting the centroids of two particles at contact at point c. It is worthy to note that the first summation is performed over all particles whereas the second summation runs over the whole set of contacts in the volume V, with each contact counting once.

Stress Tensor at Grain Scale

In order to calculate the hydrostatic stress in each particle, a definition of the stress tensor at the grain scale is needed. The definition in Eq. (1.10) of the stress tensor at the continuum scale is developed at the grain scale to provide the stress tensor  p in an individual particle:

1 p p c c ij i j cV p fr V     (1.11)
where V p is the volume of the particle. The summation in Eq. (1.11) is done over the whole set of contacts along the particle perimeter.

Molecular Dynamics (MD) Simulation

Due to the dramatic development in microcomputer, computer simulations have been applied for over thirty years to investigate the behavior of granular materials. Why computer simulations are required for granular materials could be discussed by several reasons as follows. The computer simulation can suddenly provide the details of what is actually happening inside the granular system at every time considered, e.g. position, displacement, contact forces. These enable us to investigate and analyze the behavior of granular media both in microscopic and macroscopic quantities. Moreover, it is possible to measure some quantities by using the computer simulations, despite impractical or difficult to determining in the experiments. Significantly, the simulations are able to repeatedly test with exactly the same specimens under different boundary conditions. Nonetheless, the main difficultly of the computer simulations is to develop these approaches achieving adequately accurate results, which compared to the real experiments and the analytical theories.

Granular materials are in fact not continuous materials but they are discrete materials. Hence, the simulation method for continuous materials likes the "Finite Element" (FE) cannot be used in granular medium. In 1971, the numerical simulation so-called "Distinct or Discrete Element Method" (DEM) was originally proposed for rock mechanics by Cundall [START_REF] Cundall | A computer model for simulating progressive large scale movements in blocky rock systems[END_REF]. In fact, the molecular dynamics simulation has been a widely used numerical simulation in the granular materials community. This method is a powerful approach for studying the mechanical behavior of granular media. The MD simulation, which belongs to the DEM, relies on an explicit algorithm. The particles are considered as rigid bodies with non-conforming surfaces [START_REF] Johnson | Contact mechanics[END_REF]. In order to integrate the Newton's equations of motion, there are several used integration schemes under the framework of finite difference method [START_REF] Leach | Molecular modeling: Principles and applications[END_REF][START_REF] Haile | Molecular dynamics simulation: Elementary methods[END_REF][START_REF] Frenkel | Understanding molecular simulation: From algorithms to applications[END_REF][START_REF] Allen | Computer simulation of liquids[END_REF]. However, the MD method applied in this dissertation is based on a predictor-corrector scheme with Gear's set of corrector coefficients [START_REF] Leach | Molecular modeling: Principles and applications[END_REF][START_REF] Haile | Molecular dynamics simulation: Elementary methods[END_REF][START_REF] Frenkel | Understanding molecular simulation: From algorithms to applications[END_REF][START_REF] Allen | Computer simulation of liquids[END_REF].

The so-called Gear predictor-corrector method for integration of the equations of motion, which is based on the finite difference methods, is discussed in the first section. The sketch of MD algorithm is also explained in the next section. The final section provides the model of normal and tangential forces for both cohesive and non-cohesive granular mediums, which are used in this dissertation.

Integration Scheme

In molecular dynamics simulations, a large collection of particles are modeled as a multi-body system interacting through either short-or longrange forces [START_REF] Frenkel | Understanding molecular simulation: From algorithms to applications[END_REF]. The translational and rotational motions of each particle with respect to time are computed by integrating Newton's equation of motion until reaching the equilibrium system [12,[START_REF] Leach | Molecular modeling: Principles and applications[END_REF][START_REF] Frenkel | Understanding molecular simulation: From algorithms to applications[END_REF]. It can be said that the MD method is simple because only the Newton's equation of motion is applied to each particle [12]:

2 2  i i i dx m F dt , i = 1, 2, 3,…, N (1.12)
where N is the number of particles in the granular system. This equation is an ordinary differential equation, which describes the motion of a particle of mass m i along a coordinate x i with the force F i exerted on the particle in that direction. Eq. (1.12) can be easily solved by a simple integration with an initial condition. However, this equation becomes more complicate when considering in realistic models of many interactions between particles.

Under such circumstances, a numerical approach based on a "finite difference methods" are tackled this problem by integrating the equations of motion.

Finite Difference Methods

In molecular dynamics simulations, the finite difference methods are used to generate the motion of each particle in the granular system. The basic idea of these methods can be explained hereinafter [START_REF] Leach | Molecular modeling: Principles and applications[END_REF][START_REF] Allen | Computer simulation of liquids[END_REF]. The integration is broken down into many small stages with a fixed time δt. The positions, velocities, and other kinematics information at time t are given.

The forces of each particle in the configuration at time t are determined. The accelerations of the particles can be then calculated from the forces. All kinematics quantities obtained at time t provide the positions, velocities, etc. at time t + δt. The new kinematics quantities at time t + 2δt, t + 3δt, and so on are step by step solved until achieving a convergent solution (the equilibrium system).

In order to integration of the Newton's equations of motion, there are many different algorithms using the finite difference methods, which have been commonly used in the MD simulations. All algorithms suppose that the positions and other kinematics quantities (velocities, accelerations, etc.) at time t + δt can be estimated as Taylor series expansions regarding time t [START_REF] Leach | Molecular modeling: Principles and applications[END_REF][START_REF] Haile | Molecular dynamics simulation: Elementary methods[END_REF][START_REF] Frenkel | Understanding molecular simulation: From algorithms to applications[END_REF][START_REF] Allen | Computer simulation of liquids[END_REF]:

                            2 3 2 1 1 2 6 1 2                             r t t

r t tv t t a t t b t v t t v t ta t t b t a t t a t tb t b t t b t

(1.13

)
where r is the positions, v denotes the velocities (the first derivative of the positions with respect to time), a indicates the accelerations (the secondorder time derivative), b stands for the third-order time derivative, etc.

Gear Predictor-Corrector Method

For molecular dynamics simulations, there is a variety of numerical integration schemes to solve the equations of motion in the literature, i.e.

Verlet algorithm [START_REF] Verlet | Computer "experiment" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules[END_REF], leap-frog algorithm [START_REF] Hockney | The potential calculation and some applications[END_REF], velocity Verlet method [START_REF] Swope | A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters[END_REF], Beeman algorithm [START_REF] Beeman | Some multistep methods for use in molecular dynamics calculations[END_REF]. Nonetheless, a higher-order algorithm could be appropriate to applied in the MD simulations [START_REF] Frenkel | Understanding molecular simulation: From algorithms to applications[END_REF], i.e. higherorder time derivatives of the particle positions. Such algorithm enables to use a longer time step without minimizing its precision. On the other hand, it is possible to get a higher accuracy for a given time step. The commonly used higher-order algorithm in the MD simulations is "Gear predictorcorrector" method [START_REF] Gear | Numerical initial value problems in ordinary differential equations[END_REF]. This is an efficiency algorithm, which has proven to be powerful for molecular dynamics simulations [START_REF] Pöschel | Computational granular dynamics: Models and algorithms[END_REF]. Importantly, a major advantage of the Gear algorithm is that only one calculation of the interaction forces is required in each time step. In fact, the most consuming computational time is the evaluation of the contact forces between particles in each time step, not integration of the equations [12,[START_REF] Leach | Molecular modeling: Principles and applications[END_REF][START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Pöschel | Computational granular dynamics: Models and algorithms[END_REF]. That is why the Gear predictor-corrector method is more efficiency than the other integration schemes.

The Gear predictor-corrector method is composed of three basic steps:

predictor, evaluation, and corrector. In the predictor step, the positions, velocities, accelerations, and other higher-order time derivatives at time t + δt are predicted by the Taylor expansions as shown in Eq. (1.13). Next, the interaction forces between particles are evaluated from the predicted positions at time t + δt, using various contact force models. The evaluated accelerations a e (t + δt) are obtain by using the Newton's equations of motion. The size of error in the predictor step can be estimated by subtracting between the evaluated and the predicted accelerations:

               e a t t a t t a t t (1.14)
In the corrector step, the size of error is used to correct all the predicted quantities by adding an amount which is proportional to the size of error [START_REF] Haile | Molecular dynamics simulation: Elementary methods[END_REF][START_REF] Pöschel | Computational granular dynamics: Models and algorithms[END_REF]. This correction is differently weighted depending on a degree of the order of time derivative:

                        0 1 2 3                                     c c c c r t t

r t t c a t t v t t v t t c a t t a t t a t t c a t t b t t b t t c a t t

(1.15)

where c 0 , c 1 , c 2 , c 3 … are the set of corrector coefficient suggested by Gear.

These coefficients depend upon the order of the differential equations and the order of the Taylor series expansion [START_REF] Leach | Molecular modeling: Principles and applications[END_REF][START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Pöschel | Computational granular dynamics: Models and algorithms[END_REF]. The corrector step behaves like a feedback mechanism that can diminish an instability obtained from the predictor step [START_REF] Haile | Molecular dynamics simulation: Elementary methods[END_REF]. Consequently, the Gear algorithm provides a good stability. It is very important to note that a fifth-order Gear predictorcorrector algorithm is used for solving the equations of motion in this dissertation.

Finally, the computational time is updated from time t to time t + δt.

The molecular dynamics simulation is then started again at the predictor step by using all the kinematics values from the previous time step. The implementation of the MD simulations will be described in the next section.

Sketch of MD Algorithm

The molecular dynamics simulations are implemented by a general scheme based on the Gear predictor-corrector algorithm, which could be summarized as shown by a flowchart in Figure 1.13. At first, the molecular dynamics simulations are initiated by getting the positions, velocities, and other higher-order time derivatives of each particle at an initial time t. These initial conditions also includes the material properties and the system parameters, e.g. contact stiffness, friction coefficient, number of iterations, position of boundaries, etc.

The positions, velocities, and other higher-order time derivatives at time t + δt are then estimated in the predictor step by using the kinematics quantities at the current time t. This step is the first part of the predictorcorrector integration scheme. Next, the predicted positions at time t + δt are employed to detect the contacts between each particle underlying on a criterion, i.e. the contact between two particles is discovered when the sum of their radius is greater than the distance between their centroids. If there is no contact in a particle, all its kinematics quantities at the next time step are predicted again by using the values at time t + δt of this particle. On the contrary, the MD algorithm continues to the step of force computation when at least one contact is discovered in the particle. For each particle, all the kinematics quantities from the predictor step (positions, velocities, accelerations, etc.) are corrected by using the size of error. It is interesting to note that the desired data for processing later can be stored at either each time step or only the final time. The simulation is then performed at the next time step and so on, which these kinematics quantities of the predictor step at the present time are used to predict the positions, velocities, and other higher-order time derivatives at the next time step. The program is terminated when a number of time step reach a given final value.

On the other hand, it is possible to exit the program by the users.

Contact Forces Calculation

From the previous section, the molecular dynamics algorithm was already presented. One of most important step in this algorithm is the force calculation. Absolutely, there are various contact force models, which are commonly used in the MD simulations of granular materials [START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF]. In this section, only used contact force models in this dissertation are explained.

In the MD method, the force F i as displayed in Eq. (1.12) can be arisen from various effects. In general, this force is composed of an interaction between the particles as well as the external forces. Thus,

,     i ij ext i ji F F F (1.16)
where F ij is the force applied on particle i by particle j. and F ext,i denotes the external force.

Let us consider Figure 1.14 which the particle i and the particle j are in contact. The contact force between these particles can be divided into two components:

  n t ij n t F F F (1.17)
where F n and F t and the normal and tangential components, while n and t stand for unit vectors in the normal and tangential directions. Generally, the normal force causes a change of translational motion of particles, while the tangential force causes the change of rotational and sliding motions of particles. Nevertheless, the rotational motion is permitted to be free in this study. The force models used for calculating F n and F t in the MD simulations of both cohesive and non-cohesive granular materials are then discussed.

Figure 1.14 Geometry of contact between particle i with radius r i and particle j with radius r j Normal Force

In mechanics of granular materials, the interaction force can be generally modeled under the determination of a behavior law and a failure criterion [12]. In the case of a frictionless contact, the normal interaction force can be decomposed into three components. The first component is a repulsive force F n e , which plays a role to resist the deformation or the particle overlap due to an external compression. The second component is relevant to an intrinsic characteristic of granular materials named "energy dissipation" [12]. It is taken into account through a viscous force F n d . The third component is a cohesive force F n c , which inhabits normal separation, shear, and rolling [12,[START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF]. Even though the cohesive interaction between particles may stem from different factors [12], only this cohesive interaction arising from a surface energy at the particle contacts is considered in this dissertation. As a consequence, the normal force can be written by

   e d c n n n n F F F F (1.18)
It is worthy to note that only the first two components of the righthand side in Eq. (1.18) is taken into account in the case of non-cohesive granular materials. For the case of cohesive granular materials, all the components are considered.

In fact, the local deformation occurs when two solid bodies are in contact. The occurrence of this deformation is produced by the interaction force. Since the particles in the MD method are considered as rigid bodies, this deformation is simply modeled by using a virtual overlap δ relying on an assumption that the particle shapes do not change. Hence, the interaction force can be calculated as a function of the virtual overlap. In Figure 1.14, this virtual overlap is given by

       i j i j x x r r (1.19)
where x i and x j are the vector position of the particle centers. This equation implies that no interaction force occurs when the two particles are not in contact (δ > 0). In contrast, the repulsive force is arisen when the particles overlap (δ < 0). Thus, this repulsive force can be simply modeled as a spring [12,[START_REF] Pöschel | Computational granular dynamics: Models and algorithms[END_REF][START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF][START_REF] Tsuji | Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe[END_REF][START_REF] Herrmann | Modeling granular media on the computer[END_REF][START_REF] Schöllmann | Simulation of a two-dimensional shear cell[END_REF]. In this manner, the force is linearly proportional to the virtual overlap:

  e n eff F k (1.20)
where k eff is the effective contact stiffness of particles. If the small virtual overlap between particles is desired, it means that the sufficiently high k eff must be selected. The effective contact stiffness when different types of particles in contact is calculated by (k i +k j )/(k i ×k j ). This is useful for the study of composite granular materials in Chapter 3. In the case of noncomposite granular materials, the effective contact stiffness is evaluated by the contact stiffness of such particles.

Considering the second component in Eq. (1.18), the viscous force is commonly modeled by a viscous damping that withstands the relative velocity for the duration of the contact [START_REF] Haile | Molecular dynamics simulation: Elementary methods[END_REF][START_REF] Pöschel | Computational granular dynamics: Models and algorithms[END_REF][START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF][START_REF] Tsuji | Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe[END_REF][START_REF] Herrmann | Modeling granular media on the computer[END_REF][START_REF] Schöllmann | Simulation of a two-dimensional shear cell[END_REF]:

  d n n n F v (1.21)
where  n is the normal viscous damping coefficient and v n is the normal relative velocity that is related to a time derivative of the virtual overlap.

The relative velocity can be obtained by v i -v j , where v i = dx i /dt is the velocity of particle i and v j = dx j /dt is the velocity of particle j. Thus, the normal relative velocity v n is the inner product between the relative velocity (v i -v j ) and the normal unit vector n.

The third component in Eq. (1.18) corresponding to the cohesive force is considered. There are two widely used models of elastic contact with adhesion. First, the JKR (Johnson-Kendall-Roberts) was developed by

Johnson et al. [START_REF] Johnson | Surface energy and the contact of elastic solids[END_REF]. Second, the DMT (Derjaguin-Müller-Toporov) model was proposed by Derjaduin et al. [START_REF] Derjaguin | Effect of contact deformations on the adhesion of particles[END_REF]. Both models are based on the Hertz elastic contact theory. In addition, the adhesive forces for both models are proportional to the surface energy  and the effective radius r eff (or the contact radius 2a). In this study, the cohesive force with a similar idea to the JKR model, which was developed by Radjaï et al. [START_REF] Radjaï | Cohesive granular texture[END_REF], are applied. This cohesive force is proportional to the contact radius between particles, which the relationship between the contact radius and the virtual overlap in the case of small strain can be provided by a 2 = r eff δ [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF][START_REF] Roux | Quasi-static contacts[END_REF]. In other words, the virtual overlap mainly affects to the cohesive force, which can be written by

   c n eff F r (1.22)
In this dissertation, the "linear spring-dashpot" or so-called "linear viscoelastic" models [START_REF] Haile | Molecular dynamics simulation: Elementary methods[END_REF][START_REF] Pöschel | Computational granular dynamics: Models and algorithms[END_REF][START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF][START_REF] Tsuji | Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe[END_REF][START_REF] Herrmann | Modeling granular media on the computer[END_REF][START_REF] Schöllmann | Simulation of a two-dimensional shear cell[END_REF] is employed to calculate the normal force between two particles in the case of non-cohesive granular media.

From Eq. (1.18), the expression of the normal force for non-cohesive granular materials is thus written by

    n eff n n F k v (1.23)
On the other hand, the "Geometrical Adhesion (GA)" model is used for calculation of the normal force in the case of cohesive granular materials. All the components in Eq. (1.18) are considered, which can be define as

        n eff n n eff F k v r (1.24)
Let us consider the GA model in the case of static or quasi-static as shown in Figure 1.15. The threshold tensile force F c appears as a characteristic force that is necessary for separation of two contacting particles [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF]. It can be said that the threshold tensile force is a maximum tensile force which the contact between particles still exists. In addition, this threshold tensile force occurs at the critical virtual overlap

2 2 4    c eff r k . The threshold tensile force is calculated by 2 4    eff c eff r F k (1.25)
The occurrence of an unstable zone in the range of δ  [0, δ c ] as illustrated in Figure 1.15 depends on the external applied factors. Namely, this overlap range will not occur if the contact between particles is separated by a force-controlled loading. On the contrary, the virtual overlap range between 0 and δ c can be arisen if the inter-particle contact is split by a displacement-controlled loading. It is interesting to remark that the unstable zone does not exist for the contact creation by non-contacting particles move to touch one another. the GA model in the static or quasi-static cases [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF] In this figure, the tensile and compressive zones can be distinguished at an equilibrium position which is represented by an equilibrium virtual overlap δ eq [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF]. There are two possibilities to achieve the equilibrium position: no external applied loading and no interaction forces occur. The virtual overlap at the equilibrium position can be estimated by 22   eq eff rk .

The compressive zone exists for δ > δ eq , while for δ < δ eq enables the tensile zone takes place. Note that more details in the GA model are given in ref. [START_REF] Preechawuttipong | Modélisation du comportment mécanique des milieu granulaires cohésifs[END_REF].

Tangential Force

In general, the tangential force corresponds to the normal force by the Coulomb's law of friction:

  t n F F (1.26)
where μ is the coefficient of friction. This equation implies that it is necessary to restrict the tangential force less than or equal to the multiplication between the friction coefficient and the normal force.

In fact, the Coulomb's law of friction is non-smooth (single-valued) function as demonstrated in Figure 1.16a. This non-smooth relationship implied that there are many possible value of the tangential force before sliding (v s = 0). Furthermore, the integration of the Newton's equations of motion requires a smooth force law in such a way that the tangential force is a mathematical function of the sliding relative velocity. In practice, it can be seen that it is not possible to directly apply the Coulomb's law in the framework of the molecular dynamics simulation. Therefore, the Coulomb's law of friction in the MD method can simply implemented by adopting a "regularized" form of the exact Coulomb's law by adding the term of viscous friction as shown in Figure 1.16b [START_REF] Cambou | Micromechanics of Granular Materials[END_REF][START_REF] Pöschel | Computational granular dynamics: Models and algorithms[END_REF][START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF][START_REF] Herrmann | Modeling granular media on the computer[END_REF][START_REF] Schöllmann | Simulation of a two-dimensional shear cell[END_REF]. The simplest widely used form of the regularized Coulomb's law can be written by

    min , sign     t s s n s F v F v (1.27)
where  s denotes tangential viscosity coefficient, and the sliding relative velocity v s is calculated by the inner product between the relative velocity (v i -v j ) and the tangential unit vector t. It is significant to note that Eq.

(1.27) is used for calculating the tangential force of non-cohesive granular materials in this study. In the case of cohesive granular media, the effect of adhesion should be taken into account. Thus, the tangential force can be computed by the following expression [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF][START_REF] Radjaï | Cohesive granular texture[END_REF]:

      min , sign      t s s n c s F v F F v (1.28)
At this point, the interesting static properties in granular materials were already reviewed including the granular texture for their analysis. The molecular dynamics algorithm, which is used as a numerical tool in this study, was also briefly explained. Moreover, the contact force calculations for non-cohesive and cohesive granular media were given. Thus, the basic principle of "Thermoelastic Stress Analysis (TSA)", which is a novel nondestructive experimental technique for the community of granular materials, will be described in the next section.

Thermoelastic Stress Analysis (TSA) Technique

TSA is a full-field, non-contact, and non-destructive experimental technique providing a field of the first invariant of the stress tensor over the surface of a material subjected to mechanical loading [START_REF] Dulieu-Barton | Development and applications of thermoelastic stress analysis[END_REF][START_REF] Pitarresi | A review of general theory of thermoelastic stress analysis[END_REF][START_REF] Greene | Thermoelastic stress analysis[END_REF][START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF]. This technique is based on infrared (IR) thermographic measurement. Indeed, the elastic strain of a body leads to small temperature changes. This is due to a variation of internal energy, i.e. of the first invariant of the strain tensor, when a material is subjected to a mechanical loading. This phenomenon is sometimes called thermoelastic coupling or isentropic coupling. Under adiabatic evolution together with an assumption of isotropic behavior, the temperature change is proportional to the change of the sum of principal stresses (which is proportional to the hydrostatic stress). This technique has been applied to various mechanical characterizations of materials and structures such as composites [START_REF] Stanley | The application of thermoelastic stress analysis techniques to composite materials[END_REF][START_REF] Emery | A generalized approach to the calibration of orthotropic materials for thermoelastic stress analysis[END_REF][START_REF] Emery | Thermoelastic stress analysis of damage mechanisms in composite materials[END_REF],

ceramics [START_REF] Balandraud | Experimental evidence of thermal effects in multiphase ceramic specimens subjected to cyclic loading[END_REF], aluminiums [START_REF] Pastor | Applying infrared thermography to study the heating of 2024-T3 aluminium specimens under fatigue loading[END_REF][START_REF] Badulescu | Applying the grid method and infrared thermography to investigate plastic deformation in aluminium multicrystal[END_REF], shape memory alloys [START_REF] Delpueyo | Applying infrared thermography to analyse martensitic microstructures in a Cu-Al-Be shapememory alloy subjected to a cyclic loading[END_REF], and sandwich structures or joints [START_REF] Boyd | Stress analysis of finger joints in pultruted GRP materials[END_REF][START_REF] Moutrille | Applying thermoelasticity to study stress relief and crack propagation in aluminium specimens patched with composite material[END_REF]. In case of granular materials, IR thermography was used for the measurement of mechanical dissipation in sands and soils [START_REF] Luong | Charateristic threshold and infrared vibrothermography of sand[END_REF][START_REF] Luong | Introducing infrared thermography in soil dynamics[END_REF] and for the characterization of heat transfer of moving spheres [START_REF] Sepúlveda | Infrared tracking from morphological image processing tools -application to heat transfer characterization in granular media[END_REF]. Nevertheless, the stress network in granular materials has never identified by TSA technique.

This section is organized as follows:

 Section 1.4.1 presents the theory of TSA for isotropic materials from thermodynamical bases.

 Section 1.4.2 gives the fundamentals of IR thermography and application to TSA: IR radiation, IR thermographic systems, and derivation of the thermoelastic relationship which is at the base of TSA.



The specifications of the IR thermographic system which is used in the present study is given in section 1.4.3, including topics related to data acquisition and post-processing.



Finally, section 1.4.4 explains how to experimentally determine the socalled thermoelastic calibration factor (A) in order to obtain reliable quantitative stress data.

Classical Theory of Thermoelastic Stress Analysis

Thermoelasticity is a general term in mechanics to describe the stress changes under elastic limit that occur from the reversible conversion between mechanical deformation and thermal energy in solid continuum bodies. This relationship is well-known as the thermoelastic effect [START_REF] Pitarresi | A review of general theory of thermoelastic stress analysis[END_REF][START_REF] Greene | Thermoelastic stress analysis[END_REF][START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF][START_REF] Stanley | Quantitative stress analysis by means of the thermoelastic effect[END_REF][START_REF] Wong | Thermoelastic constant or thermoelastic parameter?[END_REF]. In fact, it is worthy to note that the heat source induced by strain can be distinguished into two sources [START_REF] Grédiac | Full-field measurements and identification in solid mechanics[END_REF][START_REF] Morabito | Analysis of heat sources accompanying the fatigue of 2024 T3 aluminium alloys[END_REF]. The first source is associated with the energy dissipation, which is arisen from an irreversible mechanical deformation process. The second source corresponds to the thermomechanical couplings. They occur from an interaction between mechanical and thermal states, e.g. between stress and temperature fields.

The latter source, which is related to the thermoelastic effect, is a basic mechanism for thermoelasticity.

Let us begin with the classical thermodynamics of an elastic continuum. Since the solid bodies are considered as an elastic continuum, their behavior can be described by using the concept of reversibility. The first law of thermodynamics simply states that energy can be changed from one form to another but it can neither be created nor destroyed. This principle of energy conservation can be expressed as

dU W Q     (1.29)
where dU is the internal energy, W is the work done on the system by external forces, and Q is the heat exchange from the surroundings to the system. The heat and work processes are dependent of path or transformation from state A to state B, thus Eq. (1.29) is convenient to rewrite corresponding to the infinitesimal change of these quantities per unit mass as

du w q     (1.30)
From the principle of virtual work, the work done on the system by external forces is equal to the strain energy produced by the deformed bodies. Considering an infinitesimal unit volume, the work done on the system by external surroundings can thus be given by the strain energy density:

v ij ij w d     (1.31a)
where  ij denotes the stress tensor and dɛ ij is the small elastic strain tensor. It is necessary to rewrite this expression into a unit mass form, dividing by the density ρ:

ij ij v d w w        (1.31b)
Under the second law of thermodynamics, the heat exchange from the surroundings to the system relating to the specific entropy s can be defined by the following expression:

q Tds   (1.32)
where T is the absolute temperature. Therefore, Eq. (1.30) becomes

ij ij d du w q Tds         (1.33)
Now, the specific Helmholtz free energy function is introduced as a thermodynamic potential [START_REF] Morabito | Analysis of heat sources accompanying the fatigue of 2024 T3 aluminium alloys[END_REF], which is defined by H = u -Ts. Its differential form can be written as

dH du Tds sdT    (1.34)
Combining Eq. (1.33) and Eq. (1.34) yields

ij ij d dH sdT     (1.35)
The Helmholtz free energy H is regarded as a function of strain and temperature [START_REF] Biot | Thermoelastictiy and irreversible thermodynamics[END_REF][START_REF] Johns | Thermal stress analyses[END_REF]. Hence, its differential form is

ij ij H H dH d dT T         (1.36)
Comparison between Eq. (1.35) and Eq. (1.36) with respect to the strain and temperature, provides

ij ij H       (1.37a) H s T    (1.37b)
Simultaneously, the specific entropy s is a state function which is also related to the strain and temperature. Hence,

ij ij s s ds d dT T         (1.38)
The specific heat capacity at zero strain, i.e. in the case of constant volume transformations, is defined as

q Tds s c T dT dT T        (1.
2 ij ij H H ds d dT T T          (1.41a) 1 ij ij q dT d c T T T           (1.41b)
From the classical thermodynamics, the basic relationship between the small temperature changes caused by a change in stress state of a linear, homogeneous, elastic solid and strain can be written in the form:

ij ij T q dT d c T c           (1.42)
During TSA, the test specimen is practically subjected to a cyclic loading at a sufficiently high frequency. As a consequence, the heat transfer between the surroundings and the system in the second term of Eq. (1.42) can be neglected under an assumption of adiabaticity. This equation can

reduces to ij ij T dT d c T        (1.43)
It can be seen that a partial derivative of the stress tensor with respect to temperature exists in Eq. (1.43). In order to develop this equation into a simple form, the constitutive stress-strain-temperature relationships for an isotropic elastic solid [START_REF] Timoshenko | Theory of elasticity[END_REF] are applied:

  2 ij ij kk ij T         (1.44a) and   3 2       (1.44b)
where ɛ kk is the first invariant of the strain tensor, the temperature change is

given by δT = T -T 0 , and α is the coefficient of linear thermal expansion.

Note that T 0 is the absolute ambient temperature at the stress-free. In addition, the Kronecker delta δ ij is equal to 1 for i = j and equal to 0 for i = j.

In general, this Eq. (1.44a) is well-known as the Duhamel-Neumann law [START_REF] Pitarresi | A review of general theory of thermoelastic stress analysis[END_REF][START_REF] Wong | Thermoelastic constant or thermoelastic parameter?[END_REF]. The Lamé constants λ and μ are corresponding to the modulus of elasticity E and Poisson's ratio ν, as follows:

  21 E     (1.44c)    1 1 2 E        (1.44d)
From Eq. (1.44a), the differential form of the stress with respect to the temperature is obtained by

2 ij ij kk ij T T T T T                            (1.45)
Assuming that, for most engineering materials at room temperature the Lamé elastic parameters do not practically vary with the temperature [START_REF] Pitarresi | A review of general theory of thermoelastic stress analysis[END_REF][START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF][START_REF] Wong | Thermoelastic constant or thermoelastic parameter?[END_REF]. So, the temperature derivatives of these properties are neglected. In the case of a stress-induced temperature field, the term of (∂/∂T)δT seems to be a small order of magnitude, which can be negligible compared to . Hence, Eq. (1.45) can be simply written by

ij ij T      (1.46) Substituting Eq. (1.46) into Eq. (1.43) yields kk T dT d c      (1.47)
Integrating Eq. (1.47) between two equilibrium states from 0 to 1, thus

    1 1 0 0 kk kk T kk T dT c d T          (1.48) provides     1 1 0 0 ln kk kk T c T                   (1.49a) 0 ln 1 kk T c T               (1.49b)
Using an infinite power series to expand the natural logarithm term in Eq.

(1.49b) gives

2 3 0 0 0 1 1 2 3 kk T T T c T T T                                                                                        (1.49c)
where T is the temperature change between the two equilibrium states. The order of magnitude of the temperature change T due to a cyclic loading during the TSA experiments is very small (order of 0.001C) compared with the ambient temperature T 0 . As a result, the higher-order terms of the power series can be neglected. Thus,

0 kk T c T        (1.50a) 0 kk c T T         (1.50b)
This relationship also appeared in the pioneer work of Biot [START_REF] Biot | Thermoelastictiy and irreversible thermodynamics[END_REF]. In order to replace the first invariant of the strain tensor ɛ kk by the first invariant of the stress tensor  kk in term of the temperature change, the constitutive stress-strain-temperature relationships for an isotropic elastic solid in Eq. (1.44a) are applied again:

  2 3 3 kk kk T            (1.51a)   3 kk kk T           (1.51b) Substituting Eq. (1.51b) into Eq. (1.50b),   0 3 kk c T T T              (1.52a) 0 3 kk c T T                (1.52b)
Changing the Lamé constants in term of the elastic parameters by their relationships in Eq. (1.44b), Eq. (1.44c), and Eq. (1.44d), yields

  2 0 0 3 12 kk ET T c T                     (1.53)
It can be seen that the bracketed term in above equation mainly consist of the material properties, which vary with the temperature during the experiments. Hence, it would be better to express this equation as only a linear function of the material properties by using the relationship between the specific heat capacity at constant volume and the specific heat capacity at constant pressure [START_REF] Pitarresi | A review of general theory of thermoelastic stress analysis[END_REF][START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF][START_REF] Stanley | Quantitative stress analysis by means of the thermoelastic effect[END_REF]:

  2 0 3 12 p ET c c        (1.54)
Finally, Eq. (1.53) can be rewritten as

0 kk kk p T T A c            (1.55)
This equation is the basic relationship of the thermoelastic stress analysis for an isotropic homogeneous material encountering elastic and adiabatic transformations. In addition, the elastic and thermal properties are assumed to be constant with respect to temperature. It is significant to notice that A is the thermoelastic calibration factor which is a function of the material properties.

Infrared Thermography

Infrared thermography is based on the measurement of the thermal radiation that emits from objects [START_REF] Kylili | Infrared thermography (IRT) applications for building diagnostics: A review[END_REF]. The thermal radiation, which is invisible to the human naked eyes, is measured by the thermal sensor of an infrared camera. Such thermal radiation is then converted into a visual image [START_REF] Clark | Application of infrared thermography to the non-destructive testing of concrete and masonry bridges[END_REF]. It should be noted that other commonly used terms for the nondestructive testing and non-contact diagnostic technology, which is related to IR thermography, are thermal imaging, IR radiometry, IR imaging, and IR condition monitoring [START_REF] Balaras | Infrared thermography for building diagnostics[END_REF]. Apart from thermomechanics concerned in this study, there are a broad range of IR thermography applications [START_REF] Kylili | Infrared thermography (IRT) applications for building diagnostics: A review[END_REF][START_REF] Clark | Application of infrared thermography to the non-destructive testing of concrete and masonry bridges[END_REF][START_REF] Balaras | Infrared thermography for building diagnostics[END_REF][START_REF] Bagavathiappan | Infrared thermography for condition monitoring -A review[END_REF][START_REF] Lahiri | Medical applications of infrared thermography: A review[END_REF]: monitoring civil structures, microscale applications, postharvest quality control, industrial applications, medical science applications, investigation of the natural environment, monitoring of electrical and electronic components, military applications, and so on. How the infrared radiation emits and how the infrared system works are then discussed in the followings.

Infrared Radiation

The thermal radiation is an energy transport mechanism that occurs by means of electromagnetic waves. Such mechanism does not need a medium for transportation. This thermal radiation is generated by a motion of charged particles (protons and electrons), which such motion is caused by the temperature of matter [START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF][START_REF] Maldague | Theory and practice of infrared technology for nondestructive testing[END_REF][START_REF] Rogalski | Infrared devices and techniques[END_REF]. In fact, the thermal radiation can originate from a solid, a liquid, or even a gas because all matters with temperature above absolute zero emit the electromagnetic radiation [START_REF] Bagavathiappan | Infrared thermography for condition monitoring -A review[END_REF][START_REF] Lahiri | Medical applications of infrared thermography: A review[END_REF][START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF][START_REF] Maldague | Theory and practice of infrared technology for nondestructive testing[END_REF]. In 1800, the infrared radiation was discovered by the experimental work of Herschel [START_REF] Herschel | Experiments on the refrangibility of the invisible rays of the sun[END_REF]. His experiment measured the temperature of each color in the spectrum of sunlight transmitting through a prism as shown in Figure 1.17. Herschel built a simple monochromator that used a thermometer as a detector for the temperature measurement. He found that the highest temperature was outside the red in the visible spectrum, which the invisible light is later known as infrared. This discovery of infrared radiation opened a new way in the field of temperature measurement.

Figure 1.17 Experimental study of Herschel [START_REF] Herschel | Experiments on the refrangibility of the invisible rays of the sun[END_REF]: A, B are the small stand, C, D are the prism at window, E is the spectrum of sunlight on the table, and 1, 2, 3 are the thermometers on the stand Figure 1.18 demonstrates the electromagnetic spectrum. Even through the thermal radiation overspreads part of the ultraviolet, all of the visible, and all of infrared band, most of energy of a body at ambient temperature is emitted only in the infrared band [START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF]. In fact, the infrared band covers all electromagnetic spectrum longer than the visible light, but shorter than the microwave. The IR wavelength is positioned in between 0.75-1000 μm.

Many ideas to divide the IR range have been proposed. However, if the division is based on limits of spectral bands of commonly used IR detectors [START_REF] Rogalski | Infrared devices and techniques[END_REF], The IR radiation can be divided into 5 wavelength ranges: Near infrared (NIR) with the wavelength range of 0.78-1 μm is a sensitivity limit windows [START_REF] Kylili | Infrared thermography (IRT) applications for building diagnostics: A review[END_REF][START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF][START_REF] Rogalski | Infrared devices and techniques[END_REF][START_REF] Rogalski | Infrared detectors: status and trends[END_REF]: the middle (MWIR) and long (LWIR)

wavelength infrared spectral bands.

Figure 1.18 Electromagnetic spectrums with wavelength unit in micrometers [START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF] In summary, a major difference between MWIR and LWIR spectral bands can be distinguished [START_REF] Clark | Application of infrared thermography to the non-destructive testing of concrete and masonry bridges[END_REF][START_REF] Rogalski | Infrared devices and techniques[END_REF]. The LWIR band are well-suited to perform with objects at the ambient temperature, while the MWIR band has a great contrast for work with objects at the high temperature. Typically, both LWIR and MWIR spectral bands will generate thermal images when the small temperature changes exist. Nevertheless, the MWIR band displays a better resolution, while the LWIR band may have a higher sensitivity for ambient temperature objects. Under a diversity of weather conditions, the LWIR spectral band has a better performance to operate in the bad weather conditions like fog, dust, and winter haze, contrary to the MWIR band.

In thermal radiation theory, blackbody is considered as a hypothetical object which is a perfect emitter and absorber of radiation [START_REF] Bagavathiappan | Infrared thermography for condition monitoring -A review[END_REF][START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF][START_REF] Maldague | Theory and practice of infrared technology for nondestructive testing[END_REF]. It can be said that the blackbody can absorb all incident radiations from any direction and wavelength. It also emits the maximum feasible quantity of radiation for a given temperature and wavelength. Importantly, no other object can emit more energy than the blackbody at the same temperature and wavelength. In order to evaluate the energy flux (energy rate per unit surface area) per wavelength, the spectral hemispherical emissive power ϕ λ,b

for the blackbody in terms of temperature and wavelength according to Planck's law [START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF][START_REF] Kylili | Infrared thermography (IRT) applications for building diagnostics: A review[END_REF][START_REF] Bagavathiappan | Infrared thermography for condition monitoring -A review[END_REF][START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF][START_REF] Maldague | Theory and practice of infrared technology for nondestructive testing[END_REF] can be obtained by

  2 1 , 5 1 b C T C e       (1.56)
where λ is the wavelength of the radiation and T is the absolute temperature.

The first C 1 and second C 2 radiation constants are given by

2 1 2 C c h   (1.57a) 2 B ch C k  (1.57b)
where c is the speed of light, h is the Planck's constant, and k B is the Boltzmann's constant. The values of these physical constants are given in First radiation constant

C 1 3.741776×10 -16 W•m 2 Second radiation constant C 2 1.438769×10 4 μm•K
The distribution of the spectral hemispherical emissive power for the blackbody as a function of wavelength for a given temperature is presented in Figure 1.19. It is evident that the distribution of spectral emissive power for each wavelength significantly increases with the temperature. Such different curves never cross each other. Furthermore, the peak position of each curve (see a dash line) shifts toward the smaller wavelength when the temperature increases. The locus of maximum wavelength λ max , which the blackbody emits the maximum spectral emissive power as a function of temperature, can be estimated by Wien's displacement law [START_REF] Bagavathiappan | Infrared thermography for condition monitoring -A review[END_REF][START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF][START_REF] Maldague | Theory and practice of infrared technology for nondestructive testing[END_REF]: The total hemispherical emissive power of blackbody at a given temperature is calculated by integrating the Planck's law in Eq. (1.56) over all wavelengths (0 ≤ λ ≤ ∞). It is well-known as Stefan-Boltzman law [START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF][START_REF] Kylili | Infrared thermography (IRT) applications for building diagnostics: A review[END_REF][START_REF] Bagavathiappan | Infrared thermography for condition monitoring -A review[END_REF][START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF][START_REF] Maldague | Theory and practice of infrared technology for nondestructive testing[END_REF], which is expressed as follows.

 

4 , 0 , b b T d BT          (1.59)
where B depends on C 1 and C 2 in term of π 4 C 1 /15(C 2 ) 5 and is known as the Stefan-Boltzmann constant. Its value is given in Table 1.1.

In fact, the radiation of real bodies differs from those of blackbody.

The incident radiation or often called irradiation is perfectly absorbed by the blackbody as shown in Figure 1.20. On the contrary, only a part of the irradiation is absorbed by real bodies, while the remaining part of the irradiation may be reflected and/or transmitted through the real bodies.

Thus,

1 ar     (1.60)
where τ, a, and r are respectively the transmissibility, the absorptivity, and the reflectivity coefficients of materials. In general, all these coefficients are functions of wavelength and direction (azimuth and altitude) [START_REF] Kylili | Infrared thermography (IRT) applications for building diagnostics: A review[END_REF]. Besides the wavelength and direction, the absorptivity coefficient is dependent of the temperature [START_REF] Maldague | Theory and practice of infrared technology for nondestructive testing[END_REF].

Figure 1.20 Incident radiations that impinge on a body [START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF] Most of engineering materials are commonly opaque in the infrared band, even if they are transparent to the visible light such as glass [START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF]. In this manner, the transmitted radiation is equal to zero. Hence, Eq. (1.60) can be rewritten as

1 ar  (1.61)
The Kirchhoff's law of thermal radiation provides the relationship between the absorptivity and emissivity ɛ of materials [START_REF] Maldague | Theory and practice of infrared technology for nondestructive testing[END_REF] as

a   (1.62)
It should be noted that the emissivity is a surface property which indicates an ability of emission energy. It is defined by the ratio of the emitted radiation of real surface to the emitted radiation of blackbody at temperature T and wavelength λ in the direction of azimuth and altitude [START_REF] Kylili | Infrared thermography (IRT) applications for building diagnostics: A review[END_REF][START_REF] Maldague | Theory and practice of infrared technology for nondestructive testing[END_REF]. Hence, it implies that the emissivity is functions of wavelength, temperature, and direction. Considering the blackbody, its emissivity is equal to one which is an upper limit for other real bodies.

In practice, the infrared camera does not only detect the energy emitted from the object but also the energy reflected by the surroundings and intercepted by the object's surface [START_REF] Kylili | Infrared thermography (IRT) applications for building diagnostics: A review[END_REF][START_REF] Fokaides | Application of infrared thermography for the determination of the overall heat transfer coefficient (U-Value) in building envelopes[END_REF]. These quantities are attenuated by the processing of scattering and absorption when the energy transmits through atmosphere [START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF][START_REF] Rogalski | Infrared devices and techniques[END_REF]. Remark that the atmospheric transmission depends on distance, wavelength, and metrological conditions [START_REF] Maldague | Theory and practice of infrared technology for nondestructive testing[END_REF]. The effect of atmospheric transmission can be assumed negligible for the case of camera-to-object distance is small [START_REF] Clark | Application of infrared thermography to the non-destructive testing of concrete and masonry bridges[END_REF]. In addition, the LWIR camera is a better solution to measure the temperature at a long distance of camera-to-object as well as in the case of bad weather conditions. The effect of surrounding reflectance can be solved by stick paper stickers, paint, or electrician's tape on the surface of reflectance object [START_REF] Balaras | Infrared thermography for building diagnostics[END_REF]. In order to achieve the emissivity to be close to one, a commonly used method is to paint a matt-black color on the surface of the object.

Infrared Systems

The devices in IR thermography experiments consist of an infrared camera, a tripod, a display device, and an image processing unit. At the present time, display and image processing are performed in a personal computer with software packages. Generally, the IR camera is composed of optical systems (lens, mirrors, etc.), IR detectors, cooling systems, and associated electronics. The optical systems of IR camera create an image of observed objects in plane of the IR detectors, which is so-called "infrared imaging systems".

Such image can be generated by two main mechanisms [START_REF] Bagavathiappan | Infrared thermography for condition monitoring -A review[END_REF][START_REF] Lahiri | Medical applications of infrared thermography: A review[END_REF][START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF][START_REF] Rogalski | Infrared devices and techniques[END_REF] relying on two types of optical elements: mirrors play a role to reflect the incident radiation and lens to refract and transmit the incident radiation. In general, the IR mirrors are covered by metallic coatings [START_REF] Rogalski | Infrared devices and techniques[END_REF]. There is no a difference in design for the lens of IR camera compared to that of a camera working in the visible band. Nonetheless, there are significantly particular materials suited for IR lens manufacture due to the restriction of spectral band [START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF]. As a result, most of visible lens are opaque in IR spectral band as same as many IR lens are opaque in the visible band. Another optical element that should be introduced is filters. They are used to suppress or transmit certain wavelengths, e.g. enabling the IR camera to either measure the high-level of temperature or measure the temperature of a semitransparent material [START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF].

The first mechanism is scanning system, which is used in the first generation of IR imaging systems. This mechanism consists of a single element IR detector and two scanning mirrors in horizontal and vertical. The scanning mechanism enables a mechanical scanner, i.e. mirrors, to focus different point of the observed objects onto the detectors. Under the scanning system, the second generation of IR imaging systems was then developed by using a large linear array or a small 2D array instead of a single detector. The modern generation of IR imaging systems replaces the scanning mechanism by using "staring array detectors". This staring array has a large 2D array detector, which is popularly known as "focal plane arrays (FPA)". The staring arrays are electronically scanned by using readout integrated circuits (ROICs) with the arrays. The optic elements of IR camera only focus the image onto the detector array.

It could be said that the detector of infrared camera is the core device of the IR system. The infrared detectors are usually classified into two main categories that their physical working principles are different: thermal detectors and photon (quantum or photonic) detectors [194, 196-198, 200, 202]. In the thermal detectors, the electrical output is generated from a result of temperature change of the detectors due to the incident radiation. The thermal detectors are generally independent of the wavelength of irradiation.

They usually operate at the room temperature. The drawback of the thermal detectors is a low sensitivity and a large response time, but they are cheap and easy to use. The most common utilities in IR community are bolometer, pyroelectric detectors, and thermopiles.

In the photon detectors, the incident radiation directly interacts with the electrons of the detector material. The electrical signal causes from the changed electronic energy distribution: the photon can excite electrons and change their energy state if their energy is large enough. The photon detectors can be sub-divided into different types depending on the nature of interaction. The most important ones are intrinsic detectors, extrinsic detectors, photoemissive (metal silicide Schottky barriers) detectors, and quantum well detectors. The photon detectors exhibit a selective wavelength depending on the response per unit incident radiation energy. In addition, they provide a perfect signal-to-noise performance and a very fast response.

However, the photon detectors require extremely cooling down to a stable low temperature, in order to achieve a better performance and reduce a thermal noise generated by detectors. The most ordinary used cooling systems for IR detectors are thermoelectric, bulk, and Stirling coolers. The cooling requirement of photon detectors produces the IR camera to be bulky, heavy, expensive, and difficult to use.

Derivation of Thermoelastic Relationship from Infrared System

Recalling the total hemispherical emissive power of blackbody at a given temperature in Eq. (1.59), this is convenient for the thermal detectors because it covers over all wavelengths. On the other hand, this equation cannot be adopted in the case of photon detectors due to their discrete manner. In this case, it is important to obtain a discrete equation for the number of photons N b emitted by objects at a given temperature by dividing by the energy carried by each photon in each wavelength range [START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF]. Eq.

(1.59) can thus be expressed in term of the total number of photons per unit area and time for the blackbody by

  , 3 0 , 0.370 b b T B N d T E k         (1.63)
where the energy of a single photon given by E = hc/λ and the quantity 0.370B/k is the Stefan-Boltzmann constant for photon detectors equal to

1.52×10 15 photon•s -1 •m -3 •sr -1 •K -3 .
As mentioned earlier, there are two typically used atmospheric windows for IR thermography. It implies that the photon detectors have a relative small working wavelength range on the IR band. As a consequence, Eq. (1.63) should be integrated over only the certain working wavelength range of the photon detectors instead of the range 0 -∞. Clearly, this is impossible to analytically solve. It must be numerically performed by a suggested approximate approach [START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF][START_REF] Dulieu-Barton | A termperature correction methodology for quantitative thermoelastic stress analysis and damage assessment[END_REF]. For the detectors where the working wavelength is less than λ max , the detector response to the temperature change follows an approximate power law, i.e. N b  T n .

Therefore, a relationship between the surface temperature of an object and the total number of photons flux emitted over a specific wavelength range can be expressed with considering the emissivity ɛ as follows:

n N B T     (1.64)
Where B' is a constant depending on the detectors and the index n can be estimated by numerical integrating Eq. (1.63) over the wavelength of interest for various temperatures to determine N b,λ .

In order to develop a relationship between the IR detector voltage output S and the first invariant of the stress tensor , it starts by differentiating Eq. (1.63) with respect to the temperature [START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF][START_REF] Stanley | Quantitative stress analysis by means of the thermoelastic effect[END_REF][START_REF] Dulieu-Barton | A termperature correction methodology for quantitative thermoelastic stress analysis and damage assessment[END_REF][START_REF] Dulieu-Barton | Introduction to thermoelastic stress analysis[END_REF]. Therefore,

1 n N nB T T        (1.65)
Substituting this expression into Eq. (1.65) provides

1 1 n N TA nB T                      (1.66)
Assuming that the detector voltage output S is linearly related to the total photon flux change N λ , i.e. S = ZN λ [START_REF] Dulieu-Barton | Development and applications of thermoelastic stress analysis[END_REF][START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF][START_REF] Stanley | Quantitative stress analysis by means of the thermoelastic effect[END_REF][START_REF] Dulieu-Barton | A termperature correction methodology for quantitative thermoelastic stress analysis and damage assessment[END_REF][START_REF] Dulieu-Barton | Introduction to thermoelastic stress analysis[END_REF], where Z is the detector response factor. Hence, Eq. (1.66) becomes

1 1 1 n S TA nB ZT                 (1.67)
Note that the detector voltage output S is sometimes called "the thermoelastic signal" [START_REF] Dulieu-Barton | A termperature correction methodology for quantitative thermoelastic stress analysis and damage assessment[END_REF][START_REF] Dulieu-Barton | Introduction to thermoelastic stress analysis[END_REF]. The first bracketed term in Eq. (1.67) is referred to as the detector responsivity D, which is specified by the manufacturer of the detector system at a given temperature [START_REF] Dulieu-Barton | Development and applications of thermoelastic stress analysis[END_REF][START_REF] Dulieu-Barton | Introduction to thermoelastic stress analysis[END_REF]. If the IR system is employed at other different temperatures from that given temperature, a temperature correction factor R must be applied. This factor is given by R = (T D /T E ) 2 , where T D is the temperature at the detector responsivity obtained and T E is the temperature of the object surface during experiment [START_REF] Dulieu-Barton | Development and applications of thermoelastic stress analysis[END_REF]. Moreover, it is necessary to introduce an amplification factor F that provides for amplification of the detector voltage output over the full available range and resolution of the analogue-to-digital converter [START_REF] Dulieu-Barton | Development and applications of thermoelastic stress analysis[END_REF][START_REF] Dulieu-Barton | Introduction to thermoelastic stress analysis[END_REF]. Therefore, Eq. (1.67) is then expressed as

DRF S TA            (1.68)
The bracketed term in Eq. (1.68) is a constant, which is dependent of the detector system, the material, and the surface condition. So, this bracketed term can be denoted by a calibration constant A'. Hence, the familiar form of thermoelastic relationship as in Eq. (1.55) is written by

AS    (1.69)
Notice that this expression is valid for a linear elastic, isotropic, homogeneous material as same as the constraint of Eq. (1.55).

Cedip Jade III-MWIR Infrared Camera

The infrared system used in this dissertation is a Cedip Jade III-MWIR, which is manufactured by Cedip Infrared Systems (now FLIR). This is a photon detector with an indium antimonide (InSB) focal plane array of 320 × 240 pixels. This material type of detector is sensitive to the IR wavelength between 3.6 -5.1 μm. The InSb detector requires cooling to about 77 K by means of an internal stirling pump. A summary of the technical characteristic of the Cedip Jade III-MWIR system is specified in

Software

In practice, there are two main different types of software used to control and operate the infrared system. First, Cirrus software is utilized to set IR camera at the beginning of measurements. It enables to focus using a lens, to select the integration time, to specify the frame rate, and to carry out a non-uniformity correction. Considering the integration time, it is similar to exposure time or known as shutter speed in photography that sensitivity to variations increases with the integration time. Significantly, the nonuniformity correction must be applied to images obtained from IR focal plane arrays. Its details will be briefly described later in the topic of camera calibrations.

Second, Altair software is employed to record the thermal data as a sequence of the image at the frame rate that can be set in either Cirrus or Altair software. It is possible to view live in the Altair software. By integrating with Cirrus software, the Altair software can also manage the video configurations and preliminarily post-process the thermal data.

Camera Calibrations

The Cedip Jade III-MWIR infrared camera is equipped with a matrix made of a multitude of independent detector elements. Every detector element sends back a measurement corresponding to a pixel in the produced image. This technology enables to powerfully capture images at high frame rates, and to ensure higher reliability i.e. no mechanical mechanism in motion and each detector are independent of each other. However, additive (offset) and multiplicative (gain) parameters of each individual detector element are different even though they capture the same IR irradiation and are fabricated using the same materials as well as the same technology. Such differences generate non-uniformity in the response of individual pixels in the acquired image, which is also known as fixed noise pattern (FPN). The fixed noise pattern degrades the image quality, the radiometry accuracy, and temperature resolution [START_REF] Mudau | Nonuniformity correction and bad pixel replacement on LWIR and MWIR images[END_REF]. As a consequence, non-uniformity correction (NUC) must be done to compensate the spatial non-uniformity and the temporal variation in the response of detectors [START_REF] Hayat | Statistical algorithm for nonuniformity correction in focal-plane arrays[END_REF]. The NUC is thus necessary to provide a coherent image.

In fact, there are two simple, accurate, and common non-uniformity correction methods that used to correct the non-uniformities: single point correction and two point correction. First, only the offset of every pixel in the IR focal plane arrays is corrected in the single point correction by placing directly a reference source in front of the lens. This reference source should be a non-reflective matt black object under an assumption of homogenous temperature [START_REF] Mudau | Nonuniformity correction and bad pixel replacement on LWIR and MWIR images[END_REF]. Second, the two point correction creates NUC in terms of the offset and gain parameters by using two reference different temperatures at around 30% and 70% of the dynamic range [START_REF]DC001U-E NUC Procedure[END_REF].

These difference temperatures are provided by a black body with extensive surface area. Alternatively, the use of the same reference temperature with two different integration times can be carried out in the case of unavailable two different reference temperatures. Note that more detail about the NUC procedure can be found in ref. [START_REF]DC001U-E NUC Procedure[END_REF].

In practice, the FPN tends to vary slowly in time due to the camera environment [START_REF] Godoy | Noise-cancellation-based nonuniformity correction algorithm for infrared focal-plane arrays[END_REF]. This is a reason that why the single point correction is essential to update every time during the operation of the camera. In the case of two point correction, it will be required to create a new NUC at every time when the value of the integration time or the optical path (lens, window, filter, etc.) changes [START_REF]DC001U-E NUC Procedure[END_REF]. For most commercial IR cameras, the NUC procedure is done during the factory calibration process. The NUC table is then saved in the firmware of IR camera. It is important to notice that the calibration file that was performed during the latest maintenance by the supplier (Cedip Infrared Systems) is used in this study. In addition, the single point correction is performed before beginning the TSA experiments every time.

Apart from non-uniformity of the detectors, a bad pixel is another one parameter that presents defect of the detectors owing to the technology used.

The performance of each pixel in the IR focal plane arrays is characterized by three parameters relating to their response in electrical measurements: its level (or offset), its sensitivity (or gain), and its noise level. These three parameters can be used to determine where the bad pixel in the detectors occurs (more detail in ref. [START_REF]DC003U-A BRP Procedure[END_REF]). In this case, a bad pixel replacement (BRP) procedure is processed, i.e. replacing a bad pixel by one of nondefective neighbors. In fact, the BPR procedure is performed by the camera in real-time after running the NUC process [START_REF] Mudau | Nonuniformity correction and bad pixel replacement on LWIR and MWIR images[END_REF][START_REF]DC003U-A BRP Procedure[END_REF]. Hence, the map of temperature changes is obtained by the T at each pixel. In order to achieve a reliable quatitative stress data at each pixel, a value of the calibration constant (A') is necessary to determine and it thus adopted to multiply the thermoelastic signal S as presented in Eq. (1.69).

Data Acquisition and Post-Processing

There are three standard calibration techniques available for derivation of the value of constant A' [START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF][START_REF] Dulieu-Smith | Alternative calibration techniques for quantitative thermoelastic stress analysis[END_REF]. The first method is a direct calibration using properties of IR detectors, system variables, emissivity of specimen surface, and the thermoelastic constant of the tested specimen (material properties). The second method is a calibration with a measured stress. The last method is a calibration with a calculated stress. All these method are explained as follows.

Firstly, the direct calibration technique is based on solving Eq. (1.68).

It also depends on both material properties and detector parameters to theoretically derive a value of constant A'. These affect this method to be inaccurate, due to an uncertainty in the material properties and the detector parameter. Consequently, the calibration constant A' that derived from the direct calibration method is very sensitive to systematic errors [START_REF] Dulieu-Smith | Alternative calibration techniques for quantitative thermoelastic stress analysis[END_REF].

Secondly, the calibration method against the measured stress requires an independent of the stress responsible for the thermoelastic signal by typical using electric resistance strain gauges. The strain gauges are affixed to the tested specimen in a region of uniform stress where the thermoelastic signal is constant. The IR camera also detects the data on the specimen surface at the same area as the strain gauges applied. In the case of a linear elastic isotropic material under an assumption of plane stress, the sum of principal stresses can be calculated by using the strain measurements obtained from the strain gauges and Hooke's law as

1 kk kk E        (1.70)
Based on this method, the calibration constant A' can be therefore determined by substituting Eq. (1.70) into Eq. (1.69):

 

1 kk E A S       (1.71)
Nevertheless, a drawback of this method is a requirement of both modulus of elasticity and Poisson's ratio of the tested specimen material, inaccuracies of the strain gauges position compared to the measurement point, and feasible errors in the gauge factor [START_REF] Dulieu-Smith | Alternative calibration techniques for quantitative thermoelastic stress analysis[END_REF].

Finally, the calibration method with the calculated stress requires a known stress field. Generally, a tested specimen is subjected to a simple tension. Moreover, the stress can be theoretically calculated by utilizing the cross-sectional area and the applied load. Eq. (1.69) is thus simply to find the calibration constant A'. Significantly, this method has the least potential error sources compared to the others [START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF][START_REF] Dulieu-Smith | Alternative calibration techniques for quantitative thermoelastic stress analysis[END_REF].

As mentioned earlier, the amplitude T of temperature oscillation in the TSA experiments is obtained by a Fourier analysis of the temperature oscillation that converted from the thermoelastic signal S under the calibration process. In practice, the amplitude of temperature oscillation can be used to determine the thermoelastic calibration factor A by solving Eq.

(1.55). The second and third calibration methods are adapted to evaluate the factor A by using the amplitude of temperature oscillation instead of the thermoelastic signal, while the first calibration method is impracticable due to the detector parameters dependence. Importantly, this factor A enables to obtain a corrected stress field as same as the use of the calibration constant A'.

Conclusion

The principles of molecular dynamics (MD) simulation and thermoelastic stress analysis (TSA) technique have been recalled in this chapter. It is important to note that

The MD simulation with a fifth-order Gear predictor-corrector algorithm and the TSA method according to Eq. (1.55) is used in this study. The behaviors of granular solids in terms of granular texture and force distributions were also presented. All these concepts and observations are useful to pave a way for the studies presented in the following chapters of this dissertation. In particular, it is noteworthy that the mechanical behaviors of composite granular materials (constituted from at least two base materials) are rarely mentioned in the literature. This case is treated in the following chapter 3, both from experiments and simulations.

defined as the ratio of the force F applied on a body to the displacement  produced by the force along the axis of compression:

F k   (2.1)
Section 2.2.2 provides some information about an experimental setup.

The results are then presented in section 2.2.3.

Experimental Setup

Principle of Measurement: The value of contact stiffness is determined by compressive testing based on an assumption of plane strain.

In the beginning, a compressive force is applied on the specimen and we measure its displacement simultaneously. After that, the experimental data is utilized to plot a relationship between force and displacement. Finally, the contact stiffness value of the specimen can be calculated from the slope of the plotting between force and displacement: see In order to ensure that the specimen did not move or slip when performing the test, a preliminary compressive loading of about -1 kN were applied on the cylinders before starting the tests.

Results and Discussion

Obtained data were utilized to determine the relationship between force and displacement for each type of cylinder. The contact stiffness value was then determined from the slope between force and displacement. These values will be used later in the MD simulation, in order to calculate the contact forces when two particles are in contact. Let us note that when POM and HDPE are in contacts, the contact stiffness will be given by (k POM k HDPE )/(k POM +k HDPE ). More details will be given later in the following chapters.

Thermoelastic Calibration Factor (A)

Introduction

In the following chapters, the granular specimens will not be painted in black, in order to be able to distinguish POM and HDPE cylinders in optical images. As a consequence, the thermal emissivity of the surfaces is not guarantee to be close to one. It means that the temperature variations measured by IR camera are a priori not correct. Nonetheless, correct stress changes can be correctly estimated by determining the linear relationship between "temperature" and stress in the framework of thermoelastic stress analysis, as shown in Eq. (1.55). Let us recall the different calibration methods in section 1.4.4 of chapter 1. There are three different methods commonly used for calibration [START_REF] Dulieu-Smith | Alternative calibration techniques for quantitative thermoelastic stress analysis[END_REF]: (i) direct calibration from the radiometric properties of the detector, the system variables, the surface emissivity, and thermoelastic constant of the tested specimen, (ii) calibration with the measured stress, and (iii) calibration with the calculated stress. The accuracy of each calibration method is different depending on sources of error as mentioned in previous chapter. It is important to notice that these calibration methods depend on the thermoelastic signal S, as given in Eq.

(1.69). In this study, the derivation of the thermoelastic calibration factor A is thus developed relying on the calculated stress method. The next section will describe the experimental preparation for the measurement of thermoelastic calibration factor A.

Experimental Setup

Principle of Measurement: The thermoelastic calibration factor A can be easily defined if the stress and the "temperature" changes are known: see Eq. (1.55). In practice, a force-controlled cyclic loading is applied to the specimen and the "temperature" changes are simultaneously measured on the specimen surface. The stress changes are basically found from the difference between maximum and minimum forces divided by cross-section area of the specimen.

Specimen preparation, Loading Condition, and

Surface Energy ()

Introduction

In the case of cohesive granular materials, the role of the surface energy is to prevent the contact separation by attractive forces act between solid bodies in intimate contact. In fact, this surface energy of solid bodies can be defined as the energy required to create a unit area of new surface [START_REF] Johnson | Surface energy and the contact of elastic solids[END_REF]. The surface energy is a major parameter to evaluate the adhesive forces in numerical simulations of cohesive granular materials. The obtained magnitude of adhesive forces depends on the models used in the simulations. In the experiments, the surface energy is related to the glue, paints, and vacuum evaporated coatings [START_REF] Johnson | Surface energy and the contact of elastic solids[END_REF]. In our TSA experiments on cohesive non-composite granular materials, the glue was applied at the surface between two particles in contact. Several reasons for glue selection are specified as follows: (i) it does not destroy the material surface, (ii) it is not expensive and easy to buy, (iii) it is transparent when it is dry, and (iv) it is easy to paint on the material surface and also to clean up. In this section, a preliminary experiment is carried out to determine the surface energy of four different types of glue. The preparation of our experiments is presented in section 2.4.2. Results are then given in section 2.4.3.

Experimental Setup

Principle of Measurement: In order to obtain the surface energy, the "separation or threshold tensile" force is required. An easy procedure to acquire the separation force can be done hereinafter under an assumption of plane strain: (i) two cylinders are adjoining with glue, (ii) a uniaxial tensile loading is applied to the cylinder and the displacement is measured at the same time, and (iii) the relationship between force and displacement is established to determine the separation force. Note that the test is finished when the cylinders are separated (disbonding of glue): see Figure 2.8. In fact, there are various 2D and 3D models of the separation force such as JKR (Johnson-Kendall-Roberts) model [START_REF] Johnson | Surface energy and the contact of elastic solids[END_REF], DMT (Derjaguin-Muller-Toporov) model [START_REF] Derjaguin | Effect of contact deformations on the adhesion of particles[END_REF], Maugis model [START_REF] Maugis | Adhesion of spheres: The JKR-DMT transition using a Dugdale model[END_REF], Geometrical Adhesion (GA) model [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF], etc. The GA model, which is the 2D approach used in our MD simulations, is adopted in this study to find the surface energy. In practice, the surface energy of such glues can be evaluated by solving the threshold tensile force in GA model as defined in Eq. (1.25).

Results and Discussion

The separation forces and the surface energy for all specimens are presented in Table 2.3. It is interesting to note that an initial crack of glue has sometimes occurred before the test begins, owing to the grip may be too tight. This situation for instance occurs in specimens 5 for SADER colle bois. The fluctuation of the results might be caused from the quantity of painting glue on the specimens as well as quality of surface of specimens. From the results, the glue named "PATTEX ni clou ni vis" is chosen because its surface energy (1.060 kN/mm) is greater than the others. It is significant to note that the difference between cohesive and non-cohesive granular media is not observed in the simulations when the contact stiffness is higher than the surface energy about fifty times [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF][START_REF] Preechawuttipong | Modélisation du comportment mécanique des milieu granulaires cohésifs[END_REF].

Conclusion

In this chapter, some material properties were obtained by some preliminary experiments:

 Firstly, uniaxial compressive tests on single cylinders provided the value of contact stiffness for the POM and HDPE cylinders employed in the study.

Three diameters of POM cylinder where tested. It was observed that the contact stiffness was nearly the same for the three tested diameters (6.2 mm, 16.4 mm, and 30.9 mm). The mean value for POM cylinders was found equal to k POM = 27.35 kN/mm, which is approximately four times greater than the contact stiffness of HDPE (mean value of k HDPE = 6.43 kN/mm for a unique diameter of 10.5 mm). Let us note that the contact stiffness will be given by (k POM k HDPE )/(k POM +k HDPE ) when POM and HDPE cylinders are in contacts.



Secondly, the thermoelastic calibration factor A for POM and HDPE was investigated by thermoelastic stress analysis experiments. Let us recall that the cylinders are not painted in black, as in the tests presented in next chapter. The value for HDPE was found to be higher than the factor of POM by about two times: A HDPE = 0.0258 K/MPa and A POM = 0.0150 K/MPa.



The last experiment was performed to estimate the surface energy at the contact between two cylinders in a cohesive case. Four types of glue were tested by uniaxial tension. The glue named "PATTEX ni clou ni vis" was selected for the experiment to create adhesive forces in Chapter 6. This glue exhibits a surface energy of 1.060 kN/mm.

All these values will be used later in the next following chapter, e.g. the contact stiffness is required when performing the MD simulations, the thermoelastic calibration factor enable to achieve a reliable quantitative stress data when processing the TSA results, and the surface energy is utilized to calculate an adhesive force in the cohesive granular media when performing the MD simulations.

CHAPTER 3 ANALYSIS OF COHESIONLESS COMPOSITE GRANULAR MATERIALS: TSA EXPERIMENTS AND MD SIMULATIONS

Introduction

Granular materials in the engineering and industrial fields are commonly constructed of various particles with different physical properties. They can be considered as "composite" materials. Most of previous experimental and numerical studies were focused on granular materials made of only one type of constitutive materials. There is clearly a lack of knowledge on composite granular materials. A few previous studies dealt with composite granular media. Some mechanical properties of composite granular materials were predicted by a mathematical modeling [START_REF] Budeshtskii | Mathematical model of granular composite materials[END_REF]. The effect of pulse trapping and disintegration of shock and solitary wave in a onedimensional composite granular protector was experimentally investigated [START_REF] Daraio | Energy trapping and shock disintegration in a composite granular medium[END_REF] and then extended to numerical simulations [START_REF] Wang | Crossover in the power-law behavior of confined energy in a composite granular chain[END_REF]. The force transmission in twodimensional composite granular layers under quasi-static conditions was studied by numerical simulation [START_REF] Yang | Quasi-static response of two-dimensional composite granular layers to a localized force[END_REF]. In the present thesis, the study on composite granular materials under static configurations is performed by both experimental and numerical approaches. On one side, a novel experimental technique for the granular materials community is applied: "Infrared (IR) thermography in a framework of Thermoelastic Stress Analysis (TSA)". On the other side, one of the widely used numerical simulations in the fields of granular materials is also applied in this study: Molecular Dynamics (MD) method.

It should be noted that TSA can provide the sum of principal stresses, which is proportional to the hydrostatic stress. On the contrary, the MD simulations are able to directly give the contact forces. In fact, the contact forces between particles are important microscopic features, which govern the macroscopic behavior of granular materials. The stresses inside the particles can be also considered as significant features at the grain level. They can be for instance employed to predict initial crack or failure occurrence in a granular packing. In the present study, the analysis is performed in terms of hydrostatic stress fields instead of contact forces. As a result, in order to compare the experimental and numerical results, it is necessary to convert the simulate contact force to hydrostatic stresses in the particles.

The aim of this chapter is first to reveal the hydrostatic stress networks in composite granular packing by using TSA technique. Experimental and numerical results are then analyzed through the granular texture. This chapter can be divided into two parts.



The first part deals with the TSA experiments: section 3.2.1 explains the experimental setup; section 3.2.2 describes the post-processing of the TSA results; section 3.2.3 reveals the normalized hydrostatic stress networks in each sample; and section 3.2.4 provides preliminary analysis of the experimental results through the probability distribution of hydrostatic stresses.



The second part of this chapter provides the comparison between TSA results and simulation results through the granular texture: section 3.3.1 gives information about the numerical setup and how to derive the hydrostatic stresses from the contact forces; the comparison between experimental and numerical results is then performed in terms of normalized hydrostatic stress fields in section 3.3.2, in terms of angular distribution of contact directions in section 3.3.3, and in terms of probability distribution of hydrostatic stresses in section 3.3.4.

Thermoelastic Stress Analysis Experiments

Composite granular materials are experimentally studied in this section. Section In practice, there were three stages of loading. First, in order to achieve a static mechanical equilibrium, each granular sample was preliminary compacted by the maximum compression of -60 kN. Second, a waiting time was imposed so that the sample was back into the thermal equilibrium. Finally, cyclic compressive loading at ambient temperature was then applied as a requirement of thermoelastic stress analysis [START_REF] Dulieu-Barton | Development and applications of thermoelastic stress analysis[END_REF][START_REF] Pitarresi | A review of general theory of thermoelastic stress analysis[END_REF][START_REF] Greene | Thermoelastic stress analysis[END_REF][START_REF] Dulieu-Barton | Thermoelastic stress analysis[END_REF].

Loading was force-controlled between -6 kN and -60 kN with a frequency f of 3 Hz. It is important to note that a quasi-static structure of the sample during the cyclic compressive stage was exhibited due to the preliminary compaction at the first stage. In other words, the cylinder displacement was absent during the cyclic loading, except for the deformation. For choosing the magnitude of loading, it was simply estimated by using the basic linear relationship of the TSA as defined in Eq. (1.55). Pragmatically, an expectation of the measurement of temperature changes in the experiments must be greater than the measurement resolution, which was equal to 1 mK based on the supplier datasheet.

A Cedip Jade III-MWIR camera (see Figure 3. After obtained the list of temperature changes, the cylinder positions must be identified because they were necessary to use later for perform the MD simulation, as well as for process the hydrostatic stress fields. In fact, it is very difficult to get the same extracted positions as the actual positions in the experiment, whether identify through thermal or optical images. In practice, the cylinder positions extracted from the thermal image were quite close to the actual positions of each cylinder in the experiments.

Furthermore, the positions extracted from one configuration can be applied to prepare the numerical samples for the MD simulation and to process the hydrostatic stress networks. Therefore, the positions of each cylinder were manually extracted from the thermal image for each configuration, while the optical image was also applied as a reference to distinguish what is type of materials.

It is important to note that cylinder cross sections were not painted in black in order to distinguish between POM and HDPE cylinders in optical image. As a result of non-black painting specimens, thermal emissivity was not guaranteed to be close to one. The thermoelastic calibration factor A, which was determined in section 2. variables should be defined. For considering the whole set of cylinders, the normalized hydrostatic stress  norm in each particle is obtained by the ratio of the mean hydrostatic stress  hyd p in each "all cylinder" as mentioned in above paragraph to the average value of the hydrostatic stress over "all cylinders" in the granular system:

    hyd p norm hyd ALL (3.2)
When only POM cylinders are considered, the definition of the normalized hydrostatic stress in Eq. (3.2) is changed: the mean hydrostatic stress in each "POM cylinder" is divided by the average value of the hydrostatic stress over "POM cylinders" in the system. Likewise, this definition of Eq. (3.2) is changed for considering only HDPE cylinders: the  hyd p in each "HDPE cylinder" is divided by the average value of the  hyd over "HDPE cylinders" in the system. It should be noticed that these definitions of the normalized hydrostatic stress will be also used in the numerical aspect. between POM and HDPE cylinders is also varying as shown in Table 3.1.

The normalized hydrostatic stress fields for each sample as well as their optical pictures before testing captured by CMOS camera are illustrated in In fact, granular materials are discontinuous materials which their quantities in terms of contact forces and particle trajectories are heterogeneous at the grain scale. In order to get the macroscopic description from these microscopic quantities, the homogenization procedure which is referred to as the "representative elementary volume" (REV) is applied [START_REF] Wellmann | Comparison of the macroscopic behavior of granular materials modeled by different constitutive equations on the microscale[END_REF]. The REV is an average area (2D) or volume (3D) for the interesting macroscopic quantities, which is sufficiently accurate to measure the macroscopic quantities [START_REF] Wellmann | Comparison of the macroscopic behavior of granular materials modeled by different constitutive equations on the microscale[END_REF][START_REF] Wellmann | Homogenization of granular material modeled by a three-dimensional discrete element method[END_REF][START_REF] Stroeven | Numerical determination of representative volumes for granular materials[END_REF]. Within the range of the REV, the macroscopic quantities do not change with the size of REV [START_REF] Zhang | Pore scale study of flow in porous media: Scale dependency, REV, and statistical REV[END_REF]. Note that the REV can be defined for each material property or mechanical variables of interest, e.g. solid fraction, coordination number, forces, stresses, strains.

In this study, the solid fraction which is defined by Eq. (1.1) were used as a quantity for discover the size of REV of the granular system. Let see in Figure 3.6a, a determination of REV for sample #11 is performed by using solid fraction. Strong fluctuated value of solid fraction exhibited in the beginning until about 20% of the system area. It seems that a convergence of the solid fraction value is then observed, excluding the system area greater than 90%. In this case, the system area of approximately 72% was selected to perform the statistical analysis. In order to precisely check whether this selected value is appropriate to use as the REV for this sample, the solid fraction was determined by randomly changing position of the selected area in the medium as shown in Figure 3.6b. It is important to mention that other samples also exhibited the same trend as found in sample #11.

In our analysis, grains at the boundary were removed about 2.5 times of the mean diameter from each the boundary side, in order to eliminate the boundary effect and remain enough statistical information. It can be said that approximately 70-80% of the system area, which are dependent on the diameter size, are employed to perform the statistical analysis. This value is sufficient for analysis as mentioned in ref. [START_REF] Radjaï | Force distributions in dense two-dimensional granular systems[END_REF][START_REF] Calvetti | Experimental micromechanical analysis of a 2D granular material: relation between structure evolution and loading path[END_REF]: a size of representative elementary volume (REV) is about ten times of particle diameters in 2D granular assembly. complementary hydrostatic stress networks can be categorized by the average stress as found in the contact forces [START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]: the hydrostatic stresses which is higher than the mean value of the considered type of cylinders called the "strong network", while the "weak network" is defined by the stresses which is lower than the mean. Only the higher stresses in the strong network were considered here in order to avoid the result of low signal-tonoise ratios in our analysis. The higher stresses are also more interesting than the lower stresses: the lower stresses just only maintain the system equilibrium, while the higher stresses play the main role for supporting an external applied loading. The cumulative probability distribution function P( norm ) is given by the ratio between the number of cylinders with a normalized stress greater than  norm and the total number of cylinders. where the coefficient  describes a characteristic of the exponential distribution. It is interesting to note that the distribution law of hydrostatic stresses was the same well-known distribution law found by contact forces [67, 69, 85, 116-118, 121, 123]. A distinction of the hydrostatic stress distributions between monodisperse and polydisperse samples were evidently noticed, expect for sample #3 which the localized hydrostatic stress distributions were exhibited. Likewise, the exponential distribution function was found for both only POM cylinders and only HDPE cylinders consideration as respectively shown in Figure 3.8b and 3.8c. It can be said that the exponential decay is a generic property for granular materials carrying the stresses greater than the mean value. It seems that this property is independent of polydispersity and material types of particle in granular media.

The percentage P 0 of cylinders number carrying the hydrostatic stresses greater than the mean value is presented in Table 3.2. It should be emphasized again that the mean value of hydrostatic stresses corresponds to a consideration of such material types. For instance, "only POM cylinders" consideration means all the stresses over the whole POM cylinders is only used for calculate the mean value. The strong and weak networks were then separated by this mean value. Interestingly, the percentage P 0 for all samples was lower than 50% when all cylinders were considered. This property was also found for all samples in the case of only POM and only HDPE considerations, expect for some cases which was a bit higher than 50%.

Comparison between Experiments and Simulations

In this section, experimental and numerical results are analyzed through the granular texture by using the distributions of the hydrostatic stress in the particles and of the contact direction between particles. Only three samples are used for this purpose: monodisperse configuration in samples #2 as well as polydisperse configurations in sample #8 and #12. There are two main reasons why only three samples were chosen, which they can efficiently represent other samples in terms of particle numbers and particle sizes. First, the effect of monodisperse granular material still exists in the case of ratio between POM and HDPE diameter (D POM /D HDPE ) equals to 1. The positions of each cylinder were manually extracted from the thermal image for each configuration as mentioned in section 3.2.2, while the optical image was also applied as a reference to distinguish what is type of materials. After that, the granular samples were prepared by these positions inside a square box made of four rigid plane walls. Under a quasi-static condition, an incremental compressive loading was then applied on the lower wall up to 50 kN, while the other walls were fixed during this stage.

Although the effect of the gravitational force can be neglected comparing to the external applied loading, the gravity effect was taken into account in this study. The simulations were completed when the granular systems reach the static equilibrium conditions. For dry granular media, the contact forces consist of normal and tangential forces, while the physico-chemical interactions are not considered. The normal and tangential forces in each contact can be calculated by use of Eq. (1.23) and Eq. (1.27), respectively.

The effective contact stiffness k eff component for normal force is dependent on the type of contact between particles. The effective contact stiffness for contact between the same types of material is defined by k POM for POM-POM contact and k HDPE for HDPE-HDPE contact, while the k eff when POM and HDPE particles are in contact is given by (k POM k HDPE )/(k POM +k HDPE ).

The value of wall contact stiffness is simply assumed as same as the value of k POM .

Figure 3.9 Scheme of numerical sample preparations

The list of contact forces at the final stage was acquired. The stress tensor  p in an individual cylinder as defined in Eq. (1.11) was then calculated by using the normal and tangential forces in each contact. Under assumption of plane stress, the two in-plane principal stresses  1 and  2 of  p are determined for each cylinder in order to calculate the hydrostatic stress or mean stress [START_REF] Cambou | Micromechanics of Granular Materials[END_REF] defined by:

  1 2 2      hyd (3.3)

Normalized Hydrostatic Stress Fields

The MD simulations provide the network of normal contact forces for sample #2, #8, and #12 as shown in Figure 3.10c, while the fields of normalized hydrostatic stresses  norm obtained from the TSA experiments and calculated from simulated contact forces are displayed in Figure 3.10b and 3.10d, respectively. The normal contact forces were represented by the red lines, whose thickness is proportional to the magnitude of the normal force. An inhomogeneous contact force network was exhibited for both monodisperse (sample #2) and polydisperse medium (sample #8 and #12) as shown in Figure 3.10c, similar to the experimental and numerical hydrostatic stress networks. It can be seen that a triangular network in both terms of the hydrostatic stresses and the contact forces is found for sample #2. This is a fact that it is a common feature for monodisperse granular materials. The highest stresses were observed in the small cylinders in sample #12, even though the large particles transmitted the strong forces. This argument can be explained by Eq. (1.11): the division by V p leads to high stress values when the diameter is small, even if the contact forces are low.

It is not possible to have precisely the same stress pattern between the experiment and the simulation at a microscopic point of view, due to a simple fact that the simulation can perfectly do everything as you desire. For example, the gravitational effect can be ignored in the simulation despite the fact that it appears in real environment. Otherwise, the particle shape in the simulation is flawless: complete roundness for a circular particle which is not feasible for real object. The statistical approach regarding textural properties of granular materials [START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF][START_REF] Cambou | Micromechanics of Granular Materials[END_REF] is thus performed in order to compare the macroscopic mechanical behaviors obtained from the TSA experiment and the MD simulation. In fact, the behavior of granular materials at the global scale is attributed to an effect of contact forces at the local scale.

Angular Distribution of Contact Directions

Let us begin to statistically analyze our experimental and numerical results through the angular distribution of contact directions, which is sometimes called "geometrical fabric" [START_REF] Andreotti | Granular media: between fluid and solid[END_REF]. In order to clearly characterize and describe a chance of contact occurrence in each direction of the space, the angular distribution of contact directions is presented. The angular distribution of contact directions provides interesting information about the spatial arrangement of the particles and enables us to reveal a buildup of anisotropic structure of the granular system. It is worthy to note that the "contact direction" is the spatial aspect of the line between the centroids of two particles in contact, which is measured with respect to the horizontal.

The contact detection between two particles in the simulations is underlying on a criterion is that the contact is discovered when the sum of their radius is greater than the distance between their centroids. In the experiments, the contact directions are determined through the stress map together with the cylinder centers under the same criterion as in the simulation. Second, let us considers the contact network between POM-POM particles in Figure 3.11b. This type of contact revealed the anisotropic distribution for both the experiments and the simulations, even though the distribution of contact directions obtained from the simulations for sample #8 was slightly more anisotropic than the experimental contact distribution.

In addition, the contact directions achieved from the experiments in sample #12 was a bit more anisotropic than the numerical result of the contact distribution. It can be observed that the distribution of the POM-POM contacts in sample #12 is inclined of about 45 from the vertical axis (Figure 3.11b), contrary to sample #8. The reason is that the cross-sections area between POM and HDPE cylinders is much different. In addition, this sample consists of a high number of HDPE cylinders compared to the number of POM cylinders. Hence, the contact directions in Figure 3.11c tend to be distributed along specific directions as in the monodisperse case.

In fact, the distribution of POM-POM contacts (the higher stiffness) attempts to arrange itself along the compression axis. Based on the effect between the orientations of contacts between POM-POM and HDPE-HDPE in this sample, that is why the contact orientations between POM-POM particles in sample #12 is inclined of about 45 from the vertical axis.

From these several comments, it can be said that the results of the development of anisotropic structure due to mechanical loading obtained from the TSA experiments and the MD simulations are good correlated in overview.

Probability Distribution of Hydrostatic Stresses

All the experimental results were already analyzed through the statistical distribution of hydrostatic stresses in section 3.3.2. However, this analysis is done again in this section. Only samples #2, #8, and #12 are here considered. The results obtained from the TSA experiments and the MD simulations were then compared. Figure 3.12 presents the distribution of normalized hydrostatic stresses in the strong network, which was plotted in semi-logarithmic scale. Three different types of cylinder were separated for consideration: all cylinders, only POM cylinders, and only HDPE cylinders.

From Figure 3.12b, all both experimental and numerical data were evidently distributed as an exponential decay predicted by Eq. (3.2). Other interesting observation was that less than 50% of cylinders carrying the hydrostatic stress greater than the mean value as shown in Figure 3.12a.

Both observations were identical for all considerations: all cylinders, only experiments and molecular dynamics method for the simulations. The composite samples were prepared using two constitutive materials: polyoxymethylene (POM) and high-density polyethylene (HDPE). The contact stiffness of the former is about four times higher than that of the latter. Experiments and simulations provide similar results in a global point of view:



The distribution law for the hydrostatic stresses exhibits an exponential decay. This is consistent with the well-known distribution law found for the contact forces [67, 69, 85, 116-118, 121, 123].



The amount of cylinders carrying a hydrostatic stress higher than the mean value is about less than 50% of the quantity of considered material (all cylinders, only POM cylinders, and only HDPE cylinders). This property partly corresponds to the bimodal character of granular media composed of only one type of base materials [START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Cambou | Micromechanics of Granular Materials[END_REF][START_REF] Azéma | Force transmission in a packing of pentagonal particles[END_REF][START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF].



Considering in terms of organization of the particles and their contacts in space, the quasi-homogeneous external applied loading affects to the orientations of the whole set of contacts becomes quasi-homogeneous (isotropic) manner. When each type of contacts is separately considered, the contact network between POM-POM cylinders is preferably oriented along the axis of compression. It can be said that an external applied loading is mainly carried by the POM-POM contacts, while the other contact networks, i.e. between HDPE-HDPE cylinders and between POM-HDPE cylinders, just only make the granular system in equilibrium by supporting the contact network between POM-POM cylinders.

In conclusion, infrared thermography is a novel full-field experimental technique applied to measure stresses in granular materials. Let us recall (see general introduction) that a few studies on granular media were based on IR thermography: measurement of mechanical dissipation in sands and soils [START_REF] Luong | Charateristic threshold and infrared vibrothermography of sand[END_REF][START_REF] Luong | Introducing infrared thermography in soil dynamics[END_REF], characterization of heat transfer in moving spheres [START_REF] Sepúlveda | Infrared tracking from morphological image processing tools -application to heat transfer characterization in granular media[END_REF], and estimation of soil surface microrelief and rill morphology [START_REF] De Lima | Can infrared thermography be used to estimate soil surface microrelief and rill morphology?[END_REF]. But IR thermography has never been used for reveal the stress network in granular assembly. Our approach enables us to reveal the stress transmission in granular materials at the grain scale and opens new possibilities for studying further these materials. Importantly, IR thermography is of course complementary to the other fullfield measurement techniques and numerical simulations also.

It can be noted that the influence of the ratio of particle numbers on the mechanical behavior of non-cohesive composite granular materials cannot be easily analyzed using our experimental technique, due to the high number of particles to be prepared. It can be also noted that the study of the influence of the contact stiffness ratio is also difficult to experimentally analyze. As a consequence, both effects (influence of the ratio of particle numbers and of the ratio of the contact stiffness) are numerically analyzed in chapter 4 and chapter 5, respectively.

randomly distributed inside a square box consisting of four rigid plane walls. The ratio between the number of HDPE and POM particles (N HDPE /N POM ) was progressively changed from 1 to 10, while the ratio between POM and HDPE diameters (D POM /D HDPE ) was fixed at 3.0 for the ten configurations. It means that the simulation must be conducted at least 10 times, because the ratio of particle number is fixed at one of given values during perform the simulation. The contact stiffness of POM was about four times greater than the contact stiffness of HDPE which was determined in chapter 2:

27.35 kN/mm and 6.43 kN/mm, respectively. The effective contact stiffness k eff component for normal force calculation relies on the contact types between particles.

The effective contact stiffness for contacts between the same types of material is defined by k POM for contact between POM and POM particles and k HDPE for contact between HDPE and HDPE particles, while the k eff when POM and HDPE particles are in contact is given by (k POM k HDPE )/(k POM +k HDPE ). The value of wall contact stiffness is simply assumed as same as the value of k POM . In order to prepare an initial granular packing, the numerical samples were subjected to a compressive loading of 100 N on the upper wall as shown in Figure 4.1a, while the other walls were fixed during this stage. The gravity effect was not considered in the stage due to desire a dense system by only an external applied loading. Our numerical samples reached a quasi-static equilibrium condition at the end of preparation stage. The compression stage was then performed by using the quasi-static samples from the preparation stage. An incremental compressive loading was applied on the lower wall up to 50 kN as illustrated in Figure 4.1b. The other walls were fixed as same as the preparing stage. The gravitational force was taken into account in this stage.

Considering the contact forces, they consist of normal and tangential forces excluding the physic-chemical reaction. The normal and tangential forces in each contact were calculated by use of eq. (1.23) and eq. (1.27), respectively. The simulations were completed when the granular systems reached quasi-static equilibrium conditions.

Results and Analysis

In order to understand the influence of the N HDPE /N POM ratio, the simulation results are here analyzed by the granular texture [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF][START_REF] Cambou | Micromechanics of Granular Materials[END_REF] in term of coordination number, probability distribution of contact forces, and angular distribution of contact directions. As mentioned in chapter 1, the coordination number z is defined as the average number of contacts per particle. Figure 4.3 shows the influence of the particle number ratios on this number z of the final of compression stage. For all contact networks, z slightly increases when N HDPE /N POM increases (see the black solid line). The minimum and maximum value is about 4.3 and 4.9, respectively. In fact, the two-dimensional granular system is quite compact if the value of z exceeds 4, except for the monodisperse case which has a value of z close to 6. The trend of z for all contacts is to reach the value of z about 6 as in a monodisperse case, when the N HDPE /N POM ratio is high enough. Considering HDPE particles (blue solid line), z greatly increases when N HDPE /N POM increases, while the value of z for POM particles tend to decrease (red solid line). This is the fact that there is higher chance of contact occurrence between HDPE and HDPE particles when the number of HDPE particle is increasing. This situation also reduces the chance of POM-POM contact occurrence. It should be noted that the coordination number for each material type can be calculated by the ratio between two times of number of contacts between considered material type and number of considered material, which is based on Eq. (1.2). Let us consider POM and HDPE particles are in contacts with any material types, which can be determined by two times of contact number between considered material type and any material types divided by number of considered material. This expression is relevant to Eq. (1.2). In this analysis, any material types mean both POM and HDPE particles. It can be observed that the coordination number z for contacts between HDPE and any particles (see the blue dash line) increase with the N HDPE /N POM ratio as same as z for contacts between POM and any particles (red dash line), excluding N HDPE /N POM = 10. For each N HDPE /N POM , the coordination number z for POM-ANY contacts is higher than z for HDPE-ANY contacts. It can be said that an occurrence of contacts between POM particles and any types of particle has higher chance than contacts between HDPE particles and any particle types. This is based on the fact that the effect of the particle size, which the diameter size of POM is about three times higher than the diameter size of HDPE.

Probability Distribution of Contact Forces

Due to an inhomogeneous distribution of contact forces, it is interesting to consider the effect of this feature by a statistic of an array of contact forces concerning their magnitude. It is important to note that the array of contact forces can be split into the weak and the strong networks [START_REF] Cambou | Micromechanics of Granular Materials[END_REF][START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]. The former consists of the contacts which transmit a force lower than the mean force <F>, contrary to the latter that the contacts carry a force higher than the mean value. In this section, the contact forces in the strong network are only focused: the forces in the strong network play a role to carry an external applied loading, while the stability of the granular packing is ensured by the weak forces. This is also observed in the POM-POM, HDPE-HDPE, and POM-HDPE contact networks with respect to the number of their own contact of each type (dash red, blue, and green line respectively). This property partially relates to the bimodal character as generally found in granular media made of only one type of constitutive material [START_REF] Cambou | Micromechanics of Granular Materials[END_REF][START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]. It can be said that the property does not depend on the N HDPE /N POM ratio or the type of contacts.

The semi-logarithmic plot of the probability distributions P of normal forces in the strong network in case of N HDPE /N POM = 1 is shown in Figure In order to clearly characterize and describe the probability of contact occurrence in each direction of the space, the angular distribution of contact directions which is sometimes called "geometrical fabric" [START_REF] Andreotti | Granular media: between fluid and solid[END_REF] is presented.

This angular distribution of contact directions provides interesting information about the spatial arrangement of the particles and enables us to reveal preferred orientation of the granular system. In fact, the direction of contact is determined by the line connecting between the centroids of two particles in contact, which is measured with respect to the horizontal. The contact detection between two particles in the simulations is based on a criterion: the contact is detected when the sum of their radius is greater than the distance between their centroids. The polar diagram of the number of contacts in the normal direction is here presented in The whole contact networks and HDPE-HDPE contact networks are oriented with a quasi-homogeneous manner for a N HDPE /N POM ratio lower than 8. On the contrary, the particular directions occur as in monodisperse granular media for a N HDPE /N POM ratio greater than 8 (Figure 4.6a and 4.6c).

Note that the POM-POM contacts in Figure 4.6b are rearranged along the compression axis. Figure 4.6d shows a quasi-homogeneous distribution, which is observed on the POM-HDPE contacts for all the ratios between particle numbers.

Conclusion

In this chapter, the influence of the ratio of the particle number on twodimensional composite granular media without cohesion has been studied using numerical simulation based on the molecular dynamics method. The results show that the contact forces greater than the mean force (referred as "strong network") corresponds to about 40% of the whole contacts. It is also observed in POM-POM, HDPE-HDPE, and POM-HDPE contact networks with respect to the number of POM-POM, HDPE-HDPE, and POM-HDPE contacts. The POM-POM contacts, which carry the external applied loading, are oriented along the compression axis. The HDPE-HDPE and POM-HDPE contacts, which make the granular system in equilibrium, exhibit an isotropy distribution (except for some cases of HDPE-HDPE contacts which the particle sizes are close to each other). It is worthy to conclude that the ratio between the particle numbers affects to the granular medium through analysis of the granular texture. The property partly correlates with the bimodal character as in non-cohesive granular materials made of only one type of constitutive material [START_REF] Cambou | Micromechanics of Granular Materials[END_REF][START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF] is also found. the ratio of the particle number (N A /N B ) was fixed at 1.0 for the ten configurations. It means that the simulation must be performed at least 10 times, because the ratio of contact stiffness is fixed at one of given values during perform the simulation. In fact, the ratio of the maximum diameter size to the minimum diameter size (D max /D min ) was set to 2.0 in order to avoid the effect of triangular contact network as generally found in monodisperse case. The effective contact stiffness k eff component for normal force calculation relies on the contact types between particles. The effective contact stiffness for contacts between the same types of material is defined by k A for contact between particles A-A and k B for contact between particles B-B, while the k eff when A and B particles are in contact is given by by (k A k B )/(k A +k B ). The value of wall contact stiffness is simply assumed as same as the value of contact stiffness for particles A.

In order to perform the simulation as shown in

Results and Analysis

The granular sample with k A /k B = 4 at the equilibrium state is presented in Figure 5.2. Particle A is represented by a black circle, while the cyan circle refers to particle B.

The normal forces are demonstrated by the solid red line connecting between the centroid of particles. The thickness of that line is proportional to the magnitude of normal force. It can be observed that the inhomogeneous contact force distributions are exhibited as shown in Figure 5.2b. In order to understand the effect of the contact stiffness ratio, the simulation results are here analyzed by the granular texture [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF][START_REF] Cambou | Micromechanics of Granular Materials[END_REF] in term of coordination number, probability distribution of contact forces, angular distribution of contact directions, and angular distribution of contact forces. A force in granular media is transmitted through inter-particle contacts, which leads to an inhomogeneous distribution of contact forces [START_REF] Radjaï | Force distributions in dense two-dimensional granular systems[END_REF][START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]. In fact, the contact force networks can be separated by the average force <F> into the two complementary networks [START_REF] Cambou | Micromechanics of Granular Materials[END_REF][START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]: the contacts carrying a force greater than the average force called the "strong network", while the "weak network" is defined by the contacts carrying a force smaller than the average force. Only the forces in the strong network were considered in this section, because the strong network is more interesting than the weak network: the forces in the weak network just only maintain the system equilibrium, while the forces in the strong network play the main role for supporting an external applied loading. Interesting information about the spatial arrangement of the particles and preferred orientation of the granular system can be investigated by the angular distribution of contact directions, which is sometimes called "geometrical fabric" [START_REF] Andreotti | Granular media: between fluid and solid[END_REF]. It is important to mention that the orientation of each contact is measured by the line connecting the centroids between two particles in contact with respect to the horizontal. In the simulations, the contact between two particles is detected by a criterion is that the contact is discovered when the sum of their radius is greater than the distance between their centroids. It should be noted that the magnitude of the forces in each contact is not considered here. It is contradictory with the "angular distribution of contact forces", which will be thoroughly described in the whole contact network is also observed in the contact networks between materials B-B in Fig. 4c and A-B in Fig. 4d, while the higher stiffness contacts (A-A contacts) in Fig. 4b appear to be slightly rearranged along the compression axis. Note that it would be quite difficult to distinguish the difference of the contact distributions for all the contact types, because preparing samples in the square box subjected to a uniaxial compressive load makes our granular systems to be homogeneity.

Angular Distribution of Contact Forces

The influence of the contact stiffness ratio on the geometrical fabric of the medium is not completely clear (see previous section). It is therefore interesting to analyze the results by separate the contact networks into two complementary zones: the strong and the weak networks as mentioned in section 5.3.2. The contact orientations of strong and weak networks are then examined by the angular distribution of contact forces which is well-known as the "mechanical fabric" [START_REF] Andreotti | Granular media: between fluid and solid[END_REF]. For the mechanical fabric, the measurement of contact directions and the detection of contacts are the same produce as the geometrical fabric, which was explained in the previous section.

However, a difference is that the magnitude of the normal contact forces is taken into account for the mechanical fabric. In order to separate the contact networks between the strong and the weak networks, the magnitude of For all the tested k A /k B ratios and all the contact types, the contact network corresponding to strong forces (strong network) was more anisotropic than the contact network corresponding to weak forces (weak network). The contact orientation of the strong network was rearranged along the axis of compression, while the isotropic distribution was obviously exhibited in the weak network. On the contrary, the isotropic distribution in the weak network was not observed for the contacts between materials A-A and B-B.

Conclusion

The effect of the contact stiffness ratio in cohesionless composite granular materials was examined. For this purpose, the two-dimensional polydisperse granular samples were studied by using molecular dynamics simulation relied on the discrete element method. The results can be described as follows. Bimodal character of the composite granular materials is exhibited as commonly found in non-composite granular media [START_REF] Cambou | Micromechanics of Granular Materials[END_REF][START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]. The bimodality of granular materials can be represented by the strong network which is corresponding to about 40% of all contact forces greater than the mean force. The strong network is also oriented along the axis of compression.

The 40% component of forces higher than the mean is also observed in the contacts between materials A-A (higher stiffness), B-B, and A-B with respect to the number of their own contact types. The force distributions in strong network are exhibited as an exponential decay. Moreover, the strong contact network is more anisotropic than the weak contact network. By separating the type of contacts, the contacts of higher stiffness materials, which carry the external applied loading, are oriented along the compression axis. The B-B and A-B contacts, which act to make the granular system reach equilibrium, exhibit an isotropic distribution.

CHAPTER 6

ADDITIONAL STUDY: COHESIVE GRANULAR MEDIA

Introduction

This chapter is aimed at preliminarily examining the hydrostatic stress network in cohesive non-composite granular media by means of infrared thermography within the thermoelastic stress analysis principle. The molecular dynamics simulation is also used as a supplementary tool for this purpose. Recalling our TSA experiments on composite granular materials in chapter 3, the hydrostatic stress networks at the contact scale for non-cohesive granular sample were evidenced. This inspires to perform the TSA experiments on cohesive granular media. In addition, the previous studies on the behavior of cohesive granular material were quite restrictive. Most of them have been performed by numerical simulations [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF][START_REF] Richefeu | Force transmission in dry and wet granular media[END_REF][START_REF] Radjaï | Cohesive granular texture[END_REF][START_REF] Preechawuttipong | Microscopic features of cohesive granular media[END_REF][START_REF] D'addetta | On the application of a discrete model to the fracture process of cohesive granular materials[END_REF][START_REF] Rognon | Dense flows of cohesive granular materials[END_REF][START_REF] Hidalgo | Granular packings of cohesive elongated particles[END_REF][START_REF] Nase | Discrete characterization tools for cohesive granular material[END_REF], while a few of experimental techniques is applied [START_REF] Nase | Discrete characterization tools for cohesive granular material[END_REF][START_REF] Alexander | Avalanching flow of cohesive powders[END_REF][START_REF] Samadani | Angle of repose and segregation in cohesive granular matter[END_REF][START_REF] Delenne | Mechanical behaviour and failure of cohesive granular materials[END_REF]. The outline of this chapter is organized as follows. The first section presents experimental and numerical preparations. The experimental and numerical results are then discussed in the second section.

Methodology

In dry or non-cohesive granular materials, the macroscopic mechanical strength of these materials is insured by a confining pressure which can be described by an elastoplastic behavior. On the other hand, no confinement is required for the mechanical integrity of cohesive granular materials due to the existence of cohesive interactions.

The Mohr-Coulomb criterion is classically and extensively used model for the mechanical strength of cohesive granular materials [12]. The cohesion of granular materials arises from attraction forces between particles, which have various physical and chemical origins. These interactions restrain the loss of contacts and elastic contact deformation up to the point of particle separation [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF][START_REF] Preechawuttipong | Microscopic features of cohesive granular media[END_REF]. It can be said that the cohesive interactions between two particles are able to resist not only normal separation, but also shear and rolling. The cohesive interactions between two particles may occur from long range interactions between the mass centers (e.g. colloidal forces for small particles or gravitational forces for asteroids), surface interactions (e.g. through surface energy or liquid bridge at the particle contacts or in the case of physico-chemical interactions with very short range), and a third body (e.g. cement for partially fills the pore space) [12]. In our experiment, glue is used for bonding particles, which allows occurrence of the contact or the adhesion between particles. In addition, the cohesive forces between two particles are irreversible: the ability to support tension of the cohesive interactions is permanently lost when the bonding of glue was failure, even if both particles return in contact again. This failure or separation only occurred in the bonded joint, which does not destroy the surface of particles. It is worth noting that the cohesive interactions in our simulation are reversible according to the forces model. We consider experimental and numerical preparations for study on cohesive granular media in the following sections.

Experimental Setup

In practice, high enough magnitude of a mechanical loading is needed to apply on granular samples in order to ensure that the stress network is evident. On the contrary, the failure of cohesive bonding between interparticle contacts owing to an external applied loading is carefully taken into account. To avoid this failure, our cohesive granular sample was placed at the center of a granular system surrounded by non-cohesive grains as shown in Figure 6.1. The 2D granular system was composed of 662 particles made of polyoxymethylene (POM) cylinder with 60 mm in length [START_REF] Calvetti | Experimental micromechanical analysis of a 2D granular material: relation between structure evolution and loading path[END_REF]. Two different sizes of diameter were used: 12.4 mm and 16.4 mm, while the ratio of particle number between small and big diameter sizes was approximately set to 2.0. The cohesive agglomerate, which consists of about 77 particles, was prepared by using "PATTEX" glue as cohesive bonding. Several reasons to select this glue were already explained in section 2.4, chapter 2.

The particles were randomly placed inside a square metallic frame. Figure 6.2 describes the preparation of the cohesive specimen. First, a cylinder was painted in glue on the black shading area. This cylinder was then attached with another cylinder. After the glue had dried, it was painted on other cylinders and attached to the collection of cylinders already glued. All previous steps were performed again until the size of cohesive agglomerate was reached. This preparation was based on the assumption that the occurrence of the cohesive bonding between two particles is complete and uniform, without any breakage before testing.

Numerical Setup

As the same procedure in chapter 3, the position of each cylinder was manually extracted from the thermal image obtained by TSA experiments.

These positions were used for prepare the numerical sample inside a square box made of four rigid plane walls. The numerical model for the sample was subjected to an incremental compressive loading up to -40 kN on the lower wall, while the other walls were fixed during this stage. The gravitational force was considered although its effect can be ignored comparing to the external applied loading. The simulations were completed when the granular systems reached the static equilibrium conditions. Considering the contact forces for cohesive granular media, they consisted of normal and tangential forces including the physico-chemical interactions. Based on the geometrical adhesion model, the normal and tangential forces in each contact can be calculated by using of Eq. (1.24) and Eq. (1.28), respectively.

The surface energy  was already determined in section 2.4, chapter 2.

Results and Analysis

The experimental and numerical results in term of hydrostatic stress field are presented in this section. A statistical analysis through probability distribution of hydrostatic stresses and angular distribution of contact directions is then performed.

Although the simulation directly provides the contact forces, the hydrostatic stresses are used in this analysis in the same way as in chapter 3. The in-plane principal stresses are first determined by the list of contact forces at the final stage, and then the hydrostatic stress is calculated by using Eq. (3.3).

Normalized Hydrostatic Stress Field

The contact force networks, which were obtained from the MD simulation, are illustrated in Figure 6.3b. The lines, which connect between the centers of the particles, represent the forces transmitted through the contacts. The line thickness is proportional to the magnitude of the normal force. The red lines indicate compressive efforts, while tensile efforts are presented by the green ones. The green lines are not observed in Figure 6.3b, because of the strong compressive applied loading. However, few contacts exhibited the tensile efforts with very small magnitude. It can be seen that a heterogeneous contact force network was evidently displayed.

Furthermore, the strong forces seemed to be transmitted particularly by the contacts between dry particles.

Under the framework of thermoelastic stress analysis, the hydrostatic stresses can be calculated by using Eq. (3.1). The normalized hydrostatic stresses  norm fields obtained from the TSA experiment on the one hand, and calculated from simulated contact forces on the other hand, are shown in Figure 6.3a and 6.3b, respectively. The normalized hydrostatic stress is obtained by the ratio of the hydrostatic stress in each particle to the mean value over the whole granular sample. For both experimental and numerical results, an inhomogeneous stress network was obviously exhibited, as in the case of polydisperse granular materials without adhesion. It can be observed that the normalized hydrostatic stress paths were principally distributed along the vertical direction, which was the same direction as an external applied load. By comparison between cohesive and non-cohesive zones, it seems that an arrangement of the main stress paths was also mostly away from the cohesive agglomerate. A micromechanical analysis of these experimental and numerical results will then be performed through the granular texture. The statistical analysis by probability distribution of hydrostatic stresses is now carried out in this section. In order to restrict the effect of low signal-to-noise ratios from the stresses lower than the mean, the stresses greater than the mean value is only considered. The probability distribution function P of the normalized hydrostatic stresses greater than the mean, which is referred to as the "strong network", is displayed in Figure 6.4a (linear scale) and in Figure 6.4b (semi-logarithmic scale). The probability distribution function P is given by the ratio between the number of cylinders carrying a stress higher than the mean and the total number of cylinders.

Results take into account in turn only non-cohesive cylinders, only cohesive cylinders, and all cylinders. It is interesting to notice that the distribution of function P can be well fitted by an exponential decay for all cases obtained from the TSA experiment and the MD simulation. This behavior was also discovered by numerical investigations of cohesive granular materials [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF][START_REF] Richefeu | Force transmission in dry and wet granular media[END_REF][START_REF] Radjaï | Cohesive granular texture[END_REF]. It can be said that the distribution law in the strong network is a general characteristic of granular media regardless of types of interactions, i.e. cohesion or cohesionless. From Table 6.1, let us consider the percentage of cylinder number in the strong network for both the experimental and the numerical results. The number of cylinders carrying stresses greater than the mean value was about less than 50% of the whole cylinders. It was also observed when only non-cohesive cylinders and only cohesive cylinders were considering. In addition, it can be seen that the percentage of cylinder number in the strong network between the whole cylinders and only dry cylinders is somewhat close to each other. It is interesting to notice that the number of cohesive cylinders would be increased for the further studies in order to achieve a better analysis. The property, which was observed here, was usually found in non-cohesive granular materials by means of numerical and experimental approaches [START_REF] Radjaï | Force distributions in dense two-dimensional granular systems[END_REF][START_REF] Mueth | Force distribution in a granular medium[END_REF][START_REF] Løvoll | Force measurements on static granular materials[END_REF]. This section provides the analysis of our experimental and numerical results through the angular distribution of the contact directions, which is sometimes called "geometrical fabric" [START_REF] Andreotti | Granular media: between fluid and solid[END_REF]. The angular distribution of contact directions enables to characterize and describe the spatial different considerations of type of contacts were performed: the whole set of contacts in Figure 6.5a, the contacts between two dry particles in Figure 6.5b, the contacts between two cohesive contacts in Figure 6.5c, and the contacts between dry and cohesive particles in Figure 6.5d. The whole set of contacts consists of dry-dry, cohesive-cohesive, and dry-cohesive contacts.

From a global point of view, it can be said that the geometrical fabric obtained from the TSA experiment and the MD simulation are in good agreement except for the case of contacts between dry and cohesive particles in Figure 6.5d. Considering the whole set of contacts in Figure 6.5a and the contacts between two dry particles in Figure 6.5b, the arrangement of their structures was appeared on two particular directions (45) with respect to the horizontal axis. The anisotropic structure also occurred in the case of contacts between two cohesive particles as shown in Figure 6.5c, even though it was not completely clear to see the two specific directions.

These specific directions were initially discovered in the granular media with strong adhesion by means of numerical simulations [START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF].

Conclusion

Infrared thermography based on thermoelastic stress analysis was applied to study on two-dimensional polydisperse cohesive granular media. The molecular dynamics simulation was also used for this purpose. In order to avoid the failure of cohesive bonding, the cohesive sample was set in the middle of the granular system surrounded by non-cohesive grains. Cohesive bonding was represented in glue. The experimental results show that approximately less than 50% of all the cylinders carrying the stresses higher than the mean value are observed. Likewise, this property is discovered when only non-cohesive cylinders and only cohesive cylinders are considered. The hydrostatic stress distributions decrease exponentially for all considered type of contacts. Considering the geometrical fabric, two specific directions of 45 with respect to the horizontal are exhibited when the whole set of contacts, the contacts between two dry particles, and the contacts between two cohesive particles are taken into account. It is noteworthy to conclude that the thermoelastic stress analysis experiment enables us to preliminarily reveal the hydrostatic stress network in cohesive granular media. This technique would open a chance to further investigate on cohesive granular materials. The numerical simulations provide the results, which are in good global agreement with the experimental results. Moreover, the property is in agreement with a part of the bimodal character of the material (as in non-cohesive granular material [START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]).

CONCLUSION AND PERSPECTIVES

In this dissertation, the macroscopic mechanical behavior of both "non-cohesive composite" and "cohesive non-composite" granular materials was investigated from a micromechanical analysis of some local variables through the granular texture. Infrared (IR) thermography was employed to reveal the hydrostatic stress networks in the materials under study, using the so-called thermoelastic stress analysis (TSA) technique.

It must be noted that it is the first application of TSA to granular materials. No such stress measurement by TSA technique is available in the literature (IR thermography was already used to study heat transfer, microrelief and rill morphology, and mechanical dissipation due to friction). The molecular dynamics (MD) method was also used in this study. This is a widely used numerical approach in the community of granular materials, enabling us to provide the contact forces between particles. The present study is based on the complementarity of these experimental and numerical approaches.

The present work first focused on "composite" granular materials, i.e. made of different constitutive materials, which is not usual in the literature. In practice, twodimensional (2D) analogue cohesionless composite granular samples were prepared from an assembly of cylinders of two constitutive materials: polyoxymethylene (POM) and high-density polyethylene (HDPE). The contact stiffness of POM was about four times greater than that of HDPE. Twelve configurations with difference in terms of ratio of diameter sizes and ratio of particle numbers were examined. Both TSA experiments and MD simulations were used for this purpose. A confined compressive loading was applied to the granular samples under quasi-static conditions. From an experimental point of view, the micromechanical analysis was performed through the hydrostatic stresses (contrary to the contact forces in simulations). The following results were obtained:

 In all tested configurations, the cylinder numbers transmitting the hydrostatic stresses greater than the mean value were less than 50% of the whole number of cylinders. It was also regardless of considered material types (all cylinders, only POM cylinders, and only HDPE cylinders).



The hydrostatic stresses were found to be distributed as an exponential decay, which was similar to the well-known distribution law found for the contact forces.



The angular distribution of contact directions for the whole set of contacts was quasi-homogeneous. In addition, the contacts between particles with the higher stiffness (POM-POM contacts) were oriented along the axis of compression. It can be said that the external applied loading was mainly carried by the POM-POM contacts, while the other contact networks, i.e.

between HDPE-HDPE cylinders and between POM-HDPE cylinders, just only made the granular system in equilibrium by supporting the contact network between POM-POM cylinders.

The same conclusions were obtained by the numerical simulations.

Next, the influence of the ratio of particle numbers as well as the ratio of contact stiffnesses was systematically studied by using only the MD simulation (still on 2D cohesionless composite granular materials). Coordination number gradually increased with the ratio of particle numbers, but it seems to be steady with the contact stiffness ratios. The change of both ratios did not affect the contact forces distribution:  About 40% of the whole contacts transmitted forces were higher than the average value.



An exponentially decreasing function was exhibited in the strong network.



These properties were also observed whatever the considered type of contact.

The contact stiffness ratio did not have an influence on the angular distribution of contact directions, contrary to the ratio of particle numbers which had the effect when high enough ratios are achieved:



The orientation of the whole set of contacts in the system was found to be quasi-isotropic, while a triangular contact network (as in monodisperse granular materials) was obtained when the ratio of particle number was high enough.



The triangular contact network was also evidenced in the contacts between lower-lower stiffness.



The contacts between higher-higher stiffness were oriented along the compression axis, while the contact networks between lower-lower stiffness and higher-lower stiffness exhibited an isotropic distribution.



It is noteworthy to say that the contacts between particles with the higher stiffness mainly carry the external applied loading, while the others just make the granular system in equilibrium.

In the last chapter, a preliminary study on cohesive non-composite granular media was presented, using TSA technique and MD simulation. A 2D polydisperse granular system was prepared. In practice, it consisted of an agglomerate of cohesive granular media surrounding with a collection of non-cohesive particles. A glue was used for creating cohesive bonding. A quasi-static confined compressive loading was then applied to the granular system. The experimental results shows that the percentage of cylinders, carrying the stresses greater than the average value, was approximately lower than 50% of the whole set of the cylinders. This property was also found when considering the contacts between particles "with cohesion", or only the contacts between particles "without cohesion". The distribution of hydrostatic stresses decayed exponentially for all considered type of contacts. Considering the angular distribution of contact directions, the arrangement of their structures was appeared on two particular directions of 45 with respect to the horizontal for all considered type of contacts. As a general comment, the numerical results were in good global agreement with the experimental results.

IR thermography under the framework of the TSA technique enabled us to measure the stress networks in the granular materials under study. This method is of some cautions in the test should be paid attention. First, the variation of thermal emissivity with the angle between the IR camera and the measurement surface, with the types of material, and with the colors of material should be concerned. Second, very small particle sizes (e.g. in sands, soils, or powders) compared to the pixel size in the measurement plane should be taken into account. This second caution could be solved by identifying mesoscopic hydrostatic stress fields.



An improvement of image processing could be concerned. First, the processing time could be greatly diminished by applying a shock load on the specimen, because of the large difference between the shock duration and the thermal diffusion characteristic time in material. Therefore, no complex image processing procedure would be required: the hydrostatic stress network would be determined by a simple multiplication of thermoelastic calibration factor A with the temperature changes after subjected to the shock load. Second, an algorithm to extract the contact forces from TSA experiments would be developed as done in photoelasticimetry.



For the case of cohesive granular media, it could be roughly said that the principle of thermoelastic stress analysis could not currently measure an attractive interaction between two particles in contact because the tensile stress levels are too low. In order to distinguish between repulsive and attractive interactions between contacting particles, a development of a new image processing would be necessary. The use of kinematic full-field measurements, i.e. digital image correlation (DIC), would be also another idea which could be used as a supplementary method for analysis.



Based on the elementary behaviors of composite granular materials found in this dissertation, the numerical simulations could be applied for the optimal design of the new granular systems in terms of constitutive material choice, proportion of particle number, and placement in the material, for instance.

However, these simulated granular systems in this study are restricted only in design under assumption of plane strain. Recalling the complexity of particle shape above, an assembly of composite granular materials with more complex shape will be also aimed at studying through the numerical simulations.

Finally, the author expects that the TSA approach would be powerful and widely used for the stress analysis in the experimental mechanics of granular materials in the near future. In addition, it would be interesting if it could be possible to extract the contact forces from the stress data provided by TSA technique.

  My beloved parents, brother, and girlfrienda ACKNOWLEDGEMENT First of all, I would like to express my sincere gratitude toward my advisors Assoc. Prof. Dr. Itthichai Preechawuttipong, Assoc. Prof. Dr. Xavier Balandraud, and Prof. Dr. Michel Grédiac for giving me the opportunity to work on this PhD and particularly for my very lucky chance to be a co-PhD student between Chiang Mai University (CMU) in Thailand and Blaise Pascal University (UBP) in France. Importantly, this thesis would have never been accomplished without their guidance, criticism, encouragement, assistance, and contribution. Next, I would like to appreciatively thank Dir. Rech. Dr. Farhang Radjaï and Asst. Prof. Dr. Supasit Rodkwan for their acceptance, time, and advice to act as the reviewers of this thesis. My acknowledgement also gives to the examination committee of this thesis Assoc. Prof. Dr. Theeraphong Wongratanaphisan, Asst. Prof. Dr. Somchai Pattana, and Dr. Wetchayan Rangsri for their time and advice. I am grateful to members of the Transfer Technology Centre (CTT) at French Institute for Advanced Mechanics (IFMA) in France Nicolas Blanchard, Authur Mota, Bernard Laracine, and Clément Weigel for the manufacture of some experimental device parts. Thanks also go to Hugues Perrin and Michel Drean for their instruction to use and setup the testing machines. I am thankful to members of computer vision research team (ComSee) at Pascal Institute in France, Assoc. Prof. Dr. Frédéric Chausse and Clément Beitone for their time, assistance, and effort to develop automatic image processing software. I would like to thank all members of granular materials and numerical simulations laboratory at Department of Mechanical Engineering, Chiang Mai University, Thailand for their time and assistance to prepare and organize a room for my PhD defense and also a place for my PhD defense after party. b I would like to gratefully acknowledge the Thailand Research Fund through the Royal Golden Jubilee Ph.D. Program (Grant No. PHD/0159/2552) and the French Embassy in Thailand for their financial supports during this research work. In particular, they enable me to obtain the worth experience in the foreign countries. Finally, I feel indebted and would like to express the biggest thanks to my parents Akom and Wipawan Jongchansitto, my younger brother Kunanon Jongchansitto, and my girlfriend Rapeepan Chaiwong for their unconditional understanding, tremendous patience, and ongoing support throughout my PhD study.

Figure 1 .1 examples of granular materials which can be found in our everyday life 9 Figure 1 . 2 13 Figure 1 . 3 14 Figure 1 . 4 15 Figure 1 . 5

 1912131314141515 Figure 1.1 examples of granular materials which can be found in our everyday life

Figure 1 . 7 Figure 1

 171 Figure 1.7 Probability density functions of normalized contact forces with the effect of: a) simulation methods, b) anisotropy, and c) particleshape[START_REF] Richefeu | Force transmission in dry and wet granular media[END_REF]

Figure 1 .

 1 Figure 1.11 Angular distributions of contact directions for three different 2D granular packings obtained from numerical simulations: a)

Figure 1 .Figure 1 .

 11 Figure 1.15 Normal forces with adhesion as a function of the virtual overlapdepth for the GA model in the static or quasi-static cases[START_REF] Preechawuttipong | Static states of cohesive granular media[END_REF]

Figure 1 .Figure 2 Figure 2 . 3 Figure 2 . 4 Figure 2 . 5 Figure 2 . 6 Figure 2 . 7 Figure 2 . 8 Figure 2 . 9 Figure 3 . 2 Figure 3 . 3 Figure 3 . 4 Figure 3 . 5 0 Figure 3 . 6 Figure 3 . 8 Figure 3 . 9 Figure 3 .Figure 3 .Figure 3 .Figure 4 . 1 Figure 4 . 2 Figure 4 . 3 Figure 4 . 4 Figure 4 . 5 Figure 4 . 6 Figure 5 . 1 Figure 5 . 2 Figure 5 . 3 Figure 5 . 4 Figure 5 . 5 Figure 5 . 6 Figure 5 . 7 Figure 6 . 1 Figure 6 Figure 6 . 3 Figure 6 . 4 132 Figure 6 . 5

 122324252627282932333435036383933341424344454651525354555657616636413265 Figure 1.23 Processing of Fourier transform at each pixel

Finally, chapter 6

 6 presents a preliminary study about the hydrostatic stress network in a cohesive granular media. The objective is actually to open prospects for further study. TSA experiments and MD simulations are also carried out. In practice, cohesive bonding is created using glue. Results obtained from the experiments and simulations are presented in the form of a map of normalized hydrostatic stress. These results are analyzed through two types of textural parameter: the distribution of hydrostatic stresses and contact directions.

cohesive (see chapter 3 to 5 )

 5 and cohesive (see chapter 6) granular materials are also presented;  Section 1.4 provides a background on infrared (IR) thermography and processing in the framework of the thermoelastic stress analysis (TSA) technique. It must be recalled (see general introduction) that TSA is a new experimental technique for the community of granular materials. This full-field, non-destructive, and non-contact technique is based on IR thermography to measure first the temperature variations at the specimen surface. In this section, two classical thermoelastic relationships are derived from different constitutive equations: classical theory of thermodynamics of elastic solids and Stefan-Boltzmann law. The fundamental of IR radiation

Figure 1 . 2

 12 Figure 1.2 Normalized shear stress as a function of cumulative shear strain under asimple shear test for different initial solid fractions[START_REF] Radjaï | Shear strength of granular materials[END_REF]

Figure 1 . 3

 13 Figure 1.3 Relationship between the coordination number and the cumulative shearstrain for pentagonal samples (S1) and circular samples (S2)[START_REF] Azéma | Force transmission in a packing of pentagonal particles[END_REF]

Figure 1 . 4

 14 Figure 1.4 Force networks of the static granular sample [118]. The black lines, which are connected between the centers of the disks, represent the amplitude of the contact forces. The line thickness is proportional to the force amplitude

Figure 1 . 5

 15 Figure 1.5 Contact forces measurement: a) stress pattern in a uniaxial compression of photoelastic disks [67], b) experimental setup of the carbon paper method [121], and c)

Figure 1 .

 1 6 presents a semi-logarithmic plot of probability distribution of contact forces obtained from uniaxial compressive simulations[START_REF] Radjaï | Force distributions in dense two-dimensional granular systems[END_REF]. On this graph, the forces F are normalized by the mean force <F>. Note that only normal contact forces are taken into account in Figure1.6.

Figure 1 . 6

 16 Figure 1.6 Probability distribution of normalized contact forces in semi-logarithmic scale [118]. Four different 2D samples in term of particle size distributions and particle numbers are simulated under a uniaxial compression

Figure 1 .

 1 7a displays the probability density functions (pdf) of normal forces in two isotropic samples consisting of spherical particles. Two distributions of contact forces simulated by molecular dynamics (MD) and contact dynamics (CD) methods have the same shape, i.e. exponential decreasing for the strong forces. The position of peak is not the same for the two distributions, but the exponent of the exponential decreasing are identical with   1.4.

Figure 1 . 7

 17 Figure 1.7 Probability density functions of normalized contact forces with the effect of: a) simulation methods, b) anisotropy,and c) particle shape[START_REF] Richefeu | Force transmission in dry and wet granular media[END_REF]

Figure 1 .

 1 Figure 1.10 Schematic for definition of contact directions θ

Figure 1 .

 1 Figure 1.11 Angular distributions of contact directions for three different 2D granular packings obtained from numerical simulations: a) deposition under gravitational force [140], b) packing under isotropic compression[START_REF] Liu | Microscopic study on stress-strain relation of granular materials[END_REF], and c) packing under biaxial compression[START_REF] Liu | Microscopic study on stress-strain relation of granular materials[END_REF]

Figure 1 .

 1 Figure 1.13 Sketch of molecular dynamics algorithm

Figure 1 .

 1 Figure 1.15 Normal forces with adhesion as a function of the virtual overlap depth for

Figure 1 .

 1 Figure 1.16 Tangential force F t as a function of sliding relative velocity v s under Coulomb's law of friction: a) exact form and b) regularized form[START_REF] Cambou | Micromechanics of Granular Materials[END_REF]

  of a Si detectors; Short wavelength infrared (SWIR) with 1-3 μm of the wavelength range is a sensitivity limit of PbS and InGaAs detectors; Middle wavelength infrared (MWIR) with the wavelength range between 3-6 μm is suited for InSb, PbSe, PtSi, and HgCdTe detectors which optimized for 3-5 μm atmospheric window; Long wavelength infrared (LWIR) with 6-15 μm of wavelength range is suitable of a HgCdTe detector which optimized for 8-14 μm atmospheric window; Very long wavelength infrared (VLWIR) is lied in the wavelength range between 15-1000 μm. Currently, most of used infrared cameras are mainly focused on the wavelength of two atmospheric

Figure 1 .

 1 Figure 1.19 Spectral hemispherical emissive power of black body in vacuum for several absolute temperatures as a function of the wavelength[START_REF] Astarita | Infrared thermography for thermo-fluiddynamics[END_REF]

Figure 1 .

 1 Figure 1.21 displays atmospheric transmission through 6000 ft of air as a function of wavelength. It indicates that water vapor (H 2 O) and carbon dioxide (CO 2 ) are significant components that absorb the transmitted radiation in the two atmospheric windows. Remark that the atmospheric

Figure 1 .

 1 Figure 1.21 Transmission of a 6000 ft horizontal thick atmospheric layer at sea level containing precipitate water of 17 mm[START_REF] Rogalski | Infrared devices and techniques[END_REF]

Figure 1 .

 1 22 presents a schematic diagram of the typical components of IR thermography systems. The thermal radiation emitted by a body enters the IR camera. This radiation is collected by the optical systems and sent to the IR detectors within the measured infrared band. An electric signal generated by the detector is then processed through image processing software packages on the computer. A map of temperature is therefore obtained.

Figure 1 .

 1 Figure 1.22 Ordinary components of infrared thermography systems

During a mechanical cyclic

  loading in the TSA experiments, the optical systems of the IR camera are focused the thermal radiation emitted by a tested specimen toward the focal plane arrays. In fact, the Cedip Jade III-MWIR infrared camera is a photon counting device. Depending on the amount of photons received by the FPA, this sends back into an electronic voltage by the readout integrated circuits (ROICs) in each pixel. The electronic voltage measured in each pixel is then expressed as a digital level (DL) by the analogue-to-digital converter. The experimental data is thus stored in the form of digital level information in a sequence of images versus time. The digital level data can be converted into the temperature data by using an experimental transfer function that obtained from the calibration process. It is worthy to note that the terms of detector voltage output S and thermoelastic signal S is an uncalibrated signal similar to the digital level DL. With Altair software, either digital level or temperature that measured by IR camera is provided in form of a video. In order to evaluate a temperature oscillation as a requirement of TSA, this temperature film is subtracted by its average temperature image that considers the film. The amplitude T of the temperature oscillation at the loading frequency is then determined by Fourier transform at each pixel as illustrated in Figure 1.23.

Figure 1 .

 1 Figure 1.23 Processing of Fourier transform at each pixel

Figure 2 Figure 2 . 1

 221 Figure 2.1 Contact stiffness measurement

Figure 2 . 2

 22 Figure 2.2 Experimental setup for contact stiffness testing

Figure

2. 3

 3 presents some results of contact stiffness for both POM and HDPE cylinders. It can be observed that the relationship between force and displacement is well-fitted by a linear regression.

Figure 2 . 3

 23 Figure 2.3 Relationship between force and displacement: a) POM cylinders with different diameters subjected to compressive loading of -20 kN and b) different compressive loading applied on HDPE samples

  Measurement: Two rectangular sheets of POM and HDPE materials were prepared with dimensions (widthheightthickness) 30.910.04.9 mm for POM and 10.46.05.0 mm for HDPE as shown in Figure 2.4. A uniaxial cyclic tensile loading with a frequency of 3 Hz and a force ratio of 0.1 was applied. In order to check the linearity, several maximum force amplitudes were tested: 200, 400, 600, and 800 N for POM and 80, 100, 150, and 200 N for HDPE. All tests were performed at the ambient temperature with a uniaxial MTS 15 kN testing machine.A Cedip Jade III-MWIR camera was used to capture the temperature fields on the front face of the sheet during each cyclic loading. As indicated above, the specimens were not painted in black in order to determine the value of the calibration factor A in this configuration.

Figure 2 .

 2 5 presents an experimental setup for POM and HDPE specimens.

Figure 2 . 4 Figure 2 . 5

 2425 Figure 2.4 Specimens used to determine the thermoelastic calibration factor A

Figure 2 . 6

 26 Figure 2.6 Determination of the amplitude of temperature oscillation at the loading frequency by Fourier transform. The inset shows zoomed image of the temperature amplitude around the loading frequency (here 3 Hz)

Figure 2 . 7

 27 Figure 2.7 Results of thermoelastic calibration factor testing: a) POM and b) HDPE

Figure 2 . 8

 28 Figure 2.8 Identification of the surface energy

Figure 2 . 9

 29 Figure 2.9 Example of specimen used in the testing of surface energy: two cylinders are glued along a generatrix

Figure 2 .

 2 Figure 2.10 Kind of glue: a) TOA adhesive latex, b) SADER colle bois, c) SADER fixer sans percer, and d) PATTEX ni clou ni vis

Figure 2 .

 2 Figure 2.12 shows a rupture point of specimens 2 glued by "PATTEX ni clou ni vis" based on a relationship between displacement and force.

3. 2 .Figure 3 .

 23 Figure 3.1b. The samples were then subjected to a confined vertical compression using a uniaxial Schenck ±250 kN testing machine as displayed in Figure 3.1a. The confined loading was provided by the lower hydraulic actuator of the testing machine. Twelve configurations were systematically tested. They differ in term of a ratio between POM and HDPE diameter and a ratio between POM and HDPE number as demonstrated in Table 3.1. In order to examine what the behaviors observed in non-composite granular materials also exhibit in the case of composite granular materials, both monodisperse and polydisperse configurations are taken into account.

Figure 3 .

 3 Figure 3.1 a) Picture of the experimental setup. Two reinforcement plates are added at the top the structure to ensure the rigidity of the frame. These plates prevent any bending of the two lateral parts of the frame under loading. Note that POM and HDPE cylinders are in black and white, respectively. b) Dimensions of the frame (in millimeter)

  1a) was here employed to capture the temperature fields on the front face of the cylinder during this cyclic loading. The integration time used for the measurements was set to 1500 s with the acquisition frequency of 148 Hz. The pixel size on the measurement plane was equal to 1.60 mm. Considering the thermal diffusivity D of POM and HDPE materials, thermal diffusion length / Df [START_REF] Offermann | Thermoelastic stress analysis under nonadiabatic conditions[END_REF] was about 0.1 mm for both materials, which was much lower than the cylinder diameter and the pixel size in the measurement plane. This criterion means that the adiabatic condition was achieved.

Figure 3 . 2

 32 Figure 3.2 Schematic of post-processing of TSA results to map of temperature changes: a) thermal image captured from a film, b) field of temperature changes, and c) applying Fourier analysis at each pixel

3 of chapter 2 ,Figure 3 . 3

 233 Figure 3.3 reveals a sample of hydrostatic stress network in a polydisperse granular medium. The bigger cylinders distinctly demonstrated specific stress distributions. It seems that some cylinders encounter boundary effects, likely owing to the pixel shifting which is located along the boundary between cylinders and voids. Concerning case of the small cylinders, the boundary effects should be carefully considered compared tothe spatial resolution of the measurement. As a consequence, the mean hydrostatic stress in each cylinder was practically measured at the center of each cylinder over 50% of all pixels inside cylinder area. In other words, the mean hydrostatic stress  hyd p in each cylinder was calculated by the ratio between sum of hydrostatic stress measured at the center of each cylinder over 50% of all pixels inside cylinder area and the number of pixels inside cylinder area.

Figure 3 . 3 (

 33 Figure 3.3 (Color) Illustration of hydrostatic stress field in polydisperse granular media

3. 2 . 3

 23 Experimental Normalized Hydrostatic Stress Networks Let us first consider the monodisperse samples which are relevant to the same diameter of samples. The normalized hydrostatic stress  norm field for sample #1, #2, and #3 are shown in Figure 3.4b. The ratio between number of POM cylinders and HDPE cylinders (N POM /N HDPE ) is distinct for each sample: 1.82 for sample #1, 1.06 for sample #2, and 0.52 for sample #3. Due to the nature of monodisperse granular materials, a triangular contact network was evidently exhibited in sample #1 and #2, which is corresponding to 30 of inclined stress paths with respect to the vertical direction. On the contrary, the field of  norm for sample #3 was separated into three zones by tilted light blue dash lines. Despite the fact that the cylinders are at random filled in the frame, a specific stress pattern was observed at the mesoscopic level for sample #3 owing to unintentionally created configuration before compression tests as illustrated in Figure 3.4a. This configuration composes of three zones, which is separated by two parallel inclined straight lines similar to grain boundaries.

Figure 3 . 4 (

 34 Figure 3.4 (Color) Experimental results obtained from monodisperse samples #1, #2, and #3 which are different in term of number of POM and HDPE cylinder ratios: a) optical images before loading and b) normalized hydrostatic stress networks

Figure 3 . 5 .

 35 Figure 3.5. All polydisperse samples clearly demonstrated more complex hydrostatic stress fields than monodisperse cases, even though both of them display an inhomogeneous stress distribution.

Figure 3 . 5 (

 35 Figure 3.5 (Color) Normalized hydrostatic stress fields for polydisperse configurations together with their optical pictures prior to loading: a) D POM /D HDPE = 1.2, b) D POM /D HDPE = 1.6, and c) D POM /D HDPE = 3.0

Figure 3 . 6

 36 Figure 3.6 Solid fraction for sample #11 determined by: a) varying size of representative elementary volume and b) randomly changing position of a specified size of representative elementary volume

Figure 3 . 7 (Figure 3 . 8 (

 3738 Figure 3.7 (Color) Normal plot of probability distribution function P of the normalized hydrostatic stresses in the strong network for all samples by considering: a) all cylinders, b) only POM cylinders, and c) only HDPE cylinders

2 .

 2 Figure 3.9 demonstrates a scheme of numerical sample preparations.

Figure 3 .

 3 Figure 3.10 (Color) Results obtained for sample #2, #8, and #12. a) Optical picture of initial configuration. b) Normalized hydrostatic stress networks  norm from TSA experiments. c) Representation of normal contact force networks from MD simulations. d) Field of  norm acquired from the simulated contact forces

Figure 3 .

 3 Figure 3.11 (Color) Experimental and numerical results of the angular distribution of contact directions for both mondisperse and polydisperse samples considered by: a) all contacts, b) only contacts between POM and POM cylinders, c) only contacts between HDPE and HDPE cylinders, and d) only contacts between POM and HDPE cylinders

Figure 4 . 1

 41 Figure 4.1 Simulation scheme for study the effect of the particle number: a) preparing state and b) compression state

Figure 4 .

 4 2 presents the beginning and final of compression state for the sample with N HDPE /N POM = 8. A big black circle represents POM particle, while HDPE particle is displayed by the small cyan circle. The normal contact forces are indicated by the solid red lines, which the line thickness was proportional to the magnitude of such normal force. The complex force networks are obviously exhibited as shown in Figure4.2b. It is interesting to note that the larger the particle size, the stronger the forces transmitted.

Figure 4 . 2

 42 Figure 4.2 Granular sample in the case of N HDPE /N POM = 8: a) beginning of compression state and b) final of compression state

Figure 4 . 3 (

 43 Figure 4.3 (Color) Influence of the N HDPE /N POM ratio on coordination number z

Figure 4 .

 4 Figure 4.4 demonstrates the influence of the N HDPE /N POM ratio on the distribution of percentage of forces in the strong network.In the analysis, it should be noted that the term of "whole contact networks" or "all contacts" refers to all types of contact are considered: POM-POM, HDPE-HDPE, and POM-HDPE. The percentage of the contact forces greater than the mean

Figure 4 . 4 (

 44 Figure 4.4 (Color) Influence of the particle number on the distribution of percentages of forces greater than the mean force. The solid line refers to the whole contact network, while the dash line refers to their contact network

4. 5 .

 5 These probability distributions P are the number of contacts which carry the forces greater than the mean value divided by the total number of contacts. The forces are normalized by their mean value in the whole normal contact forces. It is interesting to see that the POM-POM contacts are able to carry forces which are higher than those carried by POM-HDPE and HDPE-HDPE contacts. It can be said that contacts between POM and POM particles are able to transmit the forces greater than the other contacts, due to the chance of contact occurrence between POM and any particles as mentioned section 4.3.1. This result is identical for all tested N HDPE /N POM ratios. All the types of contacts follow an exponential decay as commonly found in granular materials made of only one type of constitutive material[START_REF] Radjaï | Force distributions in dense two-dimensional granular systems[END_REF][START_REF] Mueth | Force distribution in a granular medium[END_REF][START_REF] Løvoll | Force measurements on static granular materials[END_REF]. It is a general characteristic for granular media, which is independent of the N HDPE /N POM ratio or the type of contacts.

Figure 4 . 5 (

 45 Figure 4.5 (Color) Semi-logarithmic plots of probability distributions of normalized normal forces F/<F>. The black square represents the result over the whole contacts. Red, blue, and green colors correspond to the results for each type of contact: POM-POM, HDPE-HDPE, and POM-HDPE, respectively

  Figure 4.6. Their directions are divided into 18 angular sectors. The data is plotted in the upper part of polar diagram as the angular range  [0,180 o ], while the range of [0,-180 o ] is the symmetrical range as in the upper part.

Figure 4 . 6 (

 46 Figure 4.6 (Color) Influence of the N HDPE /N POM on the direction of contacts of: a) the whole, b) POM-POM, c) HDPE-HDPE, and d) POM-HDPE

Figure 5 . 1 ,

 51 the samples were subjected to a compressive loading of 100 N on the upper wall, while the other walls were fixed during this stage. The gravity effect was not considered. Considering the contact forces, they consist of normal and tangential forces excluding the physicchemical reaction. The normal and tangential forces in each contact are calculated by use of eq. (1.23) and eq. (1.27), respectively. The simulations were completely performed when the systems reached quasi-static equilibrium conditions.

Figure 5 . 1

 51 Figure 5.1 Scheme of granular packing preparation: a) beginning of simulation and b) end of simulation

Figure 5 . 2

 52 Figure 5.2 Granular sample in the case of k A /k B = 4 at the equilibrium state: a) display without contact forces and b) display with contact forces

Figure 5 . 3 (

 53 Figure 5.3 (Color) Influence of the k A /k B ratio on coordination number z

Figure 5 .

 5 Figure 5.4 displays the influence of the k A /k B ratio on the force distribution percentages that are higher than the mean force. For the whole contact network, the percentage of contacts, which carries a force greater than the mean force, tends to be steady and transmits forces (see solid black line) at the level of approximately 40%. This is also observed in the contact networks between materials A-A (Higher stiffness), B-B, and A-B with respect to the number of their own contacts of each type (dash red, blue, and green line respectively). This property is consistent with a generic property found in granular media made of one type of constitutive material only. It seems that the property is not dependent on the ratio of contact stiffness or the type of contacts. In fact, the 40% component of the whole contact network consists of three contact networks: contacts between materials A-A, B-B and A-B. Most of contact networks occur between materials A-B (see the solid green line), while the networks between A-A (solid red line) and B-B (solid blue line) exhibit the same proportion about one-fourth of the contacts transmitting force higher than the mean force.

Figure 5 . 4 (

 54 Figure 5.4 (Color) Influence of contact stiffness on the distribution of percentages of forces greater than the mean force <F>. The solid line refers to the percentages separated by the whole contact network, while the dash line refers to the percentages separated by their contact network

Figure 5 . 5 (

 55 Figure 5.5 (Color) Semi-logarithmic plots of probability distributions of normalized normal forces F/<F>. The black square represents the result over all types of contacts. Red, blue, and green colors correspond to the results for each type of contact: material A-A, B-B, and A-B, respectively

  next section. In order to clearly characterize and describe the geometrical fabric, the polar diagram of the number of contacts N c in the normal directions is here presented in Figure 5.6. Their directions are divided into 18 angular sectors. The data is plotted in the upper part of polar diagram as the angular range  [0,180 o ], while the bottom part of the graph corresponding to the range of [0,-180 o ] is the symmetrical range as in the upper part.

Figure 5 . 6 (

 56 Figure 5.6 (Color) Influence of the k A /k B on the angular distribution of contact directions of: a) the all types of contacts, b) A-A, c) B-B, and d) A-B

  normal forces in each contact is necessary to consider. The number of contacts N c in the strong (see solid line) and weak (dash line) networks are then separately plotted by the polar diagram as illustrated in Figure 5.7. Their directions are split into 18 angular sectors. The upper part of the polar diagram presents the data in the angular range  [0,180 o ], while the range of [0,-180 o ] in the bottom part of this graph is the symmetrical range as in the upper part.

Figure 5 . 7 (

 57 Figure 5.7 (Color) Influence of the k A /k B on the angular distribution of contact forces of: a) the all types of contacts, b) A-A, c) B-B, and d) A-B. The solid line refers to the strong contact network, while the dash line refers to the weak contact network

Figure 6 . 1 Figure 6 . 2

 6162 Figure 6.1 Granular system made of cohesive and non-cohesive grains. The cohesive agglomerate is placed inside the red solid boundary

Figure 6 . 3 (

 63 Figure 6.3 (Color) Results for the case of cohesive granular media: a) Normalized hydrostatic stress networks  norm from the TSA experiment, b) Contact force networks obtained by the MD simulation, and c) fields of  norm derived by the simulated contact forces. The zone contained by the light blue zone represents the cohesive agglomerate

Figure 6 . 4 (

 64 Figure 6.4 (Color) Plots of probability distribution of normalized hydrostatic stresses in a) normal and b) semi-logarithmic scale. The dot light blue and dash light orange lines represent the linear regression of all data obtained from the experiment and the simulation, respectively

  arrangement of the particles. It also indicates a development of an anisotropic structure of the granular packing. The detection and occurrence of contacts for both experiment and simulation are based on the same procedure as that used in chapter 3.

Figure 6 . 5 (

 65 Figure 6.5 (Color) Polar diagram of the statistical distribution of normal contact directions for experimental and numerical results. Considerations of a) the whole contacts, b) the contacts between two dry particles, c) the contacts between two cohesive particles, and d) the contacts between dry and cohesive particles

Table 2 .

 2 2 Contact stiffness of HDPE for each configuration (in kN/mm)

Table 2 .

 2 3 Separation force in unit of N and surface energy in unit of kN/mm

Table 3 .

 3 1 All configuration of tested composite granular samples

Table 3 .

 3 2 Percentage P of cylinders number exhibiting hydrostatic stresses greater than the mean value for all experiments. The results are considered the whole cylinders, only POM cylinders, and HDPE

	cylinders, respectively

Table 6 .

 6 1 Percentage of cylinders carrying the stresses greater than the mean value obtained over the whole cylinders, only non-cohesive cylinders, and only cohesive cylinders f

  The full-field measurements are powerful and useful tools compared to classical measurement techniques which depends on the use of thermocouples, strain gauges, or displacement transducers. These classical methods enable us to characterize the material response only at a single value of the local scale, while a field of local characterization that provides a global information about the material under study can be performed by using full-field approaches. As far as experiments on granular media are concerned, the use of full-field measurement techniques in community of granular materials is in fact not completely new. In the late

1950s, photoelasticimetry was already applied to granular media by Dantu, as mentioned above. A number of literatures on granular materials studied by full-field measurement techniques have increased since that time. Most of them focus on the kinematic information in the granular medium: in-plane measurement techniques such as digital image correlation

Table 1 .

 1 1.

Table 1 .

 1 [START_REF] Jaeger | Physics of the granular state[END_REF] Values of physical constants of interest in thermal radiation theory[START_REF] Maldague | Theory and practice of infrared technology for nondestructive testing[END_REF]

	Quantity	Symbol	Value	Units
	Speed of light	c	2.997925×10 8	m/s
	Planck's constant	h	6.626076×10 -34	J•s
	Boltzmann's constant	k B	1.380658×10 -23	J/K
	Stefan-Boltzmann's constant	B	5.67051×10 -8 W/m 2 •K 4

Table 1 .

 1 2.

Table 1 .

 1 2 Technical specification of Cedip Jade III-MWIR infrared system[START_REF]DC001U-B Jade User Manual[END_REF]

	Detector materials	InSb
	Wavelength band	3.6 -5.1 μm
	Number of pixels	320 × 240 pixels
	Pixel size	30 × 30 μm
	Frame rate	5 -170 Hz full frame
	Noise equivalent difference temperature (NETD)	< 25 mK at 25C
	Cooling type	Integral stirling cooler
	Cooling time	< 10 minutes at 25C
	Integration time	10 -20000 μs programmable, in 1 μs step

Table 2 .

 2

		D = 6.2 mm D = 16.4 mm D = 30.9 mm
	F = -6 kN	24.97	24.97	22.93
	F = -8 kN	29.94	25.52	25.07
	F = -10 kN	30.74	24.91	27.21
	F = -12 kN	34.55	26.62	28.20
	F = -15 kN	30.35	26.62	29.71
	F = -20 kN	25.07	25.14	30.89
	Average in each size	29.27	25.45	27.33
	Total average		27.35

1 Contact stiffness of POM for each configuration (in kN/mm)

Table 2 .

 2 1 presents the results of POM's stiffness for each sample, while the results of HDPE's stiffness are given in Table2.2. For the POM cylinders, it can be noted that the contact stiffness is not very sensible to the diameter for the three values of diameter considered. A mean value will be used in the following of the study: 27.35 kN/mm. For the HDPE cylinders, it can be seen that the values are quite variable as a function of the test.

However, a mean value is here also considered: 6.43 kN/mm. It can be noted that this value is about 4 times lower than the value found for POM cylinders.

Table 2

 2

	.2 Contact stiffness of HDPE for each configuration (in kN/mm)
	F = -3 kN (1 st ) F = -3 kN (2 nd ) F = -5 kN (1 st ) F = -5 kN (2 nd ) F = -8 kN (1 st ) F = -8 kN (2 nd )	D = 10.5 mm 5.18 5.20 6.96 6.35 9.76 5.05
	Average	6.43

Table 2 .

 2 3 Separation force in unit of N and surface energy in unit of kN/mm

		TOA adhesive	SADER colle	SADER fixer	PATTEX ni
		latex	bois	sans percer	clou ni vis
	Specimens 1	6.42	36.00	75.91	96.75
	Specimens 2	8.38	33.36	49.88	119.04
	Specimens 3	9.28	31.26	47.84	112.94
	Specimens 4	6.14	27.62	63.33	164.76
	Specimens 5	10.49	19.76	77.70	144.19
	Separation
	force in	8.14	29.60	62.93	127.54
	average, N
	Surface energy
	from Eq.(1.40),	0.268	0.511	0.745	1.060
	kN/mm

Figure 2.12 Plotting of force versus displacement for specimen 2 glued by PATTEX ni clou ni vis

Table 3 .

 3 1 All configuration of tested composite granular samples

		Diameter ratio Number of cylinders ratio
		POM : HDPE	POM : HDPE
	#1	1.0	799 : 436 (1 : 0.55)
	#2		638 : 597 (1 : 0.94)
	#3		435 : 807 (1 : 1.86)
	#4	1.2	614 : 320 (1 : 0.52)
	#5		477 : 509 (1 : 1.07)
	#6		345 : 705 (1 : 2.04)
	#7	1.6	365 : 295 (1 : 0.81)
	#8		334 : 371 (1 : 1.12)
	#9		239 : 608 (1 : 2.54)
	#10	3.0	118 : 137 (1 : 1.16)
	#11		106 : 289 (1 : 2.73)
	#12		86 : 466 (1 : 5.42)

Table 3 .

 3 2 Percentage P 0 of cylinders number exhibiting hydrostatic stresses greater than the mean value for all experiments. The results are considered the whole cylinders, only POM cylinders, and HDPE cylinders, respectively

		Diameter	Number of	P 0 for	P 0 for	P 0 for
		ratio	Cylinders ratio	all	only POM only HDPE
		POM : HDPE	POM : HDPE	cylinders cylinders	cylinders
	#1		799 : 436 (1 : 0.55)	44.5%	42.8%	45.9%
	#2	1.0	638 : 597 (1 : 0.94)	47.4%	50.2%	50.9%
	#3		435 : 807 (1 : 1.86)	39.4%	40.9%	43.2%
	#4		614 : 320 (1 : 0.52)	36.2%	40.0%	32.7%
	#5	1.2	477 : 509 (1 : 1.07)	41.3%	43.6%	39.0%
	#6		345 : 705 (1 : 2.04)	42.1%	47.4%	40.2%
	#7		365 : 295 (1 : 0.81)	39.0%	38.4%	40.3%
	#8	1.6	334 : 371 (1 : 1.12)	46.3%	47.2%	45.1%
	#9		239 : 608 (1 : 2.54)	47.0%	50.9%	46.3%
	#10		118 : 137 (1 : 1.16)	28.2%	46.2%	39.0%
	#11	3.0	106 : 289 (1 : 2.73)	33.5%	52.2%	36.6%
	#12		86 : 466 (1 : 5.42)	45.6%	50.9%	44.3%

Table 6 .

 6 1 Percentage of cylinders carrying the stresses greater than the mean value obtained over the whole cylinders, only non-cohesive cylinders, and only cohesive

	cylinders
	Considering zones	Percentage of the stresses higher than the own mean
		Experiment	Simulation
	All	40.7%	48.0%
	Without adhesion	41.1%	48.0%
	With adhesion	50.1%	49.4%
	6.3.3 Angular Distribution of Contact Directions

CHAPTER 2 PRELIMINARY EXPERIMENTS

Introduction

The aim of the study is to analyze the mechanical response of two-dimensional (2D) analogue granular assembly in two cases: cohesionless composite and cohesive non-composite. The non-cohesive composite granular materials are constituted from two different polymeric materials: polyoxymethylene (POM) and high-density polyethylene (HDPE), while glue is employed to create the cohesive bonding.

The purpose of this chapter is to perform preliminary experiments on cylinders in order to determine some elementary properties of materials, namely contact stiffness (k) of POM and HDPE, thermoelastic calibration factor (A) of POM and HDPE, and surface energy () associated to the glue. These data will be used later in the study.

In section 2.2, Preliminary experiments are conducted to evaluate the value of contact stiffness of POM and HDPE cylinders. The thermoelastic calibration factor of the two constitutive materials is determined in section 2.3. The surface energy in the case of the use of glue is determined in section 2.4.

Contact Stiffness (k)

Introduction

The contact forces between two dry grains can be decomposed into two components: normal and tangential forces. In MD simulation, these interaction forces are simply modeled as a function of the virtual overlap.

The repulsive term in Eq. (1.23) is represented by a linear model, contrary to the Hertz's theory which is a non-linear model [START_REF] Brilliantov | Model for collisions in granular gases[END_REF]. According to the repulsive component of Eq. (1.23), the contact stiffness in this study is POM cylinders, and only HDPE cylinders. From these results, it can be seen that the experiments and the simulations are in good agreement.

Conclusion

In this chapter, the analysis of two-dimensional composite granular materials without cohesive interactions was performed through the granular texture. Both experimental and numerical approaches were used: thermoelastic stress analysis for the CHAPTER 4

ADDITIONAL STUDY: EFFECT OF THE PARTICLE NUMBER

Introduction

The objective of this chapter is to investigate the influence of the particle number on the granular texture of two-dimensional non-cohesive composite granular media by means of the molecular dynamics (MD) method. From the previous chapter, the composite granular materials with a difference in term of ratio of particle number and ratio of diameter size were systematically studied by thermoelastic stress analysis experiments based on infrared thermography. However, the experiments cannot directly provide the contact forces, contrary to numerical simulations. It can be noted that the macroscopic mechanical behavior of granular materials are governed by contact forces.

Moreover, there are many limitations of experiments such as particle number, particle size, measurement methods, testing machine, preparation, etc. In order to extend much more information and to eliminate the constraints from our experiments, molecular dynamics simulations are carried out in this chapter. One interesting parameter is the particle number because of a lack of data from our experiments: there are only three different ratio of the particle number in such diameter ratio. The diameter size ratio of 3 is aimed to study, owing to most a lack of statistical information compared with the others. In order to achieve much more statistic information, a number of particles in the system will be increased. The chapter is organized as follows. Sample preparations and simulation procedure are first presented. Numerical results are then analyzed by the granular texture.

Numerical Preparation

Ten two-dimensional composite granular samples were prepared from two types of constitutive materials: polyoxymethylene (POM) and high-density polyethylene (HDPE). The granular samples were composed of about 4000 particles which were CHAPTER 5

ADDITIONAL STUDY: EFFECT OF THE CONTACT STIFFNESS

Introduction

This chapter aims to investigate the influence of the contact stiffness on the textural properties of two-dimensional cohesionless composite granular materials by means of molecular dynamics simulations. It should be emphasized again that this research principally deals with the non-cohesive composite granular media in twodimensional. This purpose is accomplished by using the full-field experimental technique named infrared thermography. Recalling the weakness of experiments as mentioned in the previous chapter, they are not able to provide the contact forces as well as possess many limitations. In contrast, performing numerical simulations constitutes a powerful approach to study the mechanical behavior of granular media at the contact scale, and to dispose the limitation of experiments. One noteworthy component for modeling of contact forces is contact stiffness. The contact stiffness is absolutely dependent on such types of material. In addition, only one ratio of the contact stiffness has been used in our experiments. That is why the effect of the contact stiffness is interesting to study in this chapter. Organization of this chapter can be presented hereinafter. The first section describes numerical setup, loading conditions, and simulation steps. The second section provides the results in term of the granular texture.

Numerical Preparation

Ten configurations of two-dimensional composite granular media were composed of about 4000 particles. The particles were randomly placed inside a square box consisting of four rigid plane walls. The particle radiuses were uniformly distributed between 4 and 8 mm with average radius equal to 6 mm. The granular samples were prepared from two types of constitutive materials: A and B. The contact stiffness ratio of both materials (k A /k B ) was progressively changed from 1 to 10, while course complementary to other full-field measurement methods and also numerical simulations. Relied on the results shown in this dissertation, it is worth to mention that the TSA method is novel and promising. It is expected to be valuable experimental technique for the research community of granular materials in the future. Several prospects could be proposed, as discussed in the following paragraphs.



Aside from granular materials, the TSA approach has a good potential to perform experiments with any types of classical engineering constitutive material. A reason is that the applicability of this technique depends on the value of A and on the camera measurement resolution, which is known as

Noise Equivalent Temperature Difference (NETD) in the supplier datasheets. The Fourier transform used to estimate the temperature amplitude T of the thermoelastic signal oscillation can be seen as a filtering process. This leads to a significant improvement of the measurement resolution. In practice, the measurement resolution of T is equal to 1 mK whereas the NETD of the current camera is equal to 20mK.

The order of magnitude of A is dependent of the material properties (the density, the specific heat capacity, and the coefficient of thermal expansion) and the ambient temperature T 0 . Let us now consider two cases, different from polymeric materials as used in this study. For granitoid rocks and aluminium alloys at ambient temperature T 0  300 K, we obtain the following values: A  1 mK/MPa and A  3 mK/MPa respectively. Hence, a stress of 100 MPa for instance leads to a temperature change of 100 mK and 300 mK respectively for granitoid rocks and aluminium alloys. These temperature changes are much higher than the measurement resolution of T. Note that the stress level depends on the considered application and these numerical values only provide orders of magnitude.



A complexity of particle shape could be considered. An assembly of analogue granular materials with more complex shape (i.e. triangle, square, and polygon) would be studied next in the both composite and noncomposite cases. In addition, it is of interest to perform these experiments with real geotechnical materials such as rocks, sands, soils, etc. Anyway,

LIST OF PUBLICATIONS