Mehdi Bouzid

Adrien Izzet

Martin Trulsson

Eric Clément

Philippe Claudin

Bruno Andreotti

Non-local rheology in dense granular flows ⋆ Revisiting the concept of fluidity

Granular media, such as sand, cereals and gravels are ubiquitous in our world. Each of their constitutive particle is solid. However, considered as a whole, they exhibit surprising behaviors: sometimes in a jammed "solid" state, they can under certain conditions flow like a liquid, or even exhibit a gaseous regime. The characterization of this transition from one state to another is still an active field of research.

This thesis aims to characterize granular media in each of these two dense states, "solid" and "liquid".

The first part is dedicated to the characterization of the flow of a granular medium. Following preliminary studies by our team, we present a model taking into account the collective aspect of granular flows, particularly the cooperative behavior of regions above and below this jamming transition. This topic is introduced by a review paper in which we discuss the concept of fluidity in granular matter. We then present an experiment of avalanche flow in a channel. We probe this model on our experimental results and exhibit the non-local contributions on this threedimensional, frictional system. We then focus our interest on the definition of the boundary condition in the vicinity of the free surface of such an avalanche flow. We proceed by the mean of two-dimensional numerical simulation of rigid grains flowing down an incline plane. This numerical chapter is divided into three sections: we present the principles of the code and we use an original custom numerical set-up in which we exhibit and quantify these non-local effects, and then we focus our interest on the characterization of the flow at the free surface of the avalanche.

The second part of this thesis investigates the mechanical behavior of the "solid" phase of granular media, in particular in the vicinity of the unjamming of the packing.

We present the acoustic anomaly which characterizes the change of phase in granular media. We then present the experiment we built to measure elastic moduli at very low confining pressures in order to get closer to the unjamming of the medium. We propose a theoretical model in order to explain the constitutive dependence of the elastic modulus on the confining pressure and thus on the sound celerity in such granular packings.

naud Delaunay et Nicolas Taberlet, merci d'avoir accepté de rapporter cette thèse.

Antoine Séguin aussi, merci d'avoir examiné ce travail avec tant d'intérêt et cette complicité de physicien-mécanicien que j'ai appréciée. Enfin, "last but not least", merci Pierre-Yves Lagrée d'avoir présider ce jury. Merci de ton intérêt à mon égard et pour ces discussions de sciences, sur le non-local notamment. Je tiens à tous vous remercier de l'attention que vous avez porté à mon travail, ainsi que pour vos questions, commentaires et discussions durant et après la soutenance. J'ai énormément apprécié échanger avec vous.

Mes remerciements vont ensuite, bien évidemment, à mes directeurs de thèse.

Eric et Bruno, merci d'avoir accepté de m'encadrer ensemble pour ce travail. L'expérience a été pour moi riche tant personnellement que sur le plan scientifique. J'ai eu la chance de profiter de vos conseils formateurs et de vos expertises respectives.

Bruno, merci de m'avoir transmis ton souci de la méthode et de la rigueur en science. J'ai apprécié recevoir cet enseignement de ta part, ainsi que les discussions scientifiques enrichissantes que nous avons pu avoir. Je suis fier d'avoir pu bénéficier de cet enseignement de "théoricien de manip" que tu as su me transmettre. Merci enfin de m'avoir permis de faire des simulations quand l'étude expérimentale sur le canal a posé de nouvelles questions auxquelles il fallait répondre. Merci du recul que tu as su me faire prendre sur bien des choses.

Eric, je te remercie de m'avoir permis de faire des manips qui me plaisaient.

Le canal, pour la rhéologie non locale avec Bruno, mais aussi les expériences en apesanteur pour l'acoustique et le jamming. J'ai eu la chance de travailler à tes côtés et de participer à cette aventure qui je l'espère se poursuivra, peut-être même plus haut encore que les vols paraboliques. J'ai apprécié mettre en place ce projet et travailler avec la "team granu". J'ai beaucoup apprécié ton optimisme, notamment dans les moments de doutes, ainsi l'effervescence dans laquelle nous avons échangé pendant ces trois ans et quelques. Merci de tes conseils et de m'avoir un peu transmis de ton oeil "Origin" pour "regarder les data dans les yeux".

Je voudrais aussi remercier Etienne Guyon. Etienne, merci pour ton intérêt à mon égard, ta bienveillance et le temps que tu as pris plusieurs fois pour parler de science avec moi. Merci d'être si accessible et disponible malgré tes multiples activités et projets. Merci pour ton goût de la science que tu sais partager et ta curiosité. Merci de ta confiance qui m'est chère.

Je voudrais aussi remercier Philippe Petitjeans, directeur du laboratoire PMMH : merci pour ton soutien et ta bienveillance constante ; Damien Vandembroucq pour avoir été mon parrain de thèse et le temps que tu as pris pour discuter de ma thèse autour d'un café ou d'un chocolat chaud. Laurette, merci de ta relecture de mon intro, ta disponibilité et ton aide sur la topologie quand j'ai eu besoin de quelques rappels pour mes enseignements. J'ai eu la chance d'effectuer ma mission d'enseignement au sein de l'UFR de mécanique de l'UPMC et de travailler avec notamment Anne Mongruel et Hélène Dumontet. Anne, merci ta disponibilité et de la liberté que tu as su me laisser dans les TD que je faisais pour ton cours. J'ai beaucoup apprécié travailler à tes côtés et échanger avec toi sur l'enseignement. Hélène, merci pour ton soutien ainsi que ton attention à me confier des enseignements que j'ai beaucoup appréciés dispenser. Merci aussi à tous les autres membres de l'UFR ainsi que de l'Institut Jean Le Rond ∂'Alembert avec qui j'ai pu échanger lors de mes quelques passages.

Je tiens ensuite à remercier Adeline Pons, qui était en postdoc la première année de mon doctorat et avec qui j'ai travaillé sur l'expérience d'acoustique en ZéroG.

Adeline, merci pour ces moments au réfectoire, en salle R1 puis dans les bureaux devant les ordis afin de regarder les données et calibrer la manip notamment, ainsi que les acrobaties et grimaces pendant les vols. Merci pour tes conseils et tous tes enseignements sur Matlab et ta capacité à faire de la bonne science dans une ambiance joyeuse et studieuse à la fois. Merci de tes encouragements constants et pour cette amitié qui perdure malgré la distance.

Mehdi, je tiens ici à te dire un grand merci pour tes conseils. Tes encouragements et ta hargne pour avancer malgré les obstacles m'ont permis de garder la tête haute dans les moments difficiles, surtout vers la fin. J'ai aussi énormément apprécié nos échanges scientifiques, sur la rhéologie notamment. Merci d'avoir su rester proche malgré la distance et de m'avoir donné du courage tout du long. Merci pour ton amitié, qui m'est chère elle aussi. Martin Trulsson aussi, pour tes conseils ainsi que les parties de grimpes durant lesquelles il nous est aussi arrivé de parler rhéologie. J'espère qu'il y en aura d'autres, accompagnées de science bien évidemment.

J'ai passé beaucoup de temps en salle R4, notre salle de manip., avec les autres étudiants de Bruno. Avec Filippo quand il finissait sa thèse, Pan Jia et Hugo Perrin.

Les gars, merci pour l'ambiance studieuse mais décontractée pendant toutes ces heures où nous étions chacun sur nos manip. J'ai beaucoup apprécié ces moments ainsi que nos discussions scientifiques et autres sujets de sociétés que nous avons pu évoquer autour d'une tasse de thé chinois ou de café.

La science est une chose, mais elle ne pourrait être sans le soutien administratif que nous apportent au quotidien Frédérique Auger et Claudette Barez. Outre les remerciements que je veux vous formuler à cet égard, je tiens à remercier mes deux "tantes" de thèse que vous êtes devenus pendant ces années. Je voudrais vous dire combien votre affection et nos échanges m'ont permis de me sentir encore un peu plus si c'est possible dans cette famille qu'est le PMMH.

Ce dernier est un laboratoire avec une forte dominante expérimentale, et je tiens par conséquent à remercier le pôle technique et d'ingénierie sans qui la majeur partie des manip ne pourraient voir le jour: Olivier, Tahar, Amaury, merci de vos conseils pour l'usinage, ainsi que les rush de fabrication que l'on vous a confiés quand il fallait que la manip CNES soit prête à temps pour embarquer. Xavier, merci pour ton aide sur la CAO du rack, le dossier CNES à Novespace ainsi que pour les bons moments que nous avons passés sur la manip au sol, et dans l'avion. Thierry, merci pour ton aide et ton soutien technique et moral quand j'avais la tête dans la "cabane à ours" pour régler la cellule.

Merci à toute l'équipe de Novespace et du CNES (le CADMOS notamment) qui nous ont permis de participer à cette expérience inoubliable et qui ont rendu cette aventure technique et scientifique possible.

Je remercie aussi le Service Informatique de l'ESPCI (Emmanuel, Jean-Jacques, Camille. . . ) d'avoir pris le temps de répondre à mes nombreuses questions et requêtes sur le réseau. Merci pour votre efficacité et surtout votre bonne humeur. Je vais continuer à suivre cette affaire de HTCPCP-TEA ! Je voudrais aussi remercier tout particulièrement les bons potes du vendredi soir alias la team "balle neuve" (Matthew avec qui j'ai aussi eu le plaisir de parler science, Ali, Lukasz, Marine, Jessica, Yassine, Vince, Eliott, Jean, Benoit, Alex. . . ) pour les moments "bières" autour du billard ou de la table de pingpong. Mention spéciale aussi à Anthony-Darith que j'ai le plaisir de co-encadrer pendant son stage avec Eric.

Florence, merci de ton soutien lors de ma rédac, quand je me suis isolé dans ton bureau. Tu es une "super nana" comme ils le disent dans la pub. . . Je tiens aussi à remercier l'équipe des Quéré, qui a toujours la pêche (et souvent un jeu de coinche en mains), et avec qui j'ai eu plaisir à échanger. Je voudrais notamment te remercier chaleureusement Martin, parce que si j'ai pu m'envoler vers d'autres horizons, c'est beaucoup grâce à toi. Merci de ta confiance et de ta bienveillance. Merci au groupe des MécaWet d'être la team de scientifiques trop cools que vous êtes et pour votre disponibilité à répondre à mes questions sur l'élasticité.

Merci aux enseignants que j'ai eu, de m'avoir fait confiance et de m'avoir inspiré pour continuer.

Je voudrais remercier Ruben, qui m'a relu pour l'anglais, sur quelques parties de cette thèse. J'espère que nous trouverons le moyen garder contact même malgré ton tour du monde d'après thèse.

Justine, nous n'avons pas eu l'occasion de travailler beaucoup ensemble mais j'ai adoré tous nos échanges. Merci de ta sensibilité et de ta délicate empathie. Merci de tes conseils et de ton amitié.

Charles, merci de m'avoir appris à parler forain et de m'avoir présenté J. Cornella et tout un bout de la toile que je devrais peut-être vouloir oublier. Merci d'avoir été le co-bureau que tu as été et de m'avoir soutenu au quotidien.

Thibault, non je ne peux pas t'oublier. Merci pour ces magistraux moments de musique et d'amitié, avec notre duo des Pullers. Merci d'être l'ami que tu es devenu.

Sylvain, tu as pris ton poste de chargé de recherche en même temps que je débutais ma thèse. J'ai été le Padawan et tu a été l'Obi Wan. Merci pour ta bienveillance perpétuelle, tes conseils de grand frère autour d'un café ou devant le tableau pour parler science ou politique. Merci de m'avoir laissé utiliser le cluster et de m'avoir tant appris, en code mais aussi sur la vie.

Granular media, such as sand, cereals and gravels are ubiquitous in our world. Each of their constitutive particle is solid. However, considered as a whole, they exhibit surprising behaviors: sometimes in a jammed "solid" state, they can under certain conditions flow like a liquid, or even exhibit a gaseous regime. The characterization of this transition from one state to another is still an active field of research.

This thesis aims to characterize granular media in each of these two dense states, "solid" and "liquid".

The first part is dedicated to the characterization of the flow of a granular medium. Following preliminary studies by our team, we present a model taking into account the collective aspect of granular flows, particularly the cooperative behavior of regions above and below this jamming transition. This topic is introduced by a review paper in which we discuss the concept of fluidity in granular matter. We then present an experiment of avalanche flow in a channel. We probe this model on our experimental results and exhibit the non-local contributions on this threedimensional, frictional system. We then focus our interest on the definition of the boundary condition in the vicinity of the free surface of such an avalanche flow. We proceed by the mean of two-dimensional numerical simulation of rigid grains flowing down an incline plane. This numerical chapter is divided into three sections: we present the principles of the code and we use an original custom numerical set-up in which we exhibit and quantify these non-local effects, and then we focus our interest on the characterization of the flow at the free surface of the avalanche.

The second part of this thesis investigates the mechanical behavior of the "solid" phase of granular media, in particular in the vicinity of the unjamming of the packing.

We present the acoustic anomaly which characterizes the change of phase in granular media. We then present the experiment we built to measure elastic moduli at very low confining pressures in order to get closer to the unjamming of the medium. We then propose a theoretical model in order to explain the constitutive dependence of the elastic modulus on the confining pressure and thus on the sound celerity in such granular packings.

States of matter

Matter is classically described by the three phase states: solid, liquid and gas. In a solid, particles are closely packed together. Inter-particle forces are strong, hence the atoms can not move freely. This is the reason why, at rest, a solid has a definite shape and volume. An ideal liquid is also a dense state, but unlike the solid, its shape can adapt to the container. Under constant pressure-temperature conditions, its volume remains constant. In contrast to of dense states of matter, a gas in a container occupies the whole volume: the dynamics of its constitutive particles is governed by the kinetic energy and the matter is therefore dilute.

State variables

In thermodynamics, a system at equilibrium is characterized by physical quantities called state variables. These can depend on the size of the system (such as the volume, the mass, the intern energy. ..), in which case they are said to be extensive, or on the contrary, intrinsic to the system (such as the temperature, the density. ..) and they are then said to be intensive. The phase in which the object of the study is, at equilibrium, is characterized by the equation of state that links all the intensive physical quantities together. A very common example of equation of state in the ideal gas law, that links the pressure P (that is intensive) to the temperature (also intensive):

P = ρ R M T
where ρ is the density (intensive, ratio of the two extensive quantities m and V respectively the mass of the sample and its volume) ; M the molar mass (grams per mole) and R the specific gas constant.

The change of state is governed by the evolution of these physical quantities that characterize the system. There are several types of transition.

Phase transitions

A phase transition occurs when one of the state variable of the system crosses a certain threshold value. This transformation is a modification of the physical properties of the system and can be described by thermodynamic tools and statistical physics. From the mechanical point of view, we will show for the dense phases (see next section) that phase changes exhibit very drastic macroscopic and mechanical changes in the behavior of the material. In thermodynamics, for common materials, changes of phase can be presented using the pressure-temperature phase diagram (see figure 1). On this diagram, we can easily see that at a given temperature T , an increase or decrease of the pressure P can lead to the crossing of a state frontier and thus to a change of phase. Similarly, at a given pressure P , a change of temperature can also lead to a change of phase.

Solid

Sometimes a material exhibits a phase transition regarding other physical properties: the Curie point 1 characterizes the temperature above which some materials 1 named after Pierre Curie 0.1. STATES OF MATTER, CHANGE OF STATES lose their permanent magnetism.

In the sample, the change of state of the material can occur in different ways.

There are two main categories for phase transitions.

Phase transitions appear when the characteristic energy of the system -the free energy F -is no longer an analytical function (i.e. is not differentiable, or its derivatives are not continuous) for all the state variables. Paul Ehrenfest2 was the first to propose (in 1933) a classification of the phase transitions 3 . It is based on the continuity of the successive derivatives of the free energy. First order transitions are characterized by a discontinuity (a "jump") of the first derivative of the free energy with respect to one of the state variables. Examples are the transitions between solid-liquid-gas states: the volume, which is the first derivative of F with respect to pressure P changes radically (discontinuously) during the phase changes.

Second order transitions, also called continuous phase transitions, are continuous regarding the first derivative, but their second derivative (with respect to one of the state variable) is discontinuous. A classical example is the para/ferromagnetic transition in iron: at Curie Temperature, the magnetic susceptibility (which is the second derivative of the free energy with respect to the applied magnetic field) change discontinuously, unlike the magnetization (which is the first derivative). In the case of change of phase between fluid and gas, note that there are certain conditions of pressure and temperature for which the transition becomes of second order. Near this critical point, the state of matter is not well defined and the density of the medium fluctuates which scatters light and causes a milky aspect (called critical opalescence). The material is said to be in a superfluidic state.

Spatially in the sample, the two kinds of phase transitions exhibit structural differences. For the first order phase transition, there is coexistence of the two phases, and the phase appearing becomes increasingly dominant with respect to the other one that gradually disappears (the reader can think of melting ice in a pan: liquid water and ice coexist). In contrast, during second order phase transitions, coexistence of the two phases is impossible: the structure of the material changes radically.

This kind of transition, correlation lengths diverge. This gives rise to critical phenomena, and thermodynamic potentials exhibit singularities characterized by critical exponents.

In the framework of phase transitions, the concept of order parameter is key. It is a physical quantity that only depends on the state variables of the system, and whose value discriminates between the phases. For instance, in the liquid vapor transition, the order parameter is the difference of densities. Providing some theoretical hypothesis, once the order parameter ψ is determined, the free energy F of the system can be written as a phenomenological function of ψ , and one can develop it as a Taylor expansion in the order parameter. This is the approach used in our formulation of the fluidity, as discussed in the first chapter of part 1. This approach is said to be "mean-field", as it supposes that ψ is uniform in the system. However, close to the transition, the fluctuations of the order parameter can be important, but this is not the object of the present work.

Thermodynamic and statistical physics give a scope to study and understand phase transitions. However at the macroscopic scale of the system, mechanical aspects are also consequences of changes of the state of matter. In the following, we will only focus on the dense states of matter.

Rigidity transition

Structural aspects

At the atomic scale, the different states of matter are characterized by their structure.

In a solid, molecules are ordered and form structures called crystals. They exhibit long range translational order, unlike liquids for which disorder is dominant.

The pair correlation function provides a useful tool to characterize the structure of 0.1. STATES OF MATTER, CHANGE OF STATES the material. It is defined as the probability of finding a particle at a given distance r away from a given particle taken as reference. In other terms, it describes the average spatial evolution of density within the sample. One can therefore see the existence of a typical length above which the structure is repeated.

1. 1. La transition de rigidité 9 très dense, alors que d'autres le représentaient comme une cristal défectueux. La physique statistique a permis de proposer une nouvelle interprétation de la structure des liquides par l'intermédiaire de la fonction de corrélation de paire g(r) [8,9] . g(r) mesure la probabilité de trouver une particule à la distance r d'une autre particule donnée. La fonction de corrélation de paire est mesurée grâce à la diffusion de neutrons, ou de rayons X [10] . Son allure, pour un liquide, présente un aspect intermédiaire entre ce qu'on peut attendre pour un gaz et pour un solide.

Pour un gaz homogène, la fonction de corrélation de paire est constante et norma- lisée à g(r) = 1, i.e. dénuée de toute structure. La position des atomes ne montre aucune corrélation spatiale. Dans un solide idéal, les atomes sont disposés sur les noeuds d'un réseau géométrique de paramètre p, de sorte qu'on retrouve un atome à chaque fois qu'on parcourt cette distance. Par conséquent, la fonction de corrélation de paire idéalisée est un peigne de Dirac (Fig. 1.6a) qui témoigne de l'ordre structurel à longue distance.

Finalement, dans le cas d'un liquide, la fonction de corrélation de paire présente des oscillations amorties tendant vers 1 (Fig. 1.6b). Le liquide possède une structure à courte portée évoquant celle d'un solide. Cet ordre disparaît au-delà de quelques distances inter-atomiques. Ce point de vue fournit la base d'une théorie fondée sur l'interaction de sphères interagissant à courte portée via un potentiel répulsif [9] , qui s'avère être une bonne approximation des liquides simples et ho- It is therefore worth mentioning that as previously explained, a variation of one of the state variables can lead to a change of phase. This transition can be related to evolution on the microscopic structure of the material. Figure 3 exhibits the disappearance of the structural order during a phase change obtained by the increase of temperature in a solid crystal.

These structural differences between the several states of matter have consequences on the mechanical behavior. And change of states in matter can therefore be studied with a mechanical approach, as we will do in this thesis.

Mechanical aspects

As previously seen, in a molecular solid, atoms are organized. Each particle's position is well defined. Under stress, solids present two types of response: when the applied stress is sufficiently small, the network is slightly deformed and in a reversible way (atoms and molecules recover their initial position) ; when the imposed stress passes this elastic limit, a deformation remains after unloading of the sample. Figure 4 plastic (irreversible)

(remnant) plastic deformation In green, the cycle is reversible. Stress does not rise over the elastic limit σ e . In red, the cycle is irreversible as the stress goes over σ e : the material is plastically deformed and the sample exhibits a remnant deformation.

presents a stress-strain curve. From the initial state, the loading increases in a reversible way: this is the elastic response. Over the elastic yield stress, a loading 0.1. STATES OF MATTER, CHANGE OF STATES cycle exhibit a remnant deformation: the material has been plastically deformed.

When the relation between the stress and the strain is linear and the material isotropic (physical properties are independent of the direction of space), the solid is said to be Hookean and its constitutive equation is of the form:

σ ij = K kk δ ij + 2G ij - 1 3 kk δ ij (1)
with K and G the elastic moduli, respectively in bulk (isotropic compression) and in shear ; δ ij the Kroenecker symbol ; while and σ are the strain and the stress tensors.

Analogous developments can be made for the liquid state of matter. Rheology is the study of the flowing characteristics of liquids. It seeks to relate stresses with strain rates. One can for instance probe the fluid by imposing a shear rate and measure the response in stress, or the opposite (at imposed shear stress, the deformation rate is measured). The constitutive relation of such materials links these two quantities. In the most simple case of incompressible Newtonian fluids, with homogeneous viscosity in the whole sample, this relation is of the form:

σ ij = -pδ ij + µ ∂u i ∂x j + ∂u j ∂x i (2)
with p the hydrostatic pressure ; µ the viscosity of the fluid4 ; u i is the fluid's velocity along axis i and x i is the i -th spatial coordinate. This constitutive relation characterizes the behavior of Newtonian fluids, for which the viscosity does not depend on the shear rate nor evolves with time.

In crystalline solids and in fluids, the mechanical aspects are relatively well known and identified. We have presented above the mechanical properties of solids (that are characterized by their bulk and shear moduli) and liquids (whose flow properties are given by the viscosity). What about the behavior of more "exotic material", such as foam, glass, and sand? These materials indeed share properties of a "liquid" and a "solid" regime. They are dense and like liquids do not present any particular molecular structure, but can present solid behavior such as elasticity. How to characterize such materials and describe their mechanical properties? 0.2 Amorphous media

Definition

The term amorphous comes from the ancient Greek amorphos, "without shape". It refers to the class of matter that is dense but without order (at mid or long range). In an amorphous material, there is structural order, in contrast with solids, which as we already explained exhibit a characteristic length in their structure. More practically, it is the wide family of materials that gathers glass, elastomeres, pastes, emulsions, foams, ... and granular media, which is the object of the present thesis.

Amorphous media are disordered assemblies of particles (atoms, molecules, drops, bubbles, grains...). They share the structural disordered aspect of liquids. However, they exhibit a rigidity transition that does not come with a change in the symmetries of their structure.

Yield-stress fluids

Amorphous media present mechanical properties that diverge from those of Newtonian fluids. One of the most common is the shear rate dependence of their viscosity:

the viscosity decreases when the shear rate increases. This behavior is called shear thinning. It is for instance the case for mayonnaise. Other fluid behaviors exist, such as the increase of the viscosity with respect to the shear rate as in cornstarch suspensions. We hence talk about shear thickening. In some cases, such as for sand, the medium exhibits a yield stress. very active [2,3]. Note that in the present work, we do not focus our interest on the time dependence of the physical parameters of our system. Some studies (e.g. [4,5]) report thixodropic behavior in the rheology of dense granular flows. As long as we need to prepare the granular matter to run an experiment (numerical or laboratory experiment), the history of the system and the thixotropic effects might come into play. However, this is not the object of the studies presented in this work.

Granular matter

The expression granular matter is quite uncommon in everyday life. An explanation for this might be their wide variety and the fact that this type of material can be found in many diverse contexts. From asteroids to sugar, including boulders in the mountains, gravels, sand, cereals, pills in the pharmaceutical industry... Grains -round or not-are everywhere. Remarkably, granular matter represents the second most common state of matter on Earth, after liquid (because of the oceans).

What characterizes this state of matter and unifies such systems under this single expression?

A large number of particles Rigorously, granular matter refers to any material constituted of a great number of particles interacting only through contact forces. By great number we mean an ensemble in which the study of the dynamics of each individual particle is prescribed.

Consider a handful of beach sand. It has a volume of order ∼ 100cm 3 . Say the grain diameter is of order ∼ 0.1mm. We can therefore evaluate the number of grains we have in hand at 100 million5 ! When dealing with such media, it is therefore justified to deal with averaged, coarse-grained physical quantities and to describe the system as continuous as we will do in this thesis. However, regarding the size of the particles, this coarse-grained approach can seem questionable: what is the characteristic size of interest under which this coarse grain approach is no longer relevant? We will see in this thesis some questions about these issues, especially when dealing with the rheology of dense granular flows: where and how can we define a pressure in the medium, such as in a typical fluid? Where can we localize the free surface of an avalanche?... are among the questions we will have to deal with in the first part of this manuscript.

AMORPHOUS MEDIA

Only contacts and mass interactions

It is also important to emphasize the importance of the nature of the interaction governing the dynamics of such system. We just mentioned that particles only interact through contacts (perfectly rigid or deformable, frictional or not...). This means that only contact interactions and inertial (mass) effects are considered, because they are dominant compared to the others (thermal fluctuations, capillary forces, electrostatic interaction between particles...).

Figure 6: Classification of particulate matter as a function of the particle size. Picture extracted from [6].

This statement differentiates between different types of heterogeneous media that are composed of an ensemble of hard particles (see figure 6). Sufficiently large grains (diameter greater than 100µm) only interacts through solid contact interactions.

Their dynamics is dominated by their inertia and the microscopic friction between the particles. Below this size, additional interactions come into play, such as capillary forces if the medium is wet, aerodynamical effects if the motion of the particles is sufficiently important, van der Waals interactions, thermal fluctuations...

For instance, considering glass beads of 1mm diameter such as the one used for the experiments of the present studies, the potential energy of each particle falling from its own size under gravity g is of order mgd (ρ gr ain d 3 )gd ∝ 10 -8 J. In contrast, energy from thermal fluctuations is of order k B T ∝ 10 -21 J in laboratory conditions 6 .

35

Similarly, when the medium is slightly wet (or captures humidity from ambient air, [7]), capillary bridges form between grains. This creates an attractive interaction due to the surface tension of the fluid and the meniscus formed by the fluid between the grain. Condensation around the grain can modify the behavior of the medium, especially in the static phase as we will see in the second part of the thesis. In contrast, in the flowing regime, capillary bridges are constantly broken, thus according to [7] we can suppose that the flowing medium we study does not evolve with time regarding this phenomena.

The medium we consider is therefore only governed by the weight of the grains and contacts (with each other and walls). This means that such systems are dependent on the microscopic friction between the beads. This also leads to erosion (especially in our experiment of flowing grains in a channel, for which we used coated glass beads). We did not focus our study on this phenomena in the present work because we assume that the flowing beads remain perfectly rigid 7 .

A dissipative medium

It is also important to notice that shocks between particles are not completely elastic (the restitution coefficient is not 1: the kinetic energy of two interacting particles is not entirely conserved during a shock). The medium is therefore dissipative. Shocks induce non-linearities in the dynamics of such systems. Another reason for which granular media are very dissipative is their structural nature: within the medium, contacts form chains. The dissipation is due to energy transfer through a growing number of contacts along these chains. In some cases, these force chains to reach the walls and the remaining part of the shock energy is transfered to the experimental container.

thermal interactions start to be comparable to potential energy for grain size under d min = k B T ρ gr ain g

1/4
which gives d min ∝ 1µm, the colloidal limit. 7 We will detail this issue in the section dedicated to the building of a rheology for dense granular flows, first part of this manuscript

LAYOUT OF THE THESIS

Jamming transition

In analogy to the transition between solids and liquids, granular media exhibit a rigidity transition when the number of contacts between the grains in the packing goes below a critical value Z iso . This is based on the Maxwell criterion which states that a system is rigid (resistant to shear and compression) when the number of contacts forces becomes larger than the number of equation characterizing the mechanical equilibrium. This counting argument yields the following exact result for frictionless packing of spheres:

Z > Z iso = 2D, (3)
and for frictional packings:

Z > Z iso = D + 1, (4)
where D is the spatially dimension. When Z = Z iso , the system is said to be isostatic.

This value appears for a lot of material properties as a critical point, either in the liquid or in the solid state.

Layout of the thesis

In this introduction, we have presented the different states of matter and given an overview of thermodynamic and mechanical properties of phase transitions. We have presented the rigidity transition between solids and liquids. We explained how the microscopic structure leads to the mechanical response of these states and we introduce another kind of material, amorphous media, which exhibit features from both the solid and dense states. We then focused and defined granular media, which present both behaviors, "liquid" and "solid".

This thesis proposes two topics of interest, each on either side of the jamming transition that separates the "liquid" and the "solid" regime of dense granular matter.

In the first part, we present a study of the rheology of dense granular flows. We particularly focus on recent models that propose to take into account collective organization of granular matter to describe the flow of granular matter in situations where both the static and the liquid regime coexist. This topic is investigated using two approaches: an experiment of an avalanche flow in a narrow channel, and numerical simulations of flowing hard particles.

The second part of this manuscript presents a study in which the elastic properties of granular jammed assemblies are probed by the means of acoustics measurements.

We built an experimental setup which, installed on board of the Airbus ZeroG, controls the confining pressure on a granular assembly and which can attain very low confinement (∼ 10P a).

Part II

Non-local rheology of dense granular flows Chapter 1

Determining a fluidity for dense granular flows

We present in this chapter a review in which we discuss the concept of fluidity in granular matter. We first introduce and define the local, µ(I)--rheology and then expose its lacks. We define the terminology non local and review different approaches that attempt to model non-locality. In this context, we present our approach in which we proceed by a gradient expansion of the constitutive equation. We explain the choice we make in identifying the fluidity parameter f as the inertial number, as opposed to the elasto-plastic approach for instance. We report the conceptual differences between these constitutive relations proposed for granular matter. In the limit of rigid particles, we present a linearization of the non-local equation in the case of homogeneous shear stress profile. This approach will be used in the following chapters in order to compute the complete flow profiles of the configurations we investigate.

Introduction

Since non-locality was introduced as an interpretive framework for dense granular flows [1][2][3][4], it has become a key concept to describe the rheology of complex fluids in soft condensed matter. However, the connections between the various contributions to this subject, their similarities and possible conflicts need clarification. In particular, among the pending questions that must be answered, a fundamental and vivid issue is the possible emergence of non-locality as the signature of a dynamical phase transition [5][6][7][8][9]. This interrogation does not only concerns granular matter but should be apprehended in the more general context of amorphous solids undergoing a rigidity transition. At present, different conceptual approaches have been put forwards to describe non-locality and several non-local constitutive relations were proposed for granular matter. It is thus fair to ask whether these approaches are equivalent and to which extent, for example, they are similar to phase field models [10][11][12][13] built on an underlying liquid/solid phase transition. Also for granular matter, shear-banding and apparent "creeping zones" are observed which are difficult to re-conciliate with a simple local rheology [14], and this has been the starting point of different propositions for non-local constitutive relations. For many, the elements of proof validating ⋆ Contribution to the Topical Issue "Multi-scale phenomena in complex flows and flowing matter" edited by Luca Biferale, Massimo Cencini, Alessandra Lanotte and Mauro Sbragaglia.

a e-mail: andreotti@pmmh.espci.fr these approaches has often been a mere "good fitting" of the velocity profiles. A legitimate question is then to ask whether this is sufficient to demonstrate the validity of a particular model and moreover, what could be other more stringent tests providing essential information on the dynamical mechanisms responsible for non-locality.

Our aim here is to propose a critical discussion of the concepts of non-local rheology and fluidity in dense granular matter, based on recent progresses as well as older results. In the next section, we first review the rheology of dense granular flows, starting from the local rheology towards evidence for non-local effects and describing nonlocal approaches. In sect. 3, we discuss the concept of fluidity. Section 4 is devoted to the differences among the nonlocal constitutive relations proposed for granular matter. We end the paper in sect. 5 with a discussion on further possible tests that must be performed to better understand the mechanisms at the origin of non-local effects.

2 On the rheology of dense granular matter

Rigidity transition

When sufficiently polydisperse to avoid crystallisation, a granular packing at rest can be considered as an amorphous solid. By definition, amorphous solids refer here to systems that may resist to a shear stress while they do not present any long-range translational order at the microscopic scale, namely the grain size for granular matter.

Let us consider, for clarity, an ideal rheometer in which the material is submitted to a homogeneous shear stress σ. The system behaves mechanically as a solid if it reaches equilibrium at a finite strain γ. It is considered as elastic if it returns to its original state, once the stress is removed and plastic, otherwise. Conversely, the system behaves mechanically as a liquid if it flows permanently at a finite strain rate γ. The system exhibits a rigidity transition if its dynamical behaviour switches from solid-like to liquidlike, when a control parameter crosses a threshold value. Most soft amorphous solids present a rigidity transition upon varying the shear stress, the threshold value being named the yield stress σ y . We hence define the yield parameter as the ratio of the shear stress to the yield stress

Y = σ σ y . (1)
We emphasise here that the existence of a yield stress is not an intrinsic material property. It depends on the other control parameters that are kept constant during the loading. As an example, a granular material, dry or suspended in a fluid, displays a yield stress when the particle-borne pressure P is imposed, whereas it does not when the volume fraction φ is imposed. This is obviously a key issue and often a source of confusion. This question was discussed recently in the context of different amorphous particulate materials [15]. The rigidity transition is intimately related to the multi-stability of the energy landscape: the system has to cross energy barriers to flow. The physical nature of the mechanisms preventing irreversible plastic deformations allows to classify the soft amorphous materials and their corresponding rigidity transition as:

-Entropic for glasses formed by thermal quenching [16][17][18]; the rigidity transition is then called "glass transition" [19].

-Enthalpic for soft elastic particles at high volume fraction [20][21][22][23][24]; the free energy may result from capillarity (foam, emulsion), from electrostatics, from particle elasticity, etc. The rigidity transition is then called elasto-plastic depinning transition [25][26][27].

-Geometric for hard grains submitted to a confining pressure [28]; the rigidity transition is then called jamming transition [29].

Although we focus here on the last case, we will frequently discuss connections with other complex fluids in soft matter.

Local rheology

Following the seminal paper by GdR MiDi [30], major improvements were obtained to provide a consistent framework to understand and model how granular matter flows.

In the rigid limit, granular matter does not present any intrinsic energy scale and confining pressure P thus provides the only relevant scale of energy per unit volume. As a consequence, P sets the yield stress as σ y = µ c P , where µ c is the "critical" friction coefficient, which depends on microscopic material properties (e.g. packing polydispersity, inter-granular friction, shape, etc). For a real granular material, the rigidity transition is actually subcritical, a property associated with the presence of inter-granular friction, and the hysteresis of the effective friction coefficient has remained unexplained up to now. For grains of mean diameter d and mass density ρ g , the confining pressure also sets the time scale T = d/ P/ρ g for plastic reorganizations at the granular level (microscopic time scale). Following [30,31], one can define the rescaled strain rate, or inertial number, as

I ≡ γT = γd P/ρ g . (2)
Writing the yield parameter as

Y = σ σ y = σ µ c P , (3)
the constitutive relation for homogeneous steady flows takes the generic form

Y = µ(I) µ c = 1 + aI n , (4)
n = 1 for grains presenting a standard friction coefficient ≃ 0.5 at contact, and n = 1/2 for frictionless grains. This relation must be complemented by a law relating the volume fraction φ to the rescaled strain rate I φ -

φ c = -bI n , (5)
with the same phenomenological exponent n; a and b are constants that depend on the microscopic details of the system. This derivation is simply based on dimensional analysis, in the rigid limit where the grain elasticity is irrelevant. Empirical measurements indeed show a frictional behaviour with the emergence of a yield stress. Relation (4) is well suited to investigate pressure-controlled situations. However, the same equations can entirely be recast to handle situations where φ is fixed, and then the yield stress disappears from the constitutive picture. Inverting eq. ( 5), one obtains

P = b 2/n ρ γ2 d 2 (φ c -φ) 2/n , (6)
σ = µ c b 2/n 1 + a b (φ c -φ) ρ γ2 d 2 (φ c -φ) 2/n . (7)
In this representation, the shear stress clearly vanishes in the limit γ → 0, i.e. if φ < φ c . From this derivation, we can see that, even in the limit of rigid particles, i.e. without any explicit elasticity, a granular material is compressible under shear [32]. The pressure therefore requires time to establish over the size of the system. Provided this time scale remains short compared to γ-1 , the pressure can be considered as defined "instantaneously" and thus can be used as a state variable instead of φ. This is indeed a central assumption to neglect compressibility (Boussinesq approximation) when using the µ(I) formulation in heterogeneous situations: it implicitly requires that pressure is established macroscopically over a very short time scale and varies slowly in time and space.

Using this close set of constitutive relations, a quantitative agreement with numerical or experimental measurements has been reached in different configurations. In particular, for avalanche flows of glass beads on a rough inclined plane -an important situation as the yield parameter Y is fixed by the inclination angle-Pouliquen [33] has derived an effective flow rule consistent with the µ(I) rheology. This has led to a three-dimensional extension of the local rheology [34], yielding the correct scaling laws characterising the chute flow geometry.

Dense granular suspensions have been shown to follow the same rheology, with T = η f /P , where η f is the viscosity of the suspending fluid, to form the so-called viscous dimensionless number T γ [35,36] -sometimes noted J. More generally, the rheology obeyed by a granular material in a homogeneous steady state takes the very same form as observed for soft material presenting a yield stress. In those more general situations, the shear rate γ can then be rescaled by a plastic time scale T , to form a dimensionless number akin to I ≡ γT . This kind of rheology takes the same generic form (4), often called a Herschel-Bulkley constitutive relation, for all these systems.

Failure of the local rheology

Consider now a heterogeneous shear flow. Its rheology is said local if the stress tensor at a given location is still a function of the shear rate at the same place. Non-locality refers to any deviation to such a local constitutive relation. Before giving precise examples, let us discuss the choice of such a name. It has been introduced in granular material to describe the distant transmission of momentum through the granular skeleton, during collisions [1][2][3][4]. In the limit of rigid particles, this transport of momentum is instantaneous so that a stress, which is a flux of momentum, can be induced by distant collisions.

In hydrodynamics, the epitome of non-locality is pressure. In a simple fluid, pressure is transmitted at the speed of sound. At low mach numbers, the time needed for a pressure signal to cross the entire flow is small compared to γ-1 . In this limit, pressure is determined by the incompressibility condition ∇ • v = 0, where v is the velocity field. Taking the divergence of the Navier-Stokes equation one therefore obtains, for a Newtonian fluid, the Laplace equation

∇ 2 P = -ρ∇ • (v • ∇v). (8)
Pressure balances the potential part of inertial terms. The Biot-Savart equation provides an explicit solution of this equation under an integral form. The fact that pressure is a function of the whole velocity field, and not only of the local strain rate then appears explicitly. As a conclusion, there are a priori two definitions of non-locality which are not equivalent: i) Momentum is transported over large distances on a time scale small in comparison to γ-1 and to the plastic rearrangement time scale T .

ii) The constitutive relation involves a second state variable, which is not a function of the strain rate, and whose evolution is controlled by an independent equation, typically involving a Laplacian operator.

For instance, the kinetic theory, which is valid for dilute and rapid granular flows, is non-local in the weak sense ii), since it introduces an independent field representing the mean squared velocity fluctuations (the socalled granular temperature), which may control the stress tensor [39]. However, the transmission of momentum remains perfectly local in the sense i).

Non-locality in the weak sense ii) manifests itself through different properties. The first one is the evidence of a creeping flow in regions below the yield condition (Y < 1) [30,[40][41][42][43]. Instead of the expected static zone (i.e. a solid), one observes an exponential spatial relaxation of the shear rate γ. A second property is the fact that the yield conditions are sensitive to the system size and to the boundary conditions. In the case of granular matter, the yield stress measured on an inclined plane depends on the thickness of the deposit [30,33]. Also, in conditions where the grains should flow according to the local rheology, jammed regions are identified below the flow [44]. Furthermore, for self-channelised flows close to jamming, quantitative departure from the local rheology predictions are explicitly shown [45]. Finally, properties pointing on the existence of non-local effects are revealed by micro-rheology experiments. For example, the force-velocity relations assessed by an intruder plunged in the material strongly depend on the presence of a distant shear flow [46][47][48], or of a vibrating boundary [49]. Examples of manifestation of non-locality for soft and granular matter are given in fig. 1 and fig. 2.

Reviewing non-local models

We review here the main approaches that were put forwards to tackle these problems.

Cosserat approach

In the framework of plastic theories developed for soil mechanics, the formation of localised shear bands is sometimes apprehended via a Cosserat extension of the continuous elasto-plastic theory (see for example [50] and references therein). To describe the quasi-static state of deformations of a granular material, new fields are introduced that couple stresses and rotations. This theory introduces a microscopic length scale describing the range of influence of the microscopic granular rotations and provides a non-local coupling for plastic deformations over this scale. However, in spite of the fact that it may provide a useful regularisation technique for numerical computation methods, its use is often seen as limited since the issue of assessing objectively the constitutive parameters and providing Local rheology would predict a linear velocity profile. The numerical data are from [37]. (c) Experimental (triangles) and numerical (squares) velocity profiles of a foam flow in a Couette cell, data from [38] and [72], respectively. Local theory would predict a localised failure at wall. Fig. 2. Evidences of non-local behaviour for granular matter. (a) Experimental measurements of the stopping (black symbols) and starting (white symbols) heights on an inclined plane with glass beads, data from [30]. The dependence of stoppage and onset of flow with the granular height are manifestations of the non-local character of the rheology. (b) Normalised creep velocity of an intruder in a Couette cell filled with glass beads as a function of the normalised pulling force on the intruder, data from [47]. Data are taken in conditions where the intruder should be blocked (in the solid phase) in a local rheology picture based on a Coulomb yield criterion. (c) Experimental velocity profile of a dry granular avalanche flow in a narrow channel with frictional lateral walls. The dashed line is the prediction of the local rheology displaying a depth z where flow stops whereas the granular continues to flow below this limit.

consistent boundary conditions for the fields, has so far remained a shortcoming of the approach.

Phase field approach

As stated above, non-locality, in the weak sense, reflects the existence of a state parameter, beyond the strain rate, determining the stress values. As this parameter measures how fluid the system is, following Derec et al. [51], we will refer to it as the "fluidity" and will note it thereafter f . We warn the reader that in different papers, fluidity may refer to different physical quantities. Here, we keep the name and its conceptual definition in relation to non-local rheology. Importantly, we consider that for our purpose, a relevant fluidity has to be selected on physical basis, amongst all state variables.

From a phenomenological point of view, fluidity plays the role of an order parameter describing the dynamical transition from solid-like to fluid-like behaviour. It was first proposed by Aranson and Tsimring to introduce a phase field f which vanishes in the static state and which tends to 1 in the fully fluidised state [10]. In this approach f is therefore dimensionless. The overall shear stress is then formally decomposed as the sum of a solid-and liquid-like contributions weighted, respectively, by 1f and f . Following Landau standard derivation, the order parameter f is controlled by a reaction diffusion equation of the form

T ḟ = I(f ) + ℓ 2 ∇ 2 f, (9)
where I is a function of f parametrised by the state variables and in particular by the rescaled shear rate I. Note that I(f ) can be designed to reflect a subcritical, hysteretic transition from solid-like to liquid-like behaviours, as generically observed for granular matter. The microscopic time T is a characteristic time for fluidization to occur and ℓ an elementary length scale. The diffusion coefficient of the fluidity is ℓ 2 /T . The Laplacian operator results from a gradient expansion, assuming spatial isotropy. This term is the transcription of a (weak) non-locality. As a matter of fact, in the steady state, f is determined by a non-linear Laplace equation, just like pressure in hydrodynamics. This approach produces an effective rheology different from µ(I) in the sense that Pouliquen's flow rule [33], valid for avalanches of spherical beads, is not recovered. However, it yields a real solid/fluid phase transition semiquantitatively close to what is observed for sandy grains, in particular the starting and stopping heights, flow rules and erosion/deposition waves [13,[52][53][54]. Nevertheless, it does not reproduce creep zones close to a shear band.

The elasto-plastic approach

The approach proposed by Kamrin et al. [55] to model granular flows is directly adapted from the Kinetic Elasto-Plastic (KEP) model introduced by Bocquet et al. [6] for soft matter. The key concept is the fluidity, which can be defined in a unified way as

f = γ Y . ( 10
)
We will devote below an entire sub-section to this approach.

Mechanically activated plastic events

An original idea to describe non-locality has been proposed by Forterre and Pouliquen [56]. It is based on an analogy with Eyring's transition state theory for the viscosity of liquids, where mechanical fluctuations -introduced here as a synonymous of heterogeneitieswould play the role of temperature in thermal systems.

Plastic rearrangements occur at a rate proportional to the strain rate γ. They are assumed to generate at random a new realisation of the forces on the contact network, allowing for the formation of new weak zones where the next rearrangement will occur. At a semi-quantitative level, this approach improves significantly the local visco-plastic approach as it proposes a physical hint for microscopic processes inducing nonlocality for granular media. Moreover, it can predict dependence of the stopping angle on flow height, shear bands extension increasing with the flow rate. However, the theoretical outcomes are more difficult to quantitatively reconciliate with Pouliquen's flow rule for the chute flow and more importantly, it does not predict a thickness dependence of the avalanche starting height. In the non-local formulation proposed by Pouliquen and Forterre [56], the shear rate γ obeys an integral equation, which involves an exponential kernel, function of the stress tensor and of the distance, interpreted as a Boltzmann-like factor. Provided the fact that for granular packing the stress fluctuations take place generically over few grain sizes, the authors assume a spatial dependence of the interaction kernel as a Lorentzian function decaying algebraically (power -2) over this granular size.

The relation between the shear stress and the strain rate depends non-locally on two fields i.e. on two fluidity parameters: the strain rate γ and the confining pressure P . The relatively fast decay of the chosen spatial kernel makes possible a long-wavelength expansion of the rearrangement rate equation with respect to γ and P . If P varies slowly at the scale of the grain size, the expansion generates at the first order a Laplacian operator, like other models. More precisely, the constitutive relation presents a dependence on ∇ 2 γ, suggesting that the fluidity f is the rate of plastic events γ.

The gradient expansion of the constitutive relation

We have ourselves followed an approach which is significantly different from the previous ones, and which suggests possible candidates for the most appropriate fluidity of dense granular flows [37,[START_REF] Bouzid | Comportement rhéologiques et effets non-locaux dans les écoulements granulaires denses[END_REF]. Imagine the problem solved and a fluidity f built, which vanishes in the solid phase and increases with the ability to flow. We then follow the standard vision of Maxwell rigidity transition as put forwards to understand jamming in granular matter. Due to the cooperative motion of particles along soft modes, flowing is facilitated when at a given point it is surrounded by a more fluid zone. Conversely, the resistive stress is larger when the point is surrounded by a more solid neighbourhood. Experimentally, the effect is particularly significant close to the jamming transition, when f vanishes. Therefore, one needs to define a relative fluidity, which compares the degree of fluidity at one point and in its vicinity. Assuming as before that the influence of fluidity is statistically isotropic and results from short-range interactions between shear zones, the relative fluidity can be defined as

κ = ℓ 2 ∇ 2 f f , ( 11
)
where ℓ is a length on the order of few grain diameters. The Laplacian is indeed the lowest-order operator appearing in a systematic expansion in a functional of f . The relative fluidity κ remains finite when f goes to 0. The constitutive relation can be expanded around the relation Y = µ(I), valid in the homogeneous case, according to

Y = µ(I)χ(κ), with χ(κ) ≃ 1 -κ + O(κ 2 ). ( 12
)
κ is positive when the point considered is surrounded by a more liquid region (higher f ). This region flows more easily than expected from the local value of f , and the corresponding shear stress is therefore lower. In this formulation, χ(κ) must thus be a decreasing function of κ, which justifies the minus sign in front of κ in eq. (12). Note that the lack of multiplicative factor in eq. ( 12) defines ℓ in a univocal way. Importantly, this phenomenological derivation does not depend on the nature of the mechanical interaction between the shear zones; the reader may think of the analogy with the van der Waals gradient expansion of the Helmholtz free energy at a liquid-vapour interface.

There are three obvious choices for the granular fluidity. One would be to introduce the (coarse grained) number of contacts per grain Z and number of sliding contacts per grain ζ. In the spirit of Maxwell rigidity transition theory, the fluidity would then be defined as the distance to isostaticity. A second possibility is to introduce the volume fraction φ, and to build the fluidity as the distance φ cφ to the value φ c that the volume fraction reaches in the limit I → 0. This poses two problems: first, φ c depends on the fraction of sliding contacts; second, one would need to introduce the mass conservation equation for φ i.e. to consider explicitely the granular fluid as compressible. Finally, there is a last quantity that can play the role of the fluidity: the inertial number I itself. It vanishes in the solid state and increases with the degree of fluidity. Ockham's razor -law of parsimony-is obviously in favour of such a choice, as the equations are closed without involving further equations, in the dense limit φ ≃ φ c . This does not constitute a deep scientific argument, except that simple models are better testable.

The constitutive relation (12) must actually be complemented by another one for the volume fraction. Assuming that φ is a local function of I, expressing the nonlocal rheology with f = I with κ = ℓ 2 ∇ 2 I/I, and with f = φ cφ with κ = ℓ 2 ∇ 2 φ/(φ cφ) are mathematically analogous. The only subtlety is that compressibility must be taken into account if the non-local constitutive relation is expressed in terms of φ while the flow can be considered as almost incompressible (i.e. in the Boussinesq approximation) if the non-local constitutive relation is expressed in terms of I. In the later case, the pressure P is instantaneously determined and becomes a state variable.

To conclude this section, we rewrite the non-local rheology (12), with the choice of f = I = T γ as a fluidity, as

0 = IY µ(I) -I + ℓ 2 ∇ 2 I. (13)
This formulation, directly derived from eq. ( 12), makes it easier to compare with the other approaches, as discussed below.

3 The fluidity concept in soft elastic material

Since fluidity was introduced to describe granular flows, it is important to review the physical basis of the recent advances made to render the complex rheology of soft matter (foams, gels, emulsions, etc) using this concept This auxiliary equation reflects the microscopic or mesoscopic processes at the heart of the physics that one seeks to describe. We shortly review the propositions made in the literature to define fluidity.

Fluidity and elastic deformation

Many complex fluids are visco-elasto-plastic. It is the case of polymer melts, micelles and lamellar surfactant phases, which are subject to shear-banding, but also the case of foams and emulsions, which are yield stress fluids. In this situation, an obvious state variable controlling the rheology is the elastic strain ϵ. The elastic deformation is an intrinsic property describing for example, the polymer chains extension or for a foam, the current bubble deformation state. Importantly, ϵ is a field coarse grained in space, in time, or averaged over realisations. The total strain rate γ, determined from the subsequent positions of the elements, can be formally written as the sum of the elastic strain rate ε and of a plastic contribution, equal to γε. For convenience, we keep a scalar description -the extension to a tensorial form does not present any conceptual difficulty. Non-locality in visco-elasto-plastic models was introduced by Olmsted [14] who has added non-local terms terms to the Johnson-Segalman model [START_REF] Johnson | [END_REF]. As in the phase field approach, the overall stress is decomposed into a fluid-borne stress, for instance Newtonian η γ, and an elastic-borne stress Gϵ, where η is a viscosity and G an elastic shear modulus. The equation governing the evolution of the deformation ϵ reads

T ε = T γ -I(ϵ) + ℓ 2 ∇ 2 ϵ, (14)
where T ∼ η/G is the relaxation time scale of the components (e.g. the polymer). The diffusion term ℓ 2 ∇ 2 ϵ is responsible for non-locality. The function I reflects the elastic properties of the components and determines the constitutive relation measured in a homogeneous steady state. For the Herschel-Bulkley relationship (4) (we recall that I = T γ), one obtains

I(ϵ) = 1 a max 0, Gϵ σ y -1 (1/n) . ( 15
)
The standard Johnson-Segalman model corresponds to the exponent n = 1. The ratio ϵ y = σ y /G is the yield strain above which plastic events nucleate. The same model has successfully been used to describe complex fluids subject to shear-banding, with a vanishing yield stress [14,59].

In the steady state, the elastic deformation obeys the non-linear Laplace equation ℓ 2 ∇ 2 ϵ -I(ϵ) = -T γ, which leads to exponential relaxations in space (see sect. 4). In the limit where T becomes much smaller than γ-1 , the same equation holds at all times and the dynamics becomes truly non-local in the sense i).

The above relation was modified by Marmottant and Graner to model dry foams [60], assuming that the elastic part does not evolve over to the internal time T but over the time scale γ-1 , which is assumed to provide the only time scale of the problem. This assumption is close to that made by Forterre and Pouliquen [56] for granular flows. Keeping the structure of eq. ( 14), the governing equation then takes the form

ε = γ 1 -I(ϵ) + ℓ 2 ∇ 2 ϵ , ( 16
)
where I is a function which is essentially 0 at low ϵ and which sharply increases and crosses 1 at the yield strain σ y /G. Note that the derivation and tests proposed in [60] actually do not take the non-local term into account (they have ℓ = 0). We introduce it here, to clarify the connections between the various models.

Fluidity as a Maxwell relaxation rate

Considering again the ideal linear Couette cell controlled at imposed shear stress σ, the mechanical response of a soft material usually presents a transient in time. Over time scales comparable to γ-1 or smaller, the structure does not have time to evolve: the time dependence is said to reflect visco-elasticity. If the rheology evolves over time scales long compared to γ-1 , the system is said thixotropic.

= -f ϵ + γ. ( 17
)
The auxiliary equation that governs the evolution of f is macroscopic and implies ageing and shear rejuvenation processes. By contrast, Olmsted model assumes a constant time for the relaxation and Marmottant and Graner a time scale inversely proportional to γ. The three models thus differ by the identification of the relevant strain relaxation time. Interestingly, in the solid phase (i.e. below the Coulomb threshold) sheared granular packing display ageing that can be described by this equation and the fluidity parameter f is directly related to the rate of localised plastic events, called hot spots that were directly visualised [61]. We will turn back to this point later.

Fluidity as the rate of plastic events

Close enough to the yield stress, flow occurs in concentrated emulsions and foams via a succession of reversible elastic deformations [62,63] and avalanches of irreversible plastic rearrangements (also called shear transformation zones in the literature [64]). Such localised plastic events induce a long-range anisotropic relaxation of the elastic stress over the system, which constitutes an obvious source of non-locality. Based on this observation, the fluidity in a class of models has been associated to the rate of plastic events, coarse grained in space and time. The physical picture assumes that the material is essentially in an elastic state under stress but due to disorder or temperature or mechanical fluctuations, localised plastic events nucleate, which induce local irreversible stress relaxation processes.

The resulting deformation is absorbed elastically by the medium. This is at the origin of a global plastic deformation rate. For a given fast and local relaxation process, other plastic events may be triggered. When there is a continuous rate of coupled plastic events spanning a significant time, this is called an avalanche [25][26][27]73]. This avalanching process occurs preferentially at larger stress.

The exact account for avalanche dynamics as a net contribution to the final plastic deformation rate is a difficult issue and has therefore been, in most models, phenomenologically modelled by means of non-linear terms entering the fluidity equation.

The auxiliary equation for f can be a stochastic Smoluchowsky equation [65] but macroscopic -essentially mean-field-models have also been derived explicitly. For example Bocquet et al. [6] have proposed a model, called KEP (kinetic elasto-plastic), which adds a spatial coupling to the probabilistic framework introduced by Hébraud and Lequeux [65]. The model describes the evolution of the probability to observe locally a certain local "shear stress" and furthermore assumes that plastic events are triggered above a non-fluctuating "local yield stress". The macroscopic rheology turns out to be controlled by the behaviour of the shear stress probability distribution function in the immediate vicinity of the local yield stress: in a steady homogeneous state, a Herschel-Bulkley rheology is recovered.

In the KEP model, plastic events lead to a noise around them that helps to trigger other plastic events. Non-locality therefore appears as a dependence of a "local stress diffusion" on the rate of plastic events around the region considered. Both the elastic strain (and therefore the stress) and the rate of plastic events stem from an integral over the same probability distribution function. In the steady state (but not during transients), they can be related to each other. One recovers the relation between the rate f of plastic events and elastic strain of eq. ( 17), for ε = 0

f = γ ϵ . (18)
The non-local equation has been derived by Bocquet et al.

in the steady state (and only in this case) and takes the form of a non-linear Poisson equation, analogous to that of Olmsted's model (14) in the steady state

0 = γ -T -1 I(f ) + ℓ 2 ∇ 2 f. (19)
We later refer to this equation as the KEP constitutive relation or KEP model. However, a mathematical analogy is not a physical equivalence: the fluidity defined as the variable that appears in the non-local Laplacian term is, in one case, the elastic deformation ϵ and in the other, its relaxation rate f . Furthermore, to fully solve the problem in association with a Laplacian term in the formulation, it is necessary to provide boundary conditions for fluidity [66]. The choice of the expression for the fluidity must then be consistent with the physical boundary conditions. We will turn later on the possible tests to determine which fluidity is the relevant one.

Importantly, the macroscopic emergence of the auxiliary fluidity equation involves a Laplacian operator which physically represents the spatial range of plastic relaxations. Mathematically, it is the source of non-locality in the constitutive relation. A possible issue is that the stress relaxation induced by plastic events can have in general an anisotropic character, even in a statistical sense, which is not reflected by the isotropic Laplacian term in (19). If, in this triggering process, anisotropy is important, higherorder terms must be included in the spatial expansion of the stress propagator in [6]. Also, close to the yield condition, long-range avalanches may take place and spatial coupling can span large distances that eventually diverge at the yield point. In this limit, the Laplacian, which is essentially a mean-field operator, is unlikely to capture nonlocality, as it is well known that Landau-like approaches generically fail in the vicinity of a critical point.

Fluidity as the inverse of viscosity

Bocquet et al. [6] have made a further, apparently innocent, step: as the elastic stress σ is proportional to the elastic strain ϵ, the structure of the equations does not change if f is defined as

f = γ σ (20)
instead of eq. ( 18). The fluidity f becomes the inverse of the particle-borne viscosity. The equation governing f remains of the form (19). We will discuss below the problems associated with the change of variable ϵ → σ. We emphasise again that the relation between the rate of plastic events and the fluidity is not trivial and must be tested.

To handle granular matter, Kamrin et al. [55] have proposed to rescale this inverse viscosity by the yield stress σ y . If σ y is constant, this does not change the shape of the equation governing the new fluidity (eq. ( 19))

f = σ y σ γ. ( 21
)
The fluidity f is then homogeneous to a strain rate. Again, the function I can be determined from the constitutive relation measured in a homogeneous steady state. For the Herschel-Bulkley relationship, one obtains the implicit equation

I(f ) 1 + aI n (f ) = T f, (22)
whose solution takes the form

I(f ) = T f + a(T f) 1+n + O(f 1+2n ). ( 23
)
For the particular case of frictional granular material, for which n = 1, one gets the analytical solution

I(f ) = T f 1 -aT f . ( 24
)
4 What are the differences between non-local constitutive relations proposed for granular matter?

4.1 Fluidity, from soft-matter to granular matter

As seen above, the fluidity in Kamrin et al.'s model for granular matter [55], derived from the results of Bocquet et al. [6], writes

f = σ y σ γ = γ Y = µ c P σ γ. (25)
The auxiliary equation for this parameter is derived from eq. ( 19) by a linearisation that we discuss below in details, as it is problematic.

The transposition from soft matter to rigid grains poses a fundamental issue. The fluidity must obviously be a state variable. By state variable, we mean that the fluidity f must be a coarse-grained field (in space and time) which can be determined from the state of the system. It can for example depend on the strain, the strain rate, the volume fraction, the mean number of contacts. For a granular material, it can also involve the fraction of sliding contacts or the orientation of the contacts. For foams or emulsions, it can reflect the elastic deformation of elementary cells. However, it cannot depend explicitly on the stress tensor, which is not a state variable itself. This directly results from Newton second's law, which tells that positions and velocities of the particles determines the state of a mechanical system, from which forces are derived. Similarly, the (non-local) constitutive relation must relate the stress tensor to the state variables -and not the opposite. f as defined in eq. ( 25) was a state variable for an elastic system, because it was fundamentally based on the elastic deformation. It is no longer the case for a granular system composed of rigid grains.

Does flowing granular matter exhibit plastic events?

The transposition of models derived for elasto-plastic material to granular systems obviously requires that granular matter behaves as hypothesised by elasto-plasticity. The elasto-plastic picture assumes that the material behaves most of the time like a solid, but presents local and shortlived plastic events [26,64]. The associated scenario is a localised rupture initiation followed by a scale-free avalanche of localised events. In order to investigate whether this picture constitutes an alternative to the jamming scenario to interpret the non-local nature of the granular rheology, we compare, by means of numerical simulations (Discrete Element Method), the dynamics of two otherwise identical systems composed of hard and soft grains.

The general numerical set-up is that used in [37]: we consider a two-dimensional system composed of ∼ 2.10 3 spherical particles of a mean diameter d, with a ±20% polydispersity. Such a choice ensures that the sample will not crystalise. The particles can interact through contact forces modelled as a viscoelastic force along the normal contact direction and as a Coulomb friction along the tangential direction. The corresponding coefficient of restitution is e ≃ 0.9. The Coulomb friction coefficient is set to µ p = 0.4 for fictional particles and µ p = 0 for the frictionless system. The particles are confined in a plane shear cell composed of two rough solid walls made by the same particles, glued together. Periodic boundary conditions are used along the shear direction x. The position of the wall is controlled in order to impose a constant normal stress P and constant and opposite velocities of the walls along x. The system is in the asymptotic rigid limit when the ratio k n /P of the normal spring constant with the pressure is sufficiently large (typically above 10 3 ).

The presence of localised plastic events is usually based on a visual inspection of different fields. The squared deviation from an affine deformation on a local scale has for instance been proposed as a field indicating plastic activity [64,67]. However, such a quantity, as well as all those based on the non-affine velocity, characterises fluctuations around the mean flow, and not the local contribution of a certain area to the mean flow. We wish here to propose a practical definition of these events, based on the quantitative criterion that they must be separated in time and localised in space. Importantly, to match their role played in elasto-plastic models [6,64,67,68], they must also contribute additively to the average shear rate γ.

In order to detect localised plastic events, we have built a coarse-grained field Γ (r, t) reflecting, at time t, the local contribution to γ of a small region around the position r. We impose that the time average of Γ must everywhere give γ. A coarse-graining method similar to that proposed for the stress tensor [69,70] is adapted here to the computation of velocity differences and we take:

Γ (r, t) = N j=1 [u i (r, t) -u j (r, t)][z i (t) -z j (t)] exp -||∆r|| 2 2δ 2 N j=1 [z i (t) -z j (t)] 2 exp -||∆r|| 2 2δ 2 , ( 26
)
where u i (r, t) is the velocity of the grain i at the time the distance between the grain i and j. δ is the coarsegraining length, typically on the order of the grain size d. We display in fig. 3 and fig. 4, for a system of rigid and soft grains respectively, the map of the local contribution Γ to the shear rate γ at different times. In fig. 5a,b we show the corresponding spatio-temporal diagrams built on the central line of the cell. We observe contrasted behaviours in the two cases. In the soft system, nothing much happens most of the time, except for short periods of intense activity, associated with a cascade of plastic events. Conversely, the hard system presents more moderate but permanent fluctuations even for asymptotically small γ.

t and ||∆r|| = [z i (t) -z j (t)] 2 + [x i (t) -x j (t)] 2 denote
To make these observations quantitative, fig. 5c shows the probability distribution function (PDF) over time of ⟨ Γ ⟩, which is the space average of Γ over the cell. In panel (d), we similarly display the PDF of the spatial standard deviation δ Γ . The hard-particle system presents a narrow Gaussian distribution of ⟨ Γ ⟩ around γ, while the PDF corresponding to the soft system shows stretched tails, which are due to a very intermittent behaviour associated with these plastic events. The PDF of δ Γ provides informations about spatial heterogeneities in the system. The peak of the black line around 10 γ in fig. 5d indicates that they are large and permanent in the hard system. For the soft system, the PDF shows an algebraic decay, which means that the field Γ is homogeneous most of the time. However, when the computation of δ Γ is restricted to the periods of time where plastic events occur (periods where |⟨ Γ ⟩| is larger than a given value, here 5 γ in fig. 5d), its PDF also presents a peak: in the soft system, plastic events are associated with a very heterogeneous field of Γ . Conversely, an assembly of rigid particles does not present local plastic events when sheared permanently. Its dynamics is not intermittent but presents spatial heterogeneities.

In conclusion, when constituted of rigid particles, sheared granular systems do not present a succession of elastic energy accumulation and sudden release. Their dynamics rather show permanent cooperative motions. As a consequence, approaches explicitly based on elasto-plasticity developed for soft systems, such as those discussed above, cannot be transposed to granular flows, where elasticity of the grains is irrelevant. The physical foundation Kamrin et al.'s approach for granular matter [55] seems already in this perspective, extremely problematic.

Linearisation in the case of a homogeneous shear stress profile

The simplest situation in which the model predictions can be tested is a shear cell inside which the shear stress, and therefore the yield parameter Y, are homogeneous.

In such a situation, a local rheology predicts a constant shear rate, which, once rescaled by T , is denoted I ∞ . This cell is driven by boundary layers on each side, whose properties are not necessarily those of the bulk. How to realise this in practice for numerical simulations is for example described in Bouzid et al. [37].

We first consider the case where I ∞ does not vanish. By definition, we have Y = µ(I ∞ ). Making profit that Y is constant, one can linearise the equations in I around I ∞ . Using the gradient expansion model (eq. ( 13)), where I is the fluidity, we obtain

ℓ 2 ∇ 2 (I -I ∞ ) = n (1 + aI n ∞ ) -1 -1 (I -I ∞ ) = n Y -1 Y (I -I ∞ ). ( 27
)
Denoting by z the axis transverse to the flow, we obtain exponential solutions of the form

I = I ∞ + A + exp(z/L) + A -exp(-z/L), (28)
where the relaxation length L is given by It is important to emphasise the status of this length L.

L 2 = 1 + anI n ∞ anI n ∞ ℓ 2 = ℓ 2 Y n(Y -1) for Y > 1. ( 29
)
As Y is the control parameter of this particular thought experiment (or numerical simulation [37]), L can be expressed as a function of Y. It does not make Y a state variable which would control another state variable L. Note also that L can equally be expressed as a function of I ∞ .

Consider now the KEP equation ( 19) with the fluidity f proposed by Kamrin et al. (25). As above, making profit that Y is homogeneous, this equation can be linearised around f ∞ as

ℓ 2 ∇ 2 (f -f ∞ ) = ℓ 2 L 2 (f -f ∞ ), ( 30
)
where L is now given by

L 2 = ℓ 2 1 Y + 1 n(Y -1) for Y > 1. (31)
In the limit Y → 1, the relaxation length takes exactly the same form L ∼ ℓ/ n(Y -1) for the two models, despite their differences.

The second case corresponds to Y < 1, so that I ∞ = 0. Again, the equations can be expanded around I = I ∞ but this linearization is completely different from the previous one, as I = 0 is not solution of the equations: one needs the non-local term to get a solution. With the KEP model (19), one obtains

ℓ 2 ∇ 2 f = (1 -Y)f, (32)
which gives exponential relaxations over a length L given by

L 2 = ℓ 2 (1 -Y) for Y < 1. ( 33
)
Using the gradient expansion model ( 13), we get

ℓ 2 ∇ 2 I = χ -1 (Y)I. ( 34
)
The shear rate therefore relaxes over a length L given by Fig. 7. Function χ(κ) numerically measured for Y < 1 for frictional (red circles) and frictionless (yellow circles) grains -the system is that described [37]. The black solid line is a fit with the empirical expression χ = √

L 2 = ℓ 2 χ -1 (Y) for Y < 1. ( 35
)
(1-κα) 2 +κβ(κα-1) 1-ακ with α = -15.95 and β = 16. 30. Inset: zoom on the small values of κ, in Lin-Lin axes.

In the vicinity of the critical conditions, κ is indeed small, so that the linear approximation χ(κ) ≃ 1-κ can be used. The divergence of L at Y → 1 is therefore given, again, by eq. (33).

As an illustration, we display in fig. 6 such diverging relaxation lengths extracted from numerical simulations of sheared layer [37]. Velocity data, such as presented in fig. 1b can be obtained systematically for Y values above and below Y = 1. The fit of the velocity profiles is made with a function of the form γb z + C sinh(z/L), where C and L are adjustable. This provides a direct measurement of a relaxation length L, which effectively diverges on both sides of the critical point Y = 1 according to the theoretical predictions (29) and (35). Figure 7 shows the shape of the function χ(κ). Its non-linear behaviour, which roughly starts when κ > 0.1, is at the origin of the asymmetry of L with respect to the yield point |Y| = 1 when sufficiently far away from this point.

In conclusion, as far as linearisation is concerned, the two models lead, despite their distinct starting point, exactly to the same predictions for constant stress conditions in the vicinity of the yield condition Y = 1. Because two equations giving the same exponential solutions are not necessarily equivalent, this shear cell configuration in which the yield parameter Y is controlled and homogeneous thus cannot be used to discriminate between the different possibilities to build a fluidity. Further tests focusing on time transients and on heterogeneous situations are needed to test the starting constitutive equations, which do not reduce to their linearised expressions.

Does the KEP rheology reflect a dynamical phase transition?

An important claim made by Bocquet et al. in [6] is that the KEP constitutive relation reflects a dynamic phase transition controlled by the stress, in relation to the divergence of the relaxation length L on both sides of Y = 1. Because the model used by Kamrin et al. [55] is derived from the KEP approach, this claim would also apply to granular matter. On the opposite, in our framework, we argue that the non-local rheology describes the same liquid phase above and below the yield conditions. Because this controversy concerns an essential point of physics, we find it important to develop in this subsection some technical but essential details of this issue.

As already mentioned, the starting KEP equation ( 19), when linearised around a homogeneous stress state corresponding to a given constant Y, leads to the generic equation:

∇ 2 f = f -f ∞ L 2 . ( 36
)
f ∞ and L depend on the value of Y. This is the Ginsburg-Landau equation used by Kamrin et al. [55] with the fluidity given by eq. ( 25). A crucial point is that this equation is used by these authors in non-homogeneous situations, i.e. taking for Y the local value, with functions f ∞ (Y) and L(Y) [68,71]. This would be a correct assumption for a slowly varying stress field if Y was a state variable. We have discussed above why it is not the case, even though it is the control parameter of the considered linearisation.

Let us further illustrate the mathematical differences between two seemingly equivalent derivations. Consider the KEP constitutive relation (19) with γ = f Y (eq. ( 25)). Associated with (22), which determines I(f ) in the homogeneous case, and whose solution is given by (24) for n = 1, we obtain

ℓ 2 ∇ 2 f = f 1 1 -aT f -Y . (37)
However, let us alternatively start from (36) and plug in expressions of f ∞ and L. f ∞ is the solution of the equa-

tion γ = T -1 I(f ∞ ), i.e. f ∞ (Y) = (Y -1)/(aT Y).
The expression for L depends on whether Y is larger (eq. ( 31)) or smaller (eq. ( 33)) than unity. To make it compact, let us focus close to the yield condition Y = 1 for which both cases can be summed up with L(Y) 2 ∼ ℓ 2 /|Y -1|. Doing so, instead of (37), we obtain

ℓ 2 ∇ 2 f = |Y -1| f - Y -1 aT Y . ( 38
)
Equations ( 37) and ( 38) are obviously not the same. In particular, their behaviour is very different when Y → 1 as the right-hand side vanishes in the second case but stays finite in the first one. In other words, the original KEP equation ( 19) and the final Ginsburg-Landau equation used by Kamrin et al. are irreducible one to the other: the transformation is neither mathematically nor physically justified.

The divergence of the relaxation length L on both sides Y = 1 has been interpreted as the signature of a dynamic phase transition controlled by Y. The region Y < 1 would correspond to the solid-like behaviour while the region Y > 1 would correspond to the liquid-like behaviour. However, the original KEP model (37) was entirely based on the description of a single liquid-like phase and the Ginsburg-Landau equation (38) used by Kamrin et al. [55] is not a controlled approximation of it. Y, which is not a state variable, cannot be the parameter controlling the phase transition. The liquid flowing phase can exist even below the yielding conditions, for Y < 1.

Fluidity and boundary conditions

An important consequence of the choice of a particular fluidity f is the underlying assumption that f and its gradient ∇f are continuous, otherwise the use of the Laplacian operator would not make any sense. This remark provides a constraint on the nature of f , which can for example be tested in a situation where the stress varies extremely rapidly in space between two states. In this spirit, we have performed numerical simulations where a secondary micro-rheometer is placed in the bulk of a shear cell as presented previously (see figs.1b and 2 and details in Bouzid et al. [START_REF] Bouzid | Andreotti[END_REF]). Shearing within the micro-rheometer is obtained by means of localised bulk forces along two lines which induce a discontinuity of the shear stress. We have measured numerically the ratio R of the absolute value of the shear rate on one side and on the other side of the stress discontinuity. The pressure remains constant and the direction of shearing is reversed ( γ changes sign at the discontinuity). Figure 8 shows that | γ| is indeed continuous (R = 1) in both frictionless and frictional cases. On the opposite, the fluidity proposed by Kamrin et al. (eq. ( 25)) is not in agreement with the data, especially when the yield parameter approaches zero.

0 -0.2 -0.4 -0.6 -0.8 -1 1 0.1
Fig. 8. Ratio R of the absolute value of the shear rate on the outer and inner sides of a stress discontinuity, as a function of the yield parameter Y m. Yellow circles are for frictionless grains and gray lozenges are for frictional grains. Thick solid line R = 1: prediction of the gradient expansion model [37] using I as fluidity. Dotted line: prediction of fluidity theory [55] using f = γ/σ as a fluidity parameter.

Further tests: the question of dynamical mechanisms

We have presented a short review of non-locality in granular flows. We have mainly focussed on the comparison of two models: the KEP model adapted for granular matter by Kamrin et al. [55,68] (see eq. ( 25)), and the gradient expansion that we have proposed [37] (see eqs. (12) and ( 13)). The difference between these models has not been recognised so far in the literature, essentially because they lead to the same predictions for the velocity profile in a situation where the stress is homogeneous. However, our conclusion is that these models are fundamentally different and that their differences can be tested.

These tests must be performed in situations that are strongly heterogeneous in space, like the one shown in fig. 8, or to unsteady situations [75]. They can of course concern the direct predictions of the model, but also the choice for the fluidity, which is not necessarily the inverse of viscosity, as well as the associated boundary conditions. Importantly, these conditions should not be fitted, but part of the physical analysis of the problem. More fundamentally, the difference between the models is to be found in the hypothesis made to derive them and in particular in the dynamical mechanisms underlying their dynamics. For instance, many different explanations have been proposed for the very same Herschel-Bulkley rheology. In the context of non-locality, let us give several examples where ingredients could be tested. Different models assume the proportionality between the decay rate of fluidity and the rate of plastic events. Such a proportionality can therefore be investigated experimentally. Other models like KEP prescribe not only the average fluidity but also its distribution. The measurement of such a quantity is a more severe test than the fit of velocity profiles. Would a model assume the existence of microscopic yield conditions for the nucleation of plastic events, it would then be necessary to determine this quantity, to show that it exists and that it is constant.

In the case of granular matter, we have shown that the main point separating the KEP model and the gradient expansion model, is the existence or not of elasto-plastic localised events in the liquid regime. The KEP model is directly adapted from soft matter and assumes that elasticity dominates the dynamics. In the test presented here, we have numerically shown that, in the rigid limit, there are no localised plastic events and the flow is dominated by non-affine collective motion along soft modes. One could argue that the Coulomb friction condition at the contacts between the grains may lead to plastic events. However, comparing frictional and frictionless grains, we do not see any difference neither on non-locality nor on the absence of plastic events. Beyond other important reasons, we have shown that the fluidity proposed by Kamrin et al. [55], as an extension of the KEP model to rigid granular packings, is not a state variable and is thus not continuous across a stress discontinuity. We acknowledge that, in spite of the fact that our choice of the fluidity parameter for dense granular flows respects the state variable requirements and quantitatively predicts some situations, it does not clarify the understanding of the actual microscopic mechanisms at work to definitely unravel the question of non-local rheology. We recently discussed some limitations on such a choice and we pursue the work of identifying the microscopic or mesoscopic processes associated with the flow of hard grains [START_REF] Bouzid | Andreotti[END_REF]75].

Finally, amongst the points that have created a confusion in the literature, is the fact that granular matter does present localised plastic events, but only in the solid regime [61,76], not in the liquid regime discussed here, were grain elasticity is irrelevant. Let us note that numerical simulations performed with the standard Coulomb model of friction at contacts, which perfectly reproduce observations in the dense liquid regime and in particular non-locality, are not able to reproduce creep in the solid regime. These two regimes (solid-like and liquid-like) must eventually be described, but the transition between these dynamical phases is known to be subcritical and to present a hysteresis, a key aspect of granular matter that remains unexplained at present. In this context, a nonlocal transition between solid and liquid was addressed by Wyart [77], based on the generic outcome of the Maxwell rigidity transition for hard granular packing. Importantly, the KEP model claims to describe this dynamical phase transition as a critical transition controlled by the stress, and the rheology both above and below the transition. The gradient expansion, on the opposite, is based on the fact that shear stress cannot be a control parameter for this transition and describes the system as a unique continuous liquid phase both above and below yield conditions. This approach is therefore perfectly compatible with a subcritical transition to the solid regime, as it does not describe the later.

MOTIVATIONS

Motivations

As explained in [8], confined geometries (such as the heap or the narrow channel flow) present certain interests as they exhibit several granular behaviors which allow one to distinguish regions in the system: at the free surface, a so-called "gazeous", almost collisionless regime where the particles have ballistic trajectories ; a flowing ("liquid"?) regime where the flow is dominated by the shear rate and a quasi-static zone also mentioned as "creeping zone". Authors of [8] report two distinct regimes in the quasistatic region: the upper one is dominated by the evolution of the volume fraction and wall friction whereas the region beneath it is referred as a glassy zone, where grains fluctuate within successive cages. On the one hand, as mentioned by previous studies [9,10,11,12], it enables to exhibit the role of wall friction and the effect of wall confinement. On the incline plane (with no walls, or at least very distant from each other [13]), the stress state is theoretically homogeneous along the whole depth of the medium. The addition of the shear due to the wall friction tends to slow down the flow in the depth, as the pressure on the walls -and therefore the tangential component of their resisting force, caused by friction-increases. This leads to the crossing of a theoretical flow limit initially predicted by the simple study of the stress configuration in the system. The existence of a "creep zone" below the flowing limit has been observed by Komatsu et al. [14]. In this chapter, we will present the experimental set-up of the inclined narrow channel we built in order to observe the granular flow in the whole range of depths, including below the theoretical limit mentioned above. We will present the flow predicted by the rheology usually used to model dense granular flows and exhibit its shortcomings, mainly due to the fact that it lies on a local description of the flow.

The latter justifies the need to implement a correction in the previous rheology in order to take into account the influence of the non-local effects in this system. We will then bring insights to experimentally investigate this system in its entirety.

The experiment 2.2.1 Experimental set-up

The set-up we use is an inclined narrow channel within which we monitor the flow with respect to the depth. The control parameters are the flow rate Q and the dimensionless width of the channel W * = W d (W and d being respectively the width of the channel and the diameter of the beads). This system is interesting because unlike the incline plane for which the yield parameter Y is fixed for the whole medium (on the whole height of grains we have Y = tan θ µc ), the (frictional) walls add a negative shear on the flow that continuously decreases Y along the depth. As a result, going from the flowing regime at the free surface to the so-called "creep" in the depth, we can investigate the evolution of the yield parameter and the crossing of the yield stress σ y . Unlike the first 2D numerical experiment conducted in [15] and Beads are all of the same kind (same reference from Silibeads) but in two different colors in order to increase black and white contrasts on the images. [16] this experiment is therefore an important step in probing the non-local rheology we propose, in a geometry within which the yield parameter Y is not homogeneous.

The reader should note that at the output of the channel, a patch imposes a certain height of grains even when there is no flow in the channel. This imposes a static zone in the bulk, typically below the flowing limit predicted by SSH [9] and where non-local effects are therefore dominant.

The granular media: For all our experiments, we use coated glass beads from Silibeads 1 of average diameter 1.15mm (with a factory polydispersity value of ±0.15mm which avoids crystallization). The coating provides better imaging conditions. We first tried non-coated glass beads but issues with the lighting and multiple reflexions within each bead 2 lead us to change the material. Also, we used yellow fast cam lighting spots we tinted in yellow. The yellow light, with the mixing of blue and yellow coated beads enable better contrasts on the image and thus more signal (see photo 2.3). The lighting comes from beneath (see picture) and from the sides. Ideally, it should be homogeneous and therefore come from all around the lens, but we had reflections on the glassy walls, hence this oblique and diffuse lightning.

All the experiments where conducted in relatively stable but not controlled atmospheric conditions: temperatures around 20 -25 • C and humidity of 30 -35%.

This might be important as we know that dry grains attract humidity from ambient air and therefore capillary bridges can form [7] between the grains and affect the cohesion of the medium, and thus the rheology we try to characterize. 2 Reflexion spots within the bead varies with the position of the particle with respect to the lamps: particles scintillates across the images and their individual images thus evolves during their passage in front of the lens 62

THE EXPERIMENT

Control parameter: the flow rate

Flow regulation is ensure by the use of two reservoirs, one feeding the other. In 1895, Janssen measured the pressure at the bottom of a silo filled with grains (see original article translated by Sperl [17]). He observed that the pressure saturates when more and more grains were added into the silo. In other words, the vertical stress in the silo saturates above a certain height of particles in the reservoir because of the force chains within the medium: the weight of the medium is supported by the walls of the container. This saturation of the stress above a certain height from the aperture of a reservoir explains the independence of the flow rate with respect to the high of grains in the reservoir, above a critical height, in contrast of what would happen with a liquid. We use this principle as a way to regulate the flow poured into the channel (see figure 2.4): the main reservoir feeds the intermediate reservoir ("buffer"). In the latter, the level of grain is kept constant at the level of the outlet of the tube3 by which grains are released from the main reservoir.

Regulation of the flow rate is controlled by changing the outlet diameter of the buffer reservoir (see figure 2.4) and is monitored by a weighting scale linked to a computer that continuously record the output weight. We therefore compute a mass flow rate. This is important for several reasons: first it enables to make sure the flow is constant through the channel. Also, controlling the flow rate by changing the outlet diameter guaranties to have a relatively reproducible experiment. This will be important when we will try to match the velocity profiles measured from both the top and the side of the channel in order to have an estimate of the 3D structure of the flow and monitor the flow from different side, at a given flowrate. Also, note that we recover the Beverloo law4 [18] with our bin (see figure 2.5). This enables

THE EXPERIMENT

Data measurements

Fast camera imaging To monitor the flow we record images with a Phantom M340 fast-camera 5 (mounted with a Nikkor 60mm lens of f = 2.8) settled on the side of the channel and which optical axis is perpendicular to the channel glass walls (see figure 2.2). The reader should therefore note that we only have access to the beads located against the walls of the channel, not in the bulk. This will be important in future developments, when comparing experimental data with the theory (see section 2.4 starting page 69).

Angle of the free surface We built the channel in a way that the tilt angle is adjustable according to the imposed flow rate, so that the flow is parallel to the rough bottom: we seek to create a stationary flow, homogeneous along the channel direction (direction x, see figure 2.2 page 60). Indeed the angle of the free surface with the horizontal axis is selected by the system itself (see section 2.5.3 and especially equation 2.32). We therefore measure this angle θ.

How to run the experiment

Set the fast camera We had issues in defining the free surface on our first datasets for which the camera axis was precisely set perpendicular to the walls and horizontal: in order to monitor the flow from the free surface to the bottom on a single acquisition, we use to set the camera axis a bit too high and we hence also monitor some grains at the free surface that were not in contact with the wall but rather at the center of the channel. We therefore adjusted the position and the orientation of the camera. First, we lowered the camera a little, and compensated this by very slightly tilting the optical axis (∼ 3 • ) in order to keep the same window on the flow, from the free surface to the bottom of the channel. Second, we increased the aperture of the lens at maximum. Doing so deteriorates the depth of field so we are sure to only focus on a single layer of grains, the one against the wall. Increasing the lens aperture also enables to have more light, hence more signal providing no pixel saturation on the sensor (we adjust the aperture speed accordingly).

Initiate the flow and set the channel inclination so that the free surface is parallel to the bottom of the channel. This ensures that we have an homogeneous flow along the axis of the channel.

Wait until the flow rate is steady. We monitor the mass of grains at the output of the narrow channel and plot it with respect to time. We can therefore see when the flow is steady.

These two steps mix as they appears to be quite dependent: when changing the tilt angle of the channel, it momentarily changes the output flow and therefore the curve mass vs time which we use to consider whether or not the flow is stationary.

Also, especially for extreme flow rates, the metastable behavior of granular media is exhibited: low flow rate means small angle θ and therefore the tilt angle gets closer to the angle of arrest. The flow hence exhibits a bistable behavior, where avalanches start from somewhere along the channel, propagate and stop. Any vibration (sound, shock on the setup, etc.) can be the opportunity for the flow to start, without any guaranty on its time life nor on its stationarity.

Calibrations of the images lies on two main measures: the calibration of the length scale and the precise measure of the orientation of the flow with respect to gravity. For the latter we use a plumb line and take an picture of it (actually, the average image of a film so we can get rid of the remnant oscillations of the pendulum). This gives the angle the camera sensor makes with the vertical axis. We can then deduce the real angle of the free surface by adding this measure to the tilt angle of the images. As we said, this measure is directly linked to the flow rate (see relation 2.32). It is therefore crucial to measure parameters of the rheology.

THE EXPERIMENT

Calibration of the length scale is made using a millimetric paper sheet stuck at the wall, on the inner surface of the side we monitor. For each flow measurement we took a picture of this millimetric grid in order to provide a scaling reference. We can therefore convert our data into SI base units.

Strategy of measurement of the velocity profiles

We want to measure the flow velocity along the main direction of the flow (i.e. along the x-axis). As we tilted the images so that the free surface (and thus the means flow) is along x, we built an algorithm to measure the velocity field with respect to the depth, i.e. with respect to each line of the image.

The reader should note that we therefore measure the velocity of the grains in contact with the wall, not in the bulk of the channel. Also, it is important to emphasize the motivation and thus the choice of such a strategy of measurement:

we suppose the flow to be steady, and invariant along its dominant direction x. The algorithm we choose to build is therefore orientated toward the measure of the mean velocity profile as we consider that all the particle at a given depth have the velocity.

We therefore neglect their velocity fluctuations (in any direction).

The core of the algorithm lies on the comparison of two horizontal lines of pixels that correspond to two different images -and thus instants -of the film. To do so we use the χ 2 test6 By adjusting the best match between line n at time t and line n at time t + ∆t shifted by a known integer number of pixels, we manage to compute a sub-pixel estimated best shift of the studied line between the instant of image i and instant of image j that follows in the film and which delay with respect to image i is known. From this we easily compute the translational velocity of this pixels line.

A schematic of the principle of this algorithm is given figure 2. 6.

As we post-treated the images so that the mean flow is horizontal on the images, Figure 2.6: Principle of the velocity profile algorithm. We compare the same pixel line at two different time (i.e. same depth/line on the image but at t and t + dt) and look for the dt that provides the best match (in this example a shift of dn = 3 images, so 3 times the opposite of the frame rate of the camera) with the pixel line at time t. We hence deduce the velocity at which this line moved between the two images.

we can assume this velocity is the average velocity of the grains at this depth, along the horizontal axis on the image (which corresponds to the main axis of the channel, axis x in the model, see section 2.4).

Post-treatment

The whole process for the data acquisition and the post-treatment can be summarized as follow: Once we have the raw data of the velocity flow rate, we need to localize the free surface in order to shift all the velocity profiles (with respect to each flow rate), so the all share the same origin (with respect to the depth, z-axis). To do so, we compute an average intensity profile of the tilted average image of the film (average image which we computed to measure the tilt angle of the free surface with respect to the camera sensor, see figure 2.7). Fitting with a phenomenological function such as i ntensity (x) = A + B tanh (C (x -D)) provides a good estimation of the localization of the inflexion point of the intensity profile. We choose to consider this inflexion point as the position of the free surface on the image. Figure 2.8 shows raw velocity profiles we obtain from the algorithm 2.2.2 and the position of the free surface velocity. We can then shift the z-axis of all the velocity profiles so that they all exhibit the same free surface z = 0 and convert the data to the International System of Units.

Hypothesis

We consider the incline narrow channel set-up presented figure 2.9 (page 70). As the dominant direction of the flow is along the axis of the channel, we choose a reference frame accordingly: x is the axis of the channel, z toward the depth and y in the direction of the width of the channel. Origin of the depth ( z axis) is at free surface and y = 0 corresponds to the middle of the channel.

The flow is supposed to be steady and invariant along x (we assume that the flow we consider is far enough from the injection in the channel).

The system, from the top view (in the ( x, y ) plane) is the classical situation of a parallel flow. As usual in such situation, we assume the flow to be unidirectional, along axis x. The velocity therefore writes:

v = v (y , z) x

CONTINUOUS ANALYTICAL MODEL FOR DENSE GRANULAR AVALANCHES

We consider the medium as a viscous isotropic fluid of viscosity η(y , z), a function a priori only of y and z (as the flow is permanent in invariant along the x direction).

The flow is supposed to be incompressible (i.e. the volume fraction φ is supposed to be homogeneous in the whole medium 7 ). The reader should note that the φ = cst implies the fact that the medium density is as well constant in the medium, as linked to the material grain density ρ p by ρ = ρ p φ. It is important to mention that some previous studies report [8,19] a dependency of the profile of volume fraction on the tilt angle θ, with respect to depth. In the following, the main effort are devoted to the description of the whole flow, especially the "creep regime" in the depth. As φ saturates relatively quickly in depth at φ = 0.6, we consider in the following -as a

first assumption-that it is constant and at its saturation value (φ = 0.6 is commonly assumed for dense granular packing at rest and with such polydispersity).

Reservoir

Figure 2.9: Schematic of the narrow channel set-up 7 This hypothesis of the homogeneity of the volume fraction in the whole granular flow can be controversial as it decreases a lot in the region near the free surface, also where the pressure tends to zero and the shear rate is important. We will discuss this point later in the development of this theoretical computation.

Momentum balance

The momentum equation easily gives:

∇ • σ + ρ • g = 0 (2.1)
from which we obtain the three equations:

                 ∂σ xy ∂y + ∂σ xz ∂z + ρg sin θ = 0 (x)
∂σ y y ∂y + ∂σ y z ∂z = 0 (y )

∂σ zz ∂z + ∂σ zy ∂y + ρg cos θ = 0 (z) (2.2)
We use the 3D formulation

σ ij = -P δ ij + τ ij = -P δ ij + η γij (2.3)
where γij = ∂ i u j + ∂ j u i is the shear rate and η = τ / γ (τ been the shear stress). We draw the reader 's attention to the isotropicity of the normal stresses we make the assumption of, here. Indeed, we identify the "pressure" P as σ ii = P (z), such as we would naturally do with a Newtonian fluid.

So, we have:

     σ xy = σ y x = η∂ y u x σ xz = σ zx = η∂ z u x (2.4)
And

σ zz = -P (z)
All this, in the system (2.2) gives :

                   ∂ ∂y σxy [η∂ y u x ] + ∂σ xz ∂z + ρg sin θ = 0 (x)
∂σ y y ∂y = 0 (y )

- ∂P ∂z + ρg cos θ = 0 (z) (2.5)

CONTINUOUS ANALYTICAL MODEL FOR DENSE GRANULAR AVALANCHES

The momentum equation in (z) gives, after integrating 8 :

P (z) = ρgz cos θ (2.6)
In equation ( 2.2) we can integrate the first and the third terms, but not the second one. In order to estimate this velocity profile in the width of the channel, we propose to model the flow by a second order polynomial function:

u x (z) a(z) + b(z) y 2 (2.7)
and define the expression of the a and b functions. From equation (2.7) we have:

∂u x ∂y 2 b(z) y (2.8)
Boundary conditions: at the walls (y = 0 and y = W ) The friction at the walls is modeled as:

     σ xy (y = W/2) = µ w σ y y σ xy (y = -W/2) = -µ w σ y y (2.9)
Into first equation of 2.4 we therefore get (at the walls):

η ∂ z u x = -2 η(z) b(z) W 2 = µ w P (z) (2.10)
and putting everything together into (x) we get:

b(z) = - µ w P (z) η(z) W (2.11)
and therefore:

u x a(z) - µ w P (z) η(z) W y 2
(2.12) 8 As we consider a dry granular medium and thus the fact that air can be also find between the grains, so that on the whole system air pressure can be assumed to be invariant, we can consider P (z = 0) = 0 at the free surface and neglect the influence of air and atmospheric pressure.

Using ∂ y σ xy = ∂ y (η∂ y u x ) = -2µw P (z) W we get:

ρg sin θ - 2 µ w z cos θ W = - ∂ ∂z [σ xz ] (2.13)
which integrates into:

σ xz = ρg µ w z 2 cos θ W -z sin θ (2.14)
The reader should notice that surprisingly, the shear stress, σ xz , so as the pressure, does not depend on y .

As a consequence, the yield parameter does not either depends on y and is only a function of z. It writes:

Y = -σ xz µ c P = tan θ µ c - µ w µ c W z (2.15)

Local rheology approximation

Following the seminal publication by GDR Midi [20] and previous studies ( [21]), we first present the solution the local rheology would provide.

Following [15] and because we assume (and checked a posteriori ) to be in the same range of stresses and deformations, we recall the expression of the local rheology for dense granular flow, in the frictional 9 case:

µ = µ(I) = µ c (1 + a I) , (2.16)
which links the stress state to the dimensionless shear rate I = | γ| d P/ρ p (d being the typical grain diameter in the medium). 9 For such range of solicitations (I between 10 -4 and 10 -1 ), the data are perfectly described by a law of the form:

µ(I) = µ c (1 + aI α )
(the residuals form a statistical noise ; error bars are ≈ 5%). In the frictional case, it has been empirically found that α 1 ; in the frictionless case, α 0.5. See also [22,23].

LOCAL RHEOLOGY APPROXIMATION

Since:

σ xz (z) = µ(I)P (z) (2.17)

The dominant contribution of I only depends on z. Thus we take | γ| = -∂ z u x and neglect the ∂ y u x contribution as by symmetry in the channel, at the center y = 0 it can be considered as a marginal term 10 . The dimensionless equation is therefore:

tan θ µ c -z = (1 + aI) , (2.18)
where the relevant length scale is

L = W µ c µ w
. Therefore I reads:

I = 1 a tan θ µ c -1 -z , (2.19) i.e.
I = 1 a (Y -1) , (2.20)

Flow thickness

I vanishes at a depth h = tan θ µ c -1 (2.21)
which coincides with the Coulomb criterion, allowing one to rewrite:

I(z) = 1 a h -z . (2.22)
Note that dimensionally, this gives:

h = W µ c µ w tan θ µ c -1 (2.23)
10 However, this approximation can be tested a posteriori from the profile v (y , z) hence obtained.

and

I(z) = 1 a µ c µ w (h -z) W .
(2.24)

Velocity profile

The velocity profile u(z) is related to I by:

I = du dz d √ g φ z cos θ . (2.25)
We define the dimensionless velocity as:

V = L 3/2 √ g φ cos θ d ,
such that:

I = z-1/2 d ũ d z . (2.26)
The force balance equation integrates into:

ũ = 1 a 2 3 h h3/2 -z3/2 - 2 5 h5/2 -z5/2 . ( 2.27)
The reader should note that this is expression of the velocity should only be valid in the plane corresponding to the middle of the channel (y = 0). Here, ũ is in fact ũ(y = 0, z). It coincides with the dimensional expression used to fit the data and extract estimates of the parameters:

u(y = 0, z) = V a L 5/2 2 3 h h 3/2 -z 3/2 - 2 5 h 5/2 -z 5/2 . ( 2

.28)

For the sake of simplicity, in the following, we choose to neglect the dependency of the velocity profile with respect to y . Doing so, we consider that the velocity profile along z is the same in the whole width of the channel (so-called "plug flow").

However, the reader should keep in mind that the 3-dimensional profile would be 2. 5. LOCAL RHEOLOGY APPROXIMATION (using 2.12 and the relation between the viscosity and the dynamic friction coefficient 11 ):

u(y , z) = u(0, z) -y 2 µ w W d I µ(I) P (z)/ρ p , (2.29)
where I is given by 2.24, so:

u(y , z) = u(0, z) -y 2 µ w W d 1 a 2 P (z)/ρ p 1 µ c + W µ w z -W tan θ .
(2.30)

Flow rate

As said previously, we consider the plug flow profile (2.27). From the integration of this velocity profile over the width of the channel and the flow height, the flow rate writes:

Q = 4 35 V W aL 5/2 h 7/2 = 4 35 V LW a h7/2 . (2.31)
Note that we get from the weighting scale the mass flow rate is ρ p φQ (where Q is the volumic flow rate). Using the expression of h, we get a relation between the flow rate and the slope of the free surface:

tan θ = µ c + B Q Q * 2/7 , (2.32)
where Q * is a characteristic flow rate:

Q * = W 7/2 d √ g , (2.33)
and:

B = µ c µ w µ c 5/7 35 4 a √ gφ cos θ 2/7
.

(2.34)

In the following, we use these relations to fit our data and extract values of the parameters of the rheology. We can then try to fit the velocity profiles with the local rheology (see figure 2.11).

In this process, the flowing height h is at the same time a fit parameter and a cursor that sets the range of depth over which we compute the fit. We therefore iterate the fit of the velocity profile in order to converge on an auto-coherent value of h. From the fit of the velocity profile with the local rheology (figure 2.11 and formula 2.27) we get an estimation of a = 2.59 ± 0.06.

Conclusion on the local rheology

The plug flow hypothesis seems quite questionable: we measure velocity profiles at the free surface (the flow is monitored from the top of the channel, in the ( x, y ) plane) that exhibit a highly three dimensional behavior (see figure 2.14). The estimation of the flow rate and the parameter values we made from it is therefore approximative as it under-estimates the flow of granular matter in the channel.

From these top view velocity profiles, we can estimate error made on the flow rate: ∆ Q Q is the relative error 13 made by computing the flow rate from the integration of the side velocity profiles and using the plug flow assumption. In this expression, Q represents the dimensionless flow rate measured from the weighting scale. We see on figure 2.13 that a nearly 40% of the flow is forgotten by the flat profile model.

Also, it is important to notice that the flow velocity is not zero under the flowing layer (of height h) contrary to what predicted by the local rheology. Our detection algorithm allows to measure velocities over nearly six order of magnitude, as seen on figure 2. 15. Clearly, one needs to take into account the non-local feature of the system. In the next section, we propose to probe the non-local model initially proposed by Bouzid et al [15] and tested in 2D numerical systems, within homogeneous stress configurations 14 . Doing so, we aim at describing the flow on the total measurable 14 Here, the term "homogeneous" refers to the spatial invariance of stress state in the bulk region, even if it is below the yield stress, Y < 1. Of course, this flow can only exist if there is a flowing regime (Y > 1) somewhere in the system, therefore the Y is not literally homogeneous in the whole system.

13 ∆ Q = Q -W Q * total height u(z)dz

Non-local rheology

As previously said, one of the main interest of the narrow channel set-up is the coexistence of the flow above and below the flowing limit predicted by the local rheology (see formula 2.21 or its dimensional form 2.23). Here, we therefore propose to find an expression to describe the whole flow. We use the non-local rheology proposed by [15] to model the non-local effects by taking into account the fluidity environment. In this section, we propose to study the flow in the narrow channel in the frame of this non-local model.

Non-local equation

We use the non-local rheology expression from [24]:

-σ xz = µ(I)χ (κ) P (2.35)
where κ =

Y = µ(I) µ c χ(κ), with χ(κ) 1 -κ + O(κ 2 ) (2.36)
The rheology therefore reads:

Y = - σ xz µ c P = (1 + a I) 1 - 2 L 2 I d 2 I d z2 (2.37)
From which we get (from the expression of the yield parameter 2.15) the non-local equation:

d 2 I d z2 = 1 - 1 + h - z 1 + a I L 2 I 2 (2.38)
The reader should therefore note that, as the yield parameter 2.15 only depends on z, the non-local equation is an ordinary differential equation on z only. According 2. 6. NON-LOCAL RHEOLOGY to our model, the flow in the depth is therefore independent of y . This equation (2.38) is non-linear and should therefore be solved numerically. We use a built-in Matlab Runge-Kutta method called ode45 15 to do so. This resolution needs to be initialized: we need to start the resolution from a finite depth zstart at which we know the value of I and its derivative dI d z . It is clear that the flow at the surface is unknown 16 . We therefore choose to initiate the numerical resolution from an arbitrary large depth (large with respect to the value of the supposedly flowing height h the local rheology predicts), at which we compute the value of the inertial number I and dI d z using a reduction of the non-local equation in the case of an infinitely small value of I and using an asymptotic estimation the corresponding solution I asy mpt its derivative dIasympt d z when z → +∞.

Asymptotic expansion

Boundary condition: in the depth, asymptotically when a I 1 (i.e. when z → +∞), equation 2.38 becomes:

d 2 I d z2 = I L 2 z -h (2.39)
As I must vanish in depth, one gets the asymptotic z → +∞:

I C Air y Ai L 2/3 z -h , (2.40)
where Ai r y Ai is the Airy function of the first kind. The expansion of the Airy function gives:

15 ode45 is a very widely used Matlab solver. It is adapted for most classical (smooth, continuous) problems. It is a Runge-Kutta method of order five (use of the Dormand-Prince integration scheme), the error being given by the fourth-order estimate. The reader is referred to the dedicated online MathWorks documentation. 16 we recall that we used a rather arbitrary argument -on the intensity profile of our images-to adjust the z = 0 of the velocity profiles. Also we noticed the free surface is not described by the local rheology as seen on figure 2.11 page 78. For all these reasons, no boundary condition at the free surface can be given up to now.

I C 2 √ π L 1/6 z-1/4 exp - 2 3 L z -h 3/2 (2.41)
where C is a constant that should be defined by a boundary condition to the system.

The expansion of the derivative of the Air y function gives:

dI d z - C 2 √ π L -5/6 z1/4 exp - 2 3 L z -h 3/2 (2.42)
And finally, the integration of expansion 2.41 provides the corresponding velocity. It reads:

ũ C 2 √ π L 7/6 z-1/4 exp - 2 3 L z -h 3/2 (2.43)
From which we get the dimensional expression of the asymptotic velocity:

u(z) C 2 √ π W µ c µ w 7/12 7/6 z -1/4 exp - 2 3 3/2 µ w W µ c (z -h) 3/2 (2.44)
In this formula, C may depend on the control parameter of the experiment (θ, or equivalently, the flow rate Q), but the parameters of the rheology (critical friction coefficient µ c and the non-local term ) are the same for all runs (i.e. any flow rate)

and should not depend on any control parameter of the experiment. Also, h = h(Q)

(or "h = h(θ)") so it is fixed and determined for each dataset 17 .

We write formula 2.44 by blocks, each of which is to be fit with respect to the experimental data:

u(z) = C 1 z -1/4 exp - 2 3 C 2 (z -h) 3/2 (2.45)
We tried to fit the tails of the velocity profiles all together, in an auto-coherent (or "optimized") way by sharing the value of the rheological parameters among all 17 This rises the question whether the fit with the local rheology can provide the measurement of h(θ) and any other parameter of the rheology. In deed, it could be tempting to fit the values of the (local) rheology and impose them in the fit of the non-local rheology. We argue here that this would be a mistake, as we know that that the local µ(I)-rheology cannot properly describe the system. Hence the values its fit would provide would obviously not match the non-local fit.

NON-LOCAL RHEOLOGY

the datasets, within the fit algorithm. Up to now, we did not manage to fit in an auto-coherent way all the data: fixing the value of h for each dataset using the values of h previously found (see fit figure 2.12 page 78) seem to constraint the fit too much and prevent its convergence.

In order to obtain self-consistent values for the fit parameter over all the experimental datasets, we follow a four steps procedure to fit the model (equation 2.45) on the data:

1. fit each dataset individually ; 2. average the value of C 2 and then fit again ; 3. get the value of h from the previous fit and fit h vs tan θ according to equation 2.23 (page 74). The slope is W/µ w and we hence obtain a good estimation of the friction coefficient at the wall µ w = 0.17 (which is realistic considering the values we can classically find in the literature for dry contacts). 4. average C 2 (again) and impose the theoretical value of h obtained from the previous fit (step 3) and then fit.

In this process, the values of h obtained from step 2 (fit of all the parameters except C 2 which is fixed at its mean value obtained from step 1) are of the order of the millimeters, which is consistent with the prediction of the local rheology. However, it appears that for the datasets at small flow rates (i.e. for reservoir outlet of diameter 10mm and smaller), the values of h computed by the fit algorithm are negatives, if not smaller than one bead diameter. This is physically inconsistent and therefore we cannot use these datasets to further compute the parameters of the rheology. This issues might be explained by the large ratio "signal" over "noise in the raw data measurements": the smaller the flow rate, the closer to the arrest condition we get.

The system becomes bistable, exhibiting random temporally distributed avalanches that onset and stop in an random intermittent manner, following stress fluctuations in the system. The origin of these latter is still controversial and is not the purpose of the current study.

We gather the fit values of the block parameters in the following table 2 From these values, we can provide an estimation of the parameters of the nonlocal rheology: µ c = 0.47; µ w = 0.17 which gives = 0.44 using the definition of C 2 (see equations 2.44 and 2.45). The values we find for µ c and µ w supports what we found in the previous section using the local, µ(I)-rheology, which also provided a value for the parameter a ( 2.6) 18 . This first experimental test of the non-local model also provides a value of the non-local parameter . It is found to be quite smaller than what computed in previous work [15] which numerically predicted (although for a 2D system) = 2.8.

Discussion

In this section, we presented an experiment in which the yield parameter continuously evolves spatially. We reviewed the technical tools we used to monitor the flow and measure velocity profiles over five order of magnitudes. We developed and probed the local rheology with our experimental data and gave some insights on the reasons of its lacks in describing the whole flow. This led us to probe the non-local model we presented in the previous chapter of this thesis, and to confront it with our experimental data, especially regarding the so-called "creep flow" regime. We showed encouraging results exhibiting the capacity of the non-local model to capture the flow beneath the so-called "flowing zone". As the non-local model entails the local µ(I) rheology, it is expected to be able to also predict the flow in the flowing zone as well.

As presented in the first section of this chapter, it has been shown [15,24] that the relaxation length of the velocity profile within the supposedly static zone is characterized by a symmetric dependence in the yield parameter Y, around the yield stress (i.e. Y = 1). In the present work, the fact that the non-local relation we developed seems to be able to capture both zones of our experiment seems to comfort this concept of a unique flowing zone.

To go further in the description of the flow in the narrow channel, a better understanding of the 3D structure of the flow is needed. Due to the high confinement and the influence of the walls, non-local effects clearly modify the flow also in the transverse direction of the channel, hence its 3D features. This system appears very complex because of the dominance of non-locality both in the depth and at the surface, and along the transverse dimension (width of the channel). One could therefore improve the theoretical model we used by implementing the dependence of the velocity profile in y (dimension of the width of the channel) such as presented earlier (see page 76 formula 2.30). Note that these velocity profiles from the top may also be fit using a cosh function (better fit than with the quadric we used to write the rheology, and still consistent with the latter at second order) such as shown As previously mentioned at the beginning of section 2.4 page 69, several hypothesis we made in the theoretical development may be questioned. First, we made the assumption of incompressibility, which has the direct consequence φ = cst (constant and uniform volume fraction). The discontinuity that the free surface of the avalanche flow seems quite obvious. This evolution of the volume fraction with respect to the depth has for instance been measured by [19]. They mainly observe a major evolution of φ in the so-called "flowing zone". The homogeneous volume fraction hypothesis may therefore be problematic if we want to continuously describe the flow from the depth to the flowing zone (and as close as possible to the free surface). However, with respect to this observation by, our development on the non-local rheology (with the asymptotic study) is valid as the volume fraction seems rather homogeneous in the depth ("creep regime").

Another hypothesis we made was by taking the linear form of the µ(I) rheology.

By doing so, we have to keep in mind that the theory we write is only valid in the limit of small values of the inertial number. We observed that the linear expression of µ(I) is phenomenologically valid up to I ∼ 10 -1 , which is basically the order of magnitude of the highest values of inertial number we estimated in our experiments, near the free surface (on grain size below the latter). Even if this assumption of small I seems questionable near the free surface, especially regarding the fact that at the free surface, the pressure vanishes, it may still be valid as the values of I reached in the region of the free surface seems rather small.

Also, because of steric effects, walls create important heterogeneities in the yield parameter profile Y(y , z). It would therefore be interesting to link the non-local behavior of this flow to the possible origins of non-locality. As we previously explained, systems of rigid particles under constant shear display a permanent collective motion, rather than localized plastic events that diffuse energy in the whole system during sudden cascades. However, it is suggested that the non-local behavior of dense granular media comes from an increase of the density of plastic rearrangements [25].

Some tools have been developed to experimentally measure the spatial distribution of local plastic events, such as Diffusion Wave Spectroscopy (DWS). In dense granular flows in channels, the non-local fluidity is therefore to be related to local plastic events which density can be monitored using DWS. It would for instance be interesting to link the nature of such "hot spots" (localized plastic rearrangements) to the non-local parameter .

An important parameter we did not vary yet is the width of the channel. Of course, the wider the channel, the flatter the profile in the transverse direction (( x, y ) plane). This would certainly provide more information in order to calibrate the friction coefficient with the walls and the influence of the confinement of the wall on the 3D structure of the flow (see [19] and references therein).

So as the width W , the imposed height of grains at rest in the channel (imposed by the patch at the outlet of the narrow channel) could also be varied in order to exhibit the influence of the bottom on the non-local features of the flow. Furthermore, removing this patch would allow one to probe the influence of the condition at a rough wall. Some experiments of the "Chute Flow" such as Pouliquen's [26] and recent DEM simulations [27] have been made but they test different models than 2.7. DISCUSSION ours. We are currently investigating this question of the definition of the boundary condition using our non-local rheology and 2D DEM simulations. This work is still under process. The whole logic to determine this theoretical issue is not so different from the one we will present in the following on the incline plane, which we use to find an a priori definition of the free surface boundary condition.

There are many other trails to investigate on this set-up: it could also be interesting to focus on the effect of the history of the preparation. In this experiment, between each run, we empty the channel, in order to renew the contacts and set the stress network in a repetitive way as it is then set by the imposed flow rate that fills in the channel.

It is also important to note that we realize after a little time that the coating of the beads erodes: this can affect the microscopic friction (friction between the beads), the stiffness of the contact forces, and also modify the asperities at the surface of the grains (see part III of the present manuscript).

Is there aging on the contact network? Does erosion of the particles have an influence on the aging of the whole medium?

For the sake of simplicity, in the following we leave these questions on the side and we focus on the definition of the boundary condition at the free surface. This question is crucial in order to be able to predict the flow profile in such systems.

Near the free surface, as a first model, we can assume the influence of the wall confinement to be negligible compared to the inertia of the moving particles, hence we choose to investigate the boundary condition by the mean of 2D Molecular Dynamics simulations on a very elementary system: the incline plane.

91

In the previous chapter we studied a dense granular flow in a 3D experimental setup. We have shown the importance of the three-dimensionality of the system.

We exhibited the lacks of the local rheology in continuously describing the flow, from the free surface to depth where we measured a creep motion of the granular media over a wide range of velocities. We proposed an asymptotic development of the non-local rheology in order to capture this creep flow in the depth. However this resolution seems incomplete as it does not yet allow one to solve the complete non-local equation on the whole domain.

In the following chapter, we investigate non-local effects in granular flows by the mean of numerical simulations based on Discrete Element Methods.

The plane shear cell and the incline plane are two elementary systems which exhibit a constant and homogeneous distribution of the yield parameter Y. In the aim of measuring the boundary conditions at the free surface using the incline plane set-up, we first probe our new numerical code on the shear-cell system previously used by [15].

As non-locality is exhibited near the jamming transition, we use a custom numerical system to calibrate the rheology in the vicinity of the transition. The system we use is a two-dimensional plane shear cell. What makes our set-up original is the possibility to impose the profile of the yield parameter Y = τ / (P µ c ), hence to impose the value of Y both above (Y > 1) and below (Y < 1) the jamming transition. The calibration of the non-local rheology is based on the measure of the relaxation length of the flow profiles, which diverges at the jamming transition.

We will then present another configuration which exhibits non-local effects. In the incline plane set-up, the stress configuration is such that the yield parameter is fixed by the inclination angle:

Y = tan θ µ c . (3.1)
It is therefore homogeneous on the whole height of grains, and independent of the value of this parameter. The local rheology thus predicts a constant profile for the 94 inertial number with respect to depth. This implies that the height of flowing grains does not have any influence on the flow. However experimental studies (see [20] and [28]) exhibit the existence of the hysteretic nature of granular flows: an initially static layer of grains of height h starts flowing when a critical tilt angle θ star t is reached. Once the flow is ongoing, the decrease of the inclination angle shows a second critical angle θ stop at which the flow stops. Reciprocally, these critical angles can be interpreted in terms of critical layer of thickness h star t (θ) and h stop (θ). The existence of such angle of arrest θ stop is a clear manifestation of non-locality.

By varying h for a given value of θ, we can therefore exhibit the non-local effects at a given value of the yield parameter Y, as it is only a function of the inclination θ.

In this study of the incline plane, we will therefore develop the non-local prediction and show how its prediction allows one to define the boundary condition at the free surface, where the non-locality is dominant. Measures of the relaxation length of the flow profile will be also compared to the calibration previously done on the plane shear cell.

Numerical simulations of granular media

Numerical simulations appear as complementary to experiments. They allow to probe a theoretical model and also to observe/measure some quantities that are not available in the experiments (for instance the number of contacts per grains in a granular medium). They also allow to study fictional experiments, by removing some physical interactions (for instance friction, of gravity). One of the advantages of numerical simulations is the possibility to vary continuously the values of some control parameters and to separate the different effects of the model parameters. This is very valuable in the understanding of the physics of material behavior, in which we try to understand macroscopic phenomena by exploring the medium at microscopic scales.

However, it is important to keep in mind that conclusions drawn from numerical methods are only valid in respect of the theory implemented in the code. These theories are based on simplifying hypothesis that characterized a conceptual representation of the medium. The comparison with the experiment is therefore only valid if the model is appropriate (coherent) with respect to the physical system.

A numerical simulation contains two main ingredients: on the one hand a solving algorithm and on the second hand a theoretical model.

Different methods (algorithms) exist to simulate the flow of granular media.

Some are based on the continuum representation of granular assemblies and solve the continuous mechanics equations, most of the time by means of the finite element method (such as the codes Gerris or Basilik developed at ∂'Alembert Institute at Pierre-et-Marie-Curie University, see [29] and references therein). These attempts are based on the µ(I)-rheology -which is a local model-and therefore cannot properly render the behavior of granular medium in the regimes where non-local effects are dominant.

NUMERICAL SIMULATIONS OF GRANULAR MEDIA

Another common way to simulate granular flows is using Discrete Element Method ("DEM", also sometimes called "Molecular Dynamics Simulations"): for each grain the fundamental principle of the dynamics is solved. This methods enables to compute different types of interaction between the particles. One of the drawbacks of such method is the computation time which increase rapidly with the number of particle in the numerical sample.In order to study non-local effects in dense granular flows, we hence limit the size of our numerical system to ≈ 1300 particles.

Numerical method for Discrete Elements simulations

The core computation Molecular Dynamics, or Discrete Elements Method is a quite simple and intuitive algorithm. It is based on solving the dynamics of each particle individually. To do so, three steps are needed:

1. for each particle, search of the neighbor particles 2. computations of the applied forces on the particle and force total resultant ;

integration of the dynamic equations

In addition to the core computation, the solving algorithm also contains the boundary conditions.For instance in a plane Couette cell, if the simulation is pressure imposed, there is a feed-back loop in the code that adjusts the height of the cell in order to maintain the pressure at the desired value.

The particle and boundary dynamics are integrated using a Verlet algorithm [30] in our code.

In molecular dynamics, implementing the model consists in defining the interaction between the elements. The three most common interactions are:

• attractive-repulsive, which most of the time is characterized by a Lennard-Jones potential 1 ;

1 Lennard-Jones potential is a model of potential inter-atomic energy. It has the form:

V LJ = ε r m r 12 -2 r m r
• Hertzian repulsive interaction (a more elaborated type of contact that takes into account the non-linearity in the contact between particles).

• harmonic repulsive interaction (the one we used in the present study) ;

Contact forces: harmonic potential For the sake of simplicity, and because this conceptual representation is sufficient to reproduce the phenomenology we want to observe [21], we chose to implement in our simulations a harmonic potential. It corresponds to the most simple contact law [31] in which two deformable particles inter-penetrate (cf. figure 3.2). This interaction is linear and based on a damped spring system.

Figure 3.2: Model of so-called "soft spheres" (see [31]): Harmonic interaction between two grains. Dry contact model: k n is the normal spring constant ; k t is the tangential spring constant ; g n is the normal viscous damping constant and µ p is the contact friction coefficient. δ stands for the inter-penetration.

the r -12 term is repulsive (Pauli repulsion) and is dominant at short range whereas the r -6 term is attractive at long range (Van der Waals or dispersion force).

NUMERICAL SIMULATIONS OF GRANULAR MEDIA

Its potential therefore takes the form:

V ij = k 2 (r i + r j -d ij ) 2 H (r i + r j -d ij ) , (3.2)
where H is the Heaviside function, so that forces are zero between the two particles if they are not in contact with each other. r i and r j are the radius of the particles and d ij the distance between their center. In this expression, r i + r j -d ij = δ is the inter-penetration when the two particles are in contact and k is the spring constant of the latter (for normal spring, k is k n ).

According to the soft sphere model represented figure 3.2, during a contact between a particle i and a particle j, the force resultant of the action of j over i writes as the sum of a normal and a tangential contribution:

F j→i = F n j→i + F t j→i . (3.3)
The normal component writes:

F n j→i = k n δ + g n δ n ij . (3.4)
In the latter, the first term of the right hand side characterizes a linear elastic contribution while the second term describes viscous dissipation. During a collision between two grains of mass m i and m j , at relative impact velocity v 0 , the inter-penetration is hence solution of:

m ef f δ + g n δ + k n δ = 0 , (3.5)
where m ef f = m i m j m i +m j is the effective mass of the system. Such an ordinary differential equation is stable for

g n ≤ 2 m ef f k n , (3.6)
and its solution writes:

δ(t) = v 0 τ c π exp - g n 2m ef f t sin π t τ c , (3.7)
where τ c is the time lapse of the contact: .8) This characteristic time of contact is crucial in the definition of the time step dt for the computation of the equations2 .

τ c = π kn m ef f -gn 2m ef f 2 . ( 3
During the contact, one can characterize the energy dissipation by defining a restitution coefficient as the ratio of the normal relative velocities after and before the collision:

= δ(τ c ) v 0 = e -gn 2m ef f . (3.9)
From 3.8 and 3.9 we get gn = -2 ln ( )

m ef f kn π 2 +(ln ) 2 .
Restitution coefficient has a value between 0 and 1, so the stability condition 3.6 is verified.

The tangential component is also given by a harmonic potential, of spring constant k t . The tangential component is therefore proportional to the tangential relative displacement δ t between the particles. In order to model friction between particles, we implement the Coulomb criterion which is based a critical force threshold for the tangential force component:

F t j→i = min (k t δ t , µ p F n j→i ) t ij , (3.10)
where µ p is the microscopic friction coefficient (or "inter-particle friction"). The maximum value of the tangential component is therefore limited to µ p F n j→i . In all our simulations, we take µ p = 0.5 when dealing with frictional particle (µ p = 0

Dimensional analysis

In order to analyze the physical quantities that characterize the granular system and its rheology, it is important to consider the microscopic and the macroscopic quantities describing the system.

The microscopic parameters are the normal and tangential spring stiffnesses k n and k t , the viscous damping g n the restitution coefficient (for hard visco-elastic spheres, = 0.9 but previous studies do not show a significant influence on the dynamics for lower values, so we always take an intermediate value = 0.4 for our computations), and the microscopic friction coefficient µ p . The length unit is the average grains diameter d (taken equal to 1) and the reference mass is the average mass of the grain, taken at m = 1.

The macroscopic system is characterized by a shear strain γ, a confining pressure P (which here is equal to the normal stress) and the compaction φ. Considering dimensionless quantities allows to better understand how the physical parameters come into play regarding the behavior of the system. Therefore they are used as control parameters. For instance, the natural scale for pressure is k n , therefore we can characterize the system by a dimensionless number: its softness S = k n /P . It reflects the role of the elasticity of the grains. In the following, we only focus on the limit of rigid non-deformable particles, which practically implies k n /P > 3.10 3 (see figure 3.3 which shows the evolution of the critical friction coefficient µ c with respect to the softness parameter S). Regarding the microscopic quantities, we take the microscopic friction µ p equal to 0 (frictionless particles) or 0.5. The order microscopic quantities (k n , k t , g n ) are set according to remain in the limit of rigid particles, for a constant restitution coefficient = 0. cially when we consider a gravity field which introduce a pressure gradient within the medium, hence a softness gradient. For our simulations of the incline plane, we set these settings so that the maximum softness, at the very bottom of the medium, is still in the regime of rigid particles.

From the point of view of the rheology, one seeks to find a relation between the shear rate (we imposed for instance) and how the fluid resists to this deformation, i.e. its shear stress. In the system, the only scale of energy is the pressure, therefore at the stationary regime, the dimensionless shear stress is τ /P . And it must be a function of the dimensionless shear rate I, hence the relation:

τ P = µ(I) (3.11)
The inertial number I = γ m P d compares the macroscopic shear rate to the ballistic time of a grains moving of its diameter under confining pressure P . As in our simulations m = 1 and d = 1, the inertial number writes I = γ/ √ P .

For similar reasons, the volume fraction φ -which is also a dimensionless quantityis also a function of the inertial number I (see first chapter and references therein).

NUMERICAL SIMULATIONS OF GRANULAR MEDIA

Preparation

The preparation process lies in generating particles which diameter that are randomly distributed according to a uniform law. In our preparation, the polydispersity is of 20%. Particles are disposed along an hexagonal pattern in the domain, with a random noise to their position in order to initiate the disorder of the preparation. Then we let grains fall under gravity against the bottom wall in order to fill the numeric cell.

After the step of the generations of the particles, we let the system relax and reach a stationary state (see next section). In the present work, we study the system once it is at stationary state, which means that all the observable quantities are constant with respect to time. Figure 3.4 presents the evolution of the kinetic energy (

Stationary state

E k = 1 2 N i=1 m i v 2 i )
and the total number of contacts in the system with respect to time. Each of these quantities reach a stationary value after a certain time. From this moment, we assume that all the average of the physical quantities characterizing the system (such as the velocity field, the stresses and the strain) are independent of time. We therefore evaluate them by computing there value averaged over time and space.

Test case: non-local rheology in a plane shear cell

In the following, we present the process by which we calibrate the parameters of the non-local rheology. This set-up has first been proposed by [15].

As non-local effects are exhibited in the vicinity of the transition, we first focus on a custom set-up of a 2D planar shear cell, which boundaries at the top and bottom walls are modified in order to be able to tune the yield parameter Y b in the bulk. The two other boundaries (left and right side of the cell) are periodic.

In this section, we first present the numerical set-up of this customized plane shear cell. We then compute the non-local rheology and present the process by which we calibrate the non-local model.

Plane shear cell: setup

The plane shear cell setup we use is a 2D system with two parallel rough walls moving in opposite directions, along axis x (see figure 3.5). Roughness of the walls is made simply by sticking the outer 3 layers of grains one to the others on each sides. All grains are of mass = 1. The distribution of diameters follows a square distribution d = 1 ± 20%. Boundary conditions on the sides are periodic (along direction x).

The system we use has a width of ≈ 35d and a height of ≈ 50d, for a total particles of 1300 particles (wall included). In the present study, we present our results for frictionless particles.

Pressure in the cell is imposed by adjusting the position of the "top" wall (the one that is also moving transversally to impose the shear rate), therefore the height of the cell is fluctuating, but these fluctuations have an amplitude smaller than an average grain diameter (see figure 3.6).

For a given stress configuration Y b , we compute the simulation until a stationary regime is reached (see section 3.1.4). We present figure 3.7 the velocity, volume fraction, stresses and µ profiles average on the duration of a simulation of a steady Figure 3.5: Setup for the plane shear cell. Pressure is imposed, so as the shear strain (wall velocity). This is a snapshot of a simulation of hard frictionless particles. Colors correspond to the instantaneous velocity along axis x. In order to exacerbate these non-local effects, [15] proposed a custom set-up for the plane shear cell that we present in the following and use as a check of our DEM code and numerical protocol. We will then present the strategy we adopt do calibrate the whole rheology, hence the non-local contribution, by fitting all the runs we computed simultaneously (runs above and below the threshold Y = 1) and therefore extract the parameters of the rheology.

Custom shear cell

Velocity profiles generated from simulations of plane shear can exhibit non-local effects, but it is not the most obvious way to observe them as they are more marginal in the flowing case (Y > 1). The idea here is therefore to study a system in which several regions with different stress state communicate with each other and mechanically interact.

By means of gravity-like forces applied to the grains located in two buffer zones located close to the walls (see figure 3.8), we manage to impose the profile of the yield parameter Y. These forces are oriented downward at the top of the cell, and upward at the bottom. In the bulk of the cell, the pressure is constant. The shear stress is not imposed in this set-up, but it is directly related to the shear rate we impose by setting the transversal (top) wall velocity: u wall .

We can therefore tune the profile of the yield parameter Y in the bulk and set it below the threshold. In the buffer zones Y > 1. Therefore, even though in the As explained previously, we tune the buffer zones so that we can observe velocity profiles even below the flow threshold and therefore exhibit non-local effects. In the following, we present the strategy we adopted to fit the non-local model.

Calibration of the rheology

We previously explained that in such configuration, the velocity profile writes 3.12:

u P/ρ g = I ∞ z + δI sinh(z/L) , (3.14
• Y > 1      I ∞ = 0 = L (Y-1) 2Y (3.15) • Y < 1      I ∞ = 0 = L χ -1 (Y)) (3.16)
using the phenomenological expression of χ(κ) previously calibrated in [24].

Recall that the yield parameter is define by Y = τ P µc . It thus depends on the value of the critical friction coefficient µ c , which is itself a parameter of the rheology we are fitting. The adjustment of the rheology must therefore be iterative, as one has to know whether the profile characterizes a configuration above or below the threshold Y = 1. The fit function is hence implemented using the conditional statements equations 3.15 and 3.16 and the expression of the velocity 3.12 and all the runs -at different values of Y-are fit simultaneously by constraining the global parameters of the model (i.e. µ c , a and ) to be the same for all profile as they all exhibit the rheological behavior of the same medium.

It the following, this fit technique by which we constrain the model to share some parameter values on several datasets will be mentioned by the expression "global fit" or "cofit".

From this fit process, we therefore extract the parameters of the rheology. We find µ c = 0.0953 ; a = 5.51 and = 4. 22. These values are very close from those found by [15]: they found µ c = 0.094, a = 5.51 and = 2. 8. The notable difference is on the calibration of the non-local term . This discrepancy may be due to a different implementation of the microscopic interaction in each code. Also, in order to better calibrate such parameter, one would run further more simulations on systems exhibiting more non-local effects and a more dominant role of non-locality, We can also extract the interim fit parameters that are I ∞ and the relaxation length L, in order to recover the µ(I) curve (as on figure 3.11) and exhibit the divergence of the relaxation length in the vicinity of the transition Y = 1 (see figure 3.12).

Note that the difference of the value of the parameter a obtained from the global fit and its value obtained from the fit of the µ(I) phenomenological law is no more than 5% of the parameter value. This shift lies in the fact that in the global fit, this parameter is fit by the simultaneous fit of all the datasets (above and below Y = 1), whereas the µ(I) curve only characterizes the runs for which I ∞ = 0, i.e.

the configurations in which the yield parameter is greater than 1.

The consistency of the model on both side of the transition is more visible when we plot the measured relaxation lengths L in respect of the yield parameter in the bulk (see figure 3.12). We would like to draw the attention of the reader to the importance of finite size effects : the height of the cell we used in our simulations is of ≈ 50d. It is therefore meaningless to consider datasets for which the apparent relaxation length is greater than the cell height. We therefore only consider in the global fit process only velocity profiles for which L/d < 50.

In order to calibrate the medium more accurately, one could run more simulations on configurations below the yield, so that we could reconsider the determination of the non-local functional χ which is a function of the relative fluidity κ.

Conclusion on the test case

We have just presented the test case of the plane shear cell set-up, customized in order to exhibit non-local effects. We showed how the velocity profiles in such system display relaxation lengths and how to relate them to the non-local term of the model.

We used a novel implementation of the fit algorithm in order to optimize the fitting of the theory on all our run simultaneously. We will now focus on the study of the incline plane setup and use a similar approach to characterize the flow profiles. From the description of the non-local model hence fit, we will present a measure of the boundary condition at the free surface of the flow.

Inclined plane: determining the boundary conditions at the free surface

In this section, we study the incline plane set-up in the framework of the non-local rheology, and we present a way to measure the boundary condition at the free surface. Gravity is inclined by an angle θ. Control parameters of the numerical experiment are θ and the microscopic friction between particles (µ p = 0 or 0.5). On left and right sides of the cell, boundary conditions are periodic, and the bottom of the cell is made of fixed grains (in black) in order to create a rough/no slip boundary condition. On this graph, colors characterize the velocity of each particle along the main flow direction x.

Setup and parametrization

We consider an inclined plane with periodic boundary conditions on the left and right and a rough wall at the bottom (see figure 3.13). Dimensions x and y are respectively along the longitudinal dimension of the inclined plane and along its height.

Gravity g is downward and tilted with an angle θ with respect to the vertical.

In this section, we will first describe how we precisely locate the bottom of the cell, from a physical point of view and how we compute the height of grains. The system is 2D and described in the ( x, y ) plane in the numerical code (hence in the raw data). As the theoretical model we develop focuses on the free surface, we

INCLINED PLANE: BC AT THE FREE SURFACE

reverse and set the vertical axis, switching from axis y which is upward and which zero does not correspond to any physical constrain, to axis z which is downward and set so that z = 0 at the free surface.

In the following, we explain these steps on a typical frictionless run (1100 free grain numbers in the cell + 64 in the wall ; cell width W = 35d ; k n = 10 5 ; k t = 5.10 4

; g n = 191.833 ; µ p = 0 ; gravity is inclined from vertical by an angle θ = 0.2).

Stress state

The momentum equation writes:

∇ • σ + ρ g = 0 (3.17)
Therefore we obtain the two equations:

     ∂ z σ xz + ρg sin θ = 0 (x) ∂ z σ zz + ρg cos θ = 0 (z) (3.18)
We choose to localize the free surface at z = 0. Therefore the integration of the previous equation assuming: these basic checks, we have verified that the stress configuration is well known.

     P (z = 0) = 0 τ (z = 0) = 0 (3.19) implies:      σ zz = P (z) = ρg cos θ z σ xz = τ (z) = ρg sin θ z (3.
The effective friction coefficient µ should therefore follow exactly the theoretical prediction. However, on figure 3.14 we see a slight deviation of the numerical result.

INCLINED PLANE: BC AT THE FREE SURFACE

Similarly, a steep variation of the volume fraction φ in the vicinity of the free surface (figure3.15) is observe. This motivates further investigations in the definition of the free surface boundary condition of the system.

In the following, we explain how we model the flow using the non-local rheology and we present a method to measure the boundary condition at the free surface. It is therefore important to notice that, from system 3.20 we get:

τ P = tan θ. (3.21)
As the yield parameter is defined by Y = τ /(µ c P ), we can therefore conclude that Y is only a function of the tilt angle θ. The latter is therefore the control parameter of this set-up.

However, it is important to remember that granular media exhibit a metastable behavior, especially near the jamming transition.This metastable behavior is very visible on the incline plane, at small heights of grains. When reaching a characteristic height h stop which depends on the tilt angle, the flow tends to arrest ( [33,34,35]).

For a given θ, decreasing h surprisingly tends to bring the system closer to arrest.

Numerical protocol

We therefore compute dense granular flows at diverse tilt angles θ, but also, for each angle, we varied the height of grains, so that for each value of the yield parameter Y, we have several runs. If we manage to adjust the flow curves for relatively high h, with the full non-local equation, then tests of lower heights, closer to jamming and in which the boundary condition at the free surface is predominant will be provide a measure of the actual boundary condition in z = 0.

For each inclination angle θ, we start from the same preparation, at rest, with ∼ 1200 particles (on a cell width W/d ≈ 35). We incline the gravity and let the flow become stationary (we call this first run, that is needed to reach the stationary flow the "transitional run"). Once we made sure it is the case (see 3.1.4), we launch CHAPTER 3. NON-LOCAL RHEOLOGY: NUMERICAL APPROACH 117 two runs. On the one hand, a long run, in the same conditions but starting from the last configuration computed by the previous run (and which we are therefore sure is in a stationary state). This run provides statistic in our flow measurements (stress state, means velocity profiles...). In the following we will refer to such run by the expression stationary run. On the other hand, we numerically extract some grains3 from the free surface in order to decrease h and then launch a new transitional run.

Doing so, little by little we decrease the height h of flowing grains, until the first premise of the metastable state (arrest of the flow during the simulation) appear. In the following, we present results only in the frictionless case. Frictionless case.

We therefore obtain velocity profiles such as the ones presented figure 3. 16.

Following the definition of the inertial number, we compute the corresponding I(z) profiles (see figure 3.17).

Figure 3.17: Incline plane. Inertial number profiles for θ = 0.2 -0.18 -0.16 -0.14 and for several flowing heights h (all runs from h ≈ 30.3 and then -50 grains between each consecutive run). Frictionless case.

Local rheology

We recall the expression of the local rheology:

τ P = µ(I) = µ c (1 + a I n ) . (3.22)
Which leads to:

I = tan θ µ c -1 1 a 1/n . (3.23)
We recover the fact that, for the inclined plane, the local rheology predicts a constant inertial number. In the following, we will refer to this quantity by I b . From the definition of the inertial number the shear rate therefore writes:

γ = - √ P d   tan θ µc -1 a   1/n . (3.24)
(we know from the experiment γ < 0)

By integration we then get the theoretical velocity profile according to this rhe-CHAPTER 3. NON-LOCAL RHEOLOGY: NUMERICAL APPROACH 119 ology:

u(y ) √ gd = 2 3d   tan θ µc -1 a   1/n 4 π φ cos θ   h d 3/2 - y d 3/2   (3.25)
As all physical quantities in this expression (3.25) are known (providing we have a calibration of the rheology, hence a value of µ c ), it therefore provides a completely determined expression of the local velocity profile. In dashed line is the prediction from the local rheology, using the calibration we previously obtained (µ c = 0.0953 and a = 5.51). We can clearly see the discrepancy of the numerics from the prediction. This is particularly remarkable in the region of the free surface. Inset figure 3.18 also shows a variance between the local prediction and the numerical result. These inconsistency of the local rheology to describe the flow on the incline plane justifies the following developments, where we probe the non-local model in order to accurately describe the flow and define the physical conditions at the boundaries of the system.

INCLINED PLANE: BC AT THE FREE SURFACE

The local rheology predicts an homogeneous value of the inertial number I = I b in the whole depth 4 . Also, we can see on figure 3.18 that the velocity profile it predicts can not recover the whole flow, especially in the neighborhood of the free surface.

In the following, we apply the non-local rheology as an attempt to describe the flow along the whole height of flowing particles.

Non-local rheology

We recall that the non-local rheology writes at linear order:

Y = µ(I) µ c (1 -κ) (3.26)
where κ = 2 ∇ 2 I I . In the current configuration the non-local equation is therefore:

1 - tan θ µ c (1 + aI n ) I 2 = d 2 I dz 2 . (3.27)
Above the critical conditions, the linearization of equation 3.27 around a base state

I = I b + δI writes: L 2 d 2 δI dz 2 -δI = 0, (3.28)
where as said previously, L is the relaxation length:

L = Y b n (Y b -1) , (3.29)
with Y b = τ P µc in the bulk. Equation 3.28 has solutions of the form (z = 0 been the free surface): Inset is the same plot but in semilog scale. Top graph is the residuals of the cofits.

δI(z) = C 1 e z/L + C 2 e -z/L . (3.30)
Therefore we fit the data by the following formula: We therefore cofit together all the inertial number profiles that are at the same tilt angle θ. For each value of θ, we therefore extract the value of the inertial number in the bulk, I b , and the relaxation length, L. On figure 3.19 we present global fit results on inertial profiles at same yield parameter (θ = 0.2) and for several heights of grains.

I(z) = I b + C 1 e z/L + C 2 e -z/L . ( 3
On figure 3.20 we present inertial profiles obtained from steady flows on the incline plane with the same height of grains but for several values of θ. The solutions of the fit using formula 3. 31.

We clearly recover the fact that all profiles of the inertial number I for a given θ However, plot of the residuals of the fit (as shown figure 3.19) show important variance between the numerical data and the fit in the vicinity of each boundary of CHAPTER 3. NON-LOCAL RHEOLOGY: NUMERICAL APPROACH 123 the domain: at the bottom and at the free surface. This suggests that the analytical expression by which me fit the profile if not valid in these regions.

We also report the cofitted values of the relaxation length. Again, we recall that The reader should note that equation 3.28 is the linearization of the non-linear differential equation 3. 27. As we cannot solve it analytically, it must be solved numerically. We use the built-in Matlab Runge-Kutta method ode45 to do so. As seen section 2.6, we need to initiate the numerical resolution. It is of course prescribed

L

INCLINED PLANE: BC AT THE FREE SURFACE

to start the code from the free surface as this is precisely what we are trying to characterize here. The reader should also note that no assumption on a boundary condition at the bottom was used in the theory. We hence emphasize the fact that the theoretical problem here does not have an explicitly boundary condition at the bottom: the only boundary condition we used is the one at the free surface, imposing

P (z = 0) = 0.
One should therefore initiate the resolution from elsewhere.

The domain over which we fit the analytical expression is located around I b . The common fit parameters are I b and L the relaxation length. Indeed, from the non-local rheology, we recall that L is a function of (from the calibration of the rheology)

and Y (which contains the rheology and the stress state as Y = τ µc P ). We note that it has an inflexion point located at depth:

z 0 = L 2 log -C2 C1 (3.32)
From the fit values of C1 and C2 for each run, and the measure of L and I b , we try to initiate the Runge Kutta numerical resolution of the complete non-local equation by using the linear approximation around the base state, located in the vicinity of the inflexion point. Up to now, we did not manage to recover the numerical results from the DEM simulation.

However, using the fit formula provided by the linearized equation, we manage to estimate the value of the inertial number at the free surface. As previously explained, this result should not depend on the height of grains in the cell. On figure 3.23 we present estimations of the value of the inertial number at the free surface with respect to the distance from the yield.

On this figure, we check that the boundary condition does not depend on the height of flowing grains. It seems that the value of the inertial number is selected with respect to Y -1. We can assume that the grains at the surface exhibit a "tac-tac" behavior, which would also explain the non-zero value of a in the linear fit.

In order to measure more accurately the boundary condition at the free surface, we now plan to used a shooting method.

Conclusion on the numerical approach and perspectives

In this chapter, we focused on two dimensional dense granular flows. We first explained how the numerical code we used is built. We then presented two numerical systems which allowed us to exhibit non-locality in different ranges of the yield parameter Y, and calibrate the rheology. On the one hand, the plane shear cell permits to calibrate the rheology in the vicinity and on both sides of the jamming transition, i.e. for Y < 1 and Y > 1. On the other hand, the incline plane also exhibit non-local effects, but a larger values of Y. We presented the process by which we quantify the model coefficients. Up to now, we show that the linear approximation of the non-local equation is non sufficient to properly measure the relaxation length L of the rheological profile. Indeed, even if this approximation seems to phenomenologically recover the flow profile in the incline plane, we still did not manage to find an Many (numerical) studies ( [36,37]) have been conducted to understand the jamming of a granular packing and the relation to elastic response. In [36], A. J. Liu, S. R. Nagel, C. S. O'Hern and collaborators have modeled granular packings by using diverse types of repulsive potentials to model their interactions at short range. In this approach, inter-particle friction is ignored. They use either an harmonic potential (α = 2) or an Hertzian potential (α = 5/2):

V ij = k 2 (r i + r j -d ij ) α H (r i + r j -d ij ) , (1.1)
where H is the Heaviside function, so that forces are zero between the two particles if they are not in contact with each other. r i and r j are the radius of the particles and d ij the distance between their center. In this expression, r i + r j -d ij = δ is the inter-penetration when the two particles are in contact and k is the spring constant of the latter.

At zero temperature and with no shear, they show the existence a critical volume fraction φ c (sometimes written φ j with "j" standing for "jamming") at which the transition between the two phases occurs. This critical volume fraction is directly related to the critical number of contacts Z iso defining the rigidity transition and above the transition there is a general relation, whatever the interaction potential and the spatial dimension):

Z -Z iso ∼ (φ -φ c ) 1/2 , (1.2)
where Z iso = 2D (D being the spatial dimension). For volume fractions below this value (φ < φ c ), particles do not overlap, and therefore contact pressure is zero. Elastic moduli of the packing are thus not defined, as it would for solids. For φ > φ c , particles overlap so contact pressure and hence elastic moduli of the packing are finite. They both increase with the distance from the jamming ∆φ = φφ c .

In the solid phase, the elastic moduli of the packing present power law dependences with the distance from jamming. Note that this is reminiscent to critical phenomena as we mentioned before. 2)). The closed symbols denote moduli calculated by forcing the particles to move affinely and the open symbols correspond to the moduli calculated after the system has relaxed. Slopes are as indicated (adapted from [2] with permission-copyright by the American Physical Society).

It is worth noting that many soft matter systems (pastes, emulsions) have shear moduli which are much smaller than compressional moduli-from an application point of view, this is a crucial property. Putting all this together, we conclude that the affine assumption gives the correct prediction for the bulk modulus (since k ∼ δ α-2 ∼ φ α-2 ), but fails for the shear modulus. This failure is due to the strongly non-affine nature of shear deformations: deviations from affine deformations set the elastic constants [2,20,30,43,62]. As we will see below, the correspondence between the bulk modulus and the affine prediction is fortuitous, since the response becomes singularly non-affine close to point J for both compressive and shear deformations (section 3.5.5).

Non-affine character of deformations.

Approaching the jamming transition, the spatial structure of the mechanical response becomes less and less similar to continuum elasticity, but instead increasingly reflects the details of the underlying disordered packing and becomes increasingly nonaffine [30]-see figure 4(a). Here we will discuss this in the light of equation (8), which expresses the changes in energy as a function of the local deformations u ∥ and u

⊥ : E = 1 2 i, j k i j (u 2 ∥,i j - δij α-1 u 2 ⊥,i j ).
To capture the degree of non-affinity of the response, Ellenbroek and co-workers have introduced the displacement angle α i j . 4 Here α i j denotes the angle between u i j and r i j , or

tan α i j = u ⊥,i j u ∥,i j . (14)
The probability distribution P(α) can probe the degree of non-affinity by comparison with the expected P(α) for affine deformations. Affine compression corresponds to a uniform shrinking of the bond vectors, i.e. u ⊥,i j = 0 while u ∥,i j = -εr i j < 0: the corresponding P(α) exhibits a delta peak at α = π. The effect of an affine shear on a bond vector depends 4 Not to be confused by the power law index of the interaction potential.

on its orientation, and for isotropic random packings P(α) is flat.

Numerical determination of P(α) shows that systems far away from the jamming point exhibit a P(α) similar to the affine prediction but that, as point J is approached, P(α) becomes increasingly peaked around α = π/2 (figures 13(b) and (c)). This is reminiscent of the P(α) of floppy deformations, where the bond length does not change and P(α) exhibits a δ peak at π/2. Hence deformations near jamming become strongly non-affine, and, at least locally, resemble those of floppy modes. Non-affinity of floppy modes and elastic response. Wyart and co-workers have given variational arguments for deriving bounds on the energies and local deformations of soft (low energy) modes starting from purely floppy (zero energy) modes [54,63]. They construct trial soft modes that are basically floppy modes, obtained by cutting bonds around a patch of size ℓ * and then modulating these trial modes with a sine function of wavelength ℓ * to make the displacements vanish at the locations of the cut bonds [30,54]. In particular, for the local deformations, they find [63]

u ∥ u ⊥ ∼ 1 ℓ * → u ∥ u ⊥ ∼ z, ( 15
)
where symbols without indices i j refer to typical or average values of the respective quantities.

The question is whether the linear response follows this prediction for the soft modes. The width w of the peak in P(α) is, close to the jamming transition, roughly u ∥ /u ⊥ because |α i jπ/2| ≈ u ∥,i j /u ⊥,i j if u ∥,i j ≪ u ⊥,i j . It turns out that the scaling behavior ( 15) is consistent with the width w of the peak of P(α) for shear deformations, but not for compression. There the peak of P(α) does not grow as much and a substantial shoulder for large α remains even close to jamming: the tendency for particles to move towards each other remains much more prominent under compression. Scaling of u ∥ and u ⊥ . The scaling of the distributions of u ∥ and u ⊥ has also been probed. The key observation is that in Figure 1.1: Bulk (K) and shear (G) modulus as a function of distance to jamming for two-dimensional bidisperse systems, according to the potential. The closed symbols denote moduli calculated by forcing the particles to move affinely and the open symbols correspond to the moduli calculated after the system has relaxed. Data from [36] and adapted by [37] According to the implemented interactions (Harmonic α = 2 ; Hertzian α = 5/2), authors report dependencies of the pressure and the elastic moduli (see figure 1.1) on ∆φ following the power laws:

             P ∼ (∆φ) α-1 K ∼ (∆φ) α-2 G ∼ (∆φ) α-3/2 (1.3)
From these considerations, we see that, whatever the potential of interaction, the exponent of the shear modulus is always larger than the exponent of the bulk modulus, so that:

G K ∝ (∆φ) 1/2 ∝ Z -Z iso . (1.4)
The shear modulus G gets progressively smaller than the bulk modulus when going

ELASTICITY OF JAMMED GRANULAR MATTER

to the jamming. This relation is called the shear anomaly and is a signature of the approach of the rigidity transition.

In this framework, for frictionless grains, pressure hence appears as good control parameter to reach the jamming transition.

However, for real granular material, the inter-particle interaction can be more complicated due to the presence of contact roughness. This can have several important consequences, like the emergence of solid friction and also changes in interaction laws. When the overlap between particles is at the scale of the contacts asperities [38,39], the interaction laws may not remain purely Hertzian.

Rigidity transition was generalized in the case of frictional contacts usually studied with either harmonic or Hertzian interactions. For infinite friction, the critical number of contacts is Z iso = D + 1. However for finite friction, it is more complicated as one has to take into account the number of contacts in the Coulomb cone.

Authors like Ellenbroek [40] have defined a critical Z µ iso , dependent of the microscopic friction coefficient µ, such that D + 1 < Z µ iso < 2D. An essential claim of this team is that the critical character of the rigidity transition remains when Z -Z iso vanishes. In particular the relation between the number of contacts and pressure can be generalized Z -Z iso ∼ P 1/3 (for Hertzian contacts). At vanishing pressures, this would imply an elastic anomaly as for frictionless particles. However, other authors, like Magnanimo [41], claim that at vanishing pressure, the number of contacts would saturate at a constant value eventually far from the generalized isostatic value. This means that even though the shear modulus could be significantly smaller than the bulk modulus, the scaling relation 1.4 would not be observed. Therefore a mean field behavior for the elastic moduli as described in the following should be observed at evanescent pressures.

Mean field approach

The Mean field approach, or effective medium theory (EMT) [42,43] is based on the affine approximation that the motion of each grain follows linearly the macroscopic applied strain. Assuming a Hertzian interactions between spherical grains, at a confining pressure P , the overlap δ scales as P 2/3 . This means that the effective inter-particle stiffness scales as P 1/3 . As a consequence, the elastic moduli scales also as P 1/3 . The full mean field calculation for a granular packing of compact fraction φ and number of contacts Z gives a dependence of the elastic moduli as:

K ∼ G ∼ E 0 (φZ) 2/3 (P/E 0 ) 1/3 , (1.5)
where E 0 is the material Young's modulus. The mean field theory thus does not predict the anomaly on the elastic moduli previously described.

More explicitly for spherical grains made of material which properties are µ g (shear modulus) and ν g (Poisson's coefficient), EMT predicts:

K = 1 3π √ 2 µ g 1 -ν g Zφ 2/3 P 1/3 , (1.6)
and

G = 3 5 1 + 3(1 -ν g ) 2 -ν g K, (1.7)
where would be 0 in the frictionless case and is 1 in the infinite friction limit.

Note, EMT would also predicts sound velocities scaling as P 1/6 .

However, Makse et al [44] have shown the failure of the mean field theory. They compare the mean field prediction with DEM results of Hertzian granular packing.

He artificially introduces a parameter α in order to continuously vary the influence of the tangential force contribution in contacts. He shows that the EMT prediction qualitatively reproduces the evolution of the bulk elastic moduli K (EMT over estimates the simulation results by ≈ 10%) but fails to predict the evolution of the 1.1. ELASTICITY OF JAMMED GRANULAR MATTER shear modulus at vanishing tangential contribution (which corresponds to frictionless packing). However, in a second series of simulations, he constrains the rearrangements of particles to be affine with respect to the applied strain rate, in order to artificially impose the affine assumption of EMT at the scale of particles. He hence manages to recover the EMT prediction, in particular for the shear modulus.

This finding is qualitatively coherent with the frictionless scaling of C. O'Hern et al [36] we previously presented about the anomaly of the shear modulus in the vicinity of the jamming transition.

On an experimental point of view, the dependence of elastic moduli with pressure has been probed by the measure of sound celerities in packing of glass beads [42,45,46,47]. These studies show that the mean field prediction is not valid and exhibit an effective power law of c ∼ P 1/4 , instead of the c ∼ P 1/6 expected from the mean field of Hertz prediction [45] (see figure 1.2).

An empirical relation between the mean number of contacts Z and the confining pressure Z(P ) in order to take into account this discrepancy.

There is however another way of thinking that introduces contact relations different from the Hertz scaling. For example, de Gennes' theory [38] takes explicitly into account the softness of an outer shell which can model an oxide layer or asperities inducing an effective reduced softness at the surface of the grains.

Noticeably, all these experiments were made under very high confining pressures and in a limited pressure range. In the following, we present an experiment we designed in order to reach very low confining pressures in a granular assembly of glass beads. We measure sound velocities of plane compressive waves at pressures down to ∼ 10P a and propose a theoretical framework in which we revisit the nature of particles contacts by taking into account the asperities at the surface of the grains as a soft superficial layer.

Experimental set-up

We took the opportunity of three CNES parabolic flight campaigns (March/October 2015 and October 2016) to study the acoustic propagation of a wave packet in a granular packing and probe the elastic response of a granular packing at low confining pressure. Because in a lab experiment the gravity induces a pressure gradient within the grain assembly, we need to be in a weightless environment in order to be able to control the confining pressure at a lower value than the one induced by the gravity gradient on Earth. The French national space agency (CNES) offers parabolic flight campaigns dedicated to science, on board of the Airbus Zero-G (Airbus A300 until late 2015 and since then A310).

Overview of the experiment

The principle of the experiment lies in propagating a sound wave through a cell containing a granular medium confined at a controlled pressure. From a physical point of view, the control parameters are the confining pressure within the medium and the characteristics of the generated wave (amplitude of the signal and frequency)

; the measurements are the characteristics of the sound wave transmitted through the cell.

In absence of gravity the confinement pressure can in principle be fixed at very low values. However, one has to keep in mind that for real zero-G flights, there is always a remnant G-jitter acceleration spanning ±5.0 10 -3 g in amplitude. To compensate for these acceleration variations, due to the ability of the pilot to adjust the plane trajectory to the desired parabola, we set up a feedback loop suited to maintain a constant pressure within the sample.

The experiment is thus composed of two racks, one containing the core experiment (i.e. the cell with its sensors and actuators, and the confining system) and the other embedding the power supplies for the electronics, the acquisition boards and the computer for user interface and the feedback program that controls the pressure within the medium (see figure 2.1). In order to monitor the effective source shape, a 1D accelerometer is fixed underneath the moving plate. Similarly, a 3D accelerometer is stuck in the lid, on the side in contact with the beads. This allows to record the transmitted wave. From the target pressure P set , one computes the corresponding force to impose on the grains in the ideal case of perfect "zero-G" (g current = 0). The g-jitter (g current ) creates a resisting force F bot = g current M t on the force sensor, due to the momentum of the moving components (granular sample, lid of the box, polymer cubes and the aluminum plate that is the the support of the force gauge). The deviation F from the target force is then converted to a command tension U set which is compared to the tension characterizing the force on the medium. The linear stage only adjusts the pressure when the following conditions are gathered:

• the deviation between the tensions characterizing respectively the force to impose on the medium U set and the output of the force sensor U real is out of a defined range (we set the tolerance at 0.1V from the target which corresponds to ±5P a.

• the force sensor do not saturate

• we set the g-target to zero (when running experiment on the ground, to probe the feedback loop when need to change this target, otherwise g current remains out of bounds).

Also, note that the experiment was also tested under terrestrial gravity (g = 1m/s 2 ), therefore there is a flag (a boolean) that enables to turn off the retroaction of the feed-back loop is the case of laboratory tests.

Sensors

The main measurements lie on the comparison of the wave characteristics, before and after it has propagated through the medium. To do so, a one-axis accelerometer is mounted on the vibrating plate (i.e. the source) and a three-axis accelerometer is embedded on the lid of the cell.

Jitter accelerometer We use a dedicated mounting of 3 low frequency variable capacitance accelerometers (MEGGITT 7290E, bandwidth 0-15kHz ; one along each axis) in order to monitor the environment apparent gravity. The z-axis measurement (which is parallel to the axis along which we propagate the wave and which corresponds to the yaw axis of the aircraft) is directly input to the feed-back loop which acts on the fine translational stage.

Force gauges Up to now we have used four Futek LSB200 of 250gr ams capacity each. They were displayed in a square, hanging on the stage plate and also screwed to the plate which is used to load the polymer cubes and therefore the cell lid. This set-up enables to reach a maximum pressure of 200P a in so-called "Zero-G" and a minimum load of 15P a with a ±5P a precision. A drawback of using several sensors is the balance of the load on the plate (one sensor can saturate and therefore the computed sum of forces does not reflect the actual pressure the stage imposes). For last campaign we chose to replace these forces sensors by a single weighting scale sensor. This choice offers a wider range of forces with the same precision and only one sensor which simplifies greatly the mechanical balance of the plates between which it is mounted, which guaranties a more homogeneous pressure on the medium. We also replace the pre-amplifier system with a compatible conditioner. This upgrade is also improved the precision of our set-up (by a factor 2). precise and accurate) calibration is embedded so that when connected to the amplifier, it reads the accelerometer data and set it up properly.

High frequency accelerometers

Control of the confining pressure

The confining pressure is monitored by force gauges that are mounted between the stage and the plate that sits on the top of the box (see figure 2.5 page 145). In the ideal case of a perfect "Zero-G" (no jitter), the sum of the measured forces over the area of the top wall in contact with the medium would give the confining pressure.

In reality the G-jitter has a double influence: on the one hand it creates a pressure gradient within the medium while we assumed the pressure to be homogeneous within the medium, and on the other hand it perturbs the feedback loop because of the inertia of its components.

Generation of a plane wave

We probe the medium by propagating a single frequency plane wave 1 through it.

This plane wave is generated by vibrating a rigid plane plate in a direction normal to

1 A plane wave is a field A( x, t) which take the form:

A( x, t) = f ( n c x -t)
with an arbitrary function f and n unit vector. The points of equal field value of A( x, t) always form a plane in space. This plane shifts with time t along direction n and with celerity c. its surface thanks to the piezo-electric pistons previously mentioned. We chose to constrain the number of carrier wave oscillations within one wave packet. Therefore the width σ of the Gaussian command is fixed by:

σ = 10 4f √ 2 log 2 (2.1)
The command signal sent to the actuators is a command in displacement, of form:

D(t) = A cos (ωt + φ) exp - 1 2 t σ 2 (2.2)
It is important to underline the fact that this command is in displacement. We will come back to this later but notice that when we will compare the emitted and received waves, our measurements will be in terms of acceleration as the sensors are accelerometers. Not displacement.

Choice of the wave carrier frequency

We propagate a single frequency (carrier)

wave in the medium and probe several frequencies to check this hypothesis.

Therefore we need to choose the frequency of the carrier wave. It is limited by

EXPERIMENTAL SET-UP

the two length scales of the system: the wavelength must be larger than the grain diameter (so the mean field hypothesis is valid) and smaller than the height of the box (which correspond to the length along which the wave propagates). Therefore:

d ≤ λ ≤ H (2.3)
where d is the diameter of the beads (distributed in a diameter range from 1mm to 1.3mm) ; H 8cm is the height of the box ; λ the wavelength. If we propagate the trend extrapolated from Jia [45] data we can estimate the celerities for the pressures we wish to test: c min 10m/s at minimum and c max 10 3 m/s for highest pressures.

With H 8cm this leads to f min 1.250kHz and f max 10kHz.

Note that we also checked that this was compatible with the actuator frequency domain. In order to make sure that the source stays away from its resonance we chose frequencies closer to the f min . In the following, we only presents results obtained at a frequency (of the carrier wave) equal to 4kHz. The minimum sound celerity we can propagate is therefore c min = f d = 4000 × 0.001 = 4m/s and the maximum c max = f H = 4000 × 0.08 = 320m/s. is simply computed from the measure from the force sensor over the lid surface.

Conducting the experiment

Pressure at the bottom is estimated by P bot = P top + M grains g current S (where S is the area corresponding to the horizontal cross section of the cell 2 ). During the flight, g current is supposedly zero, therefore P top = P bottom . When testing the experiment on 2 it is assumed to be equal to the surface of the lid.

DATA PROCESSING AND SELECTION

the ground, we can simulate this condition by turning off the "keep g value" setting.

At the top of the window, the raw (tension) output of the force gauges are monitored 3 , as well as their mean value in real time. We also this value with time, and below in the window we monitor the raw output of the source and the reception accelerometers in order to visualize in real time if we propagate wave packets. Note that in the case of small amplitudes, the ratio signal over noise sometimes prevents us from distinguishing the transmitted wave packet, by the filtering and the posttreatment process we make enable to recover some usefull signal as we will see in the next section.

Data processing and selection

Calibration of the wave sensors

The main question here is to rigorously define what physical quantities are measured in the experiment, especially regarding the accelerometer embedded in the lid within an envelope of silicon. This is a key question as we want to compare with a theoretical model.

As shown on figure 2.9 the source signal comes from an accelerometer which is fixed under the vibrating plate, at the bottom of the cell. Therefore it gives a measurement of the acceleration of the latter, i.e. the second derivative of its displacement.

At the opposite, the reception accelerometer (see figure 2.10) which is embedded in the lid is isolated by a thin layer of silicon in its compartment. This was made to prevent an acoustic bridge with the whole lid: we wanted to make sure we measure the acceleration of the grains, not any other signal. The lid accelerometer is thus 3 as we changed the force measurement setup by replacing the four Futek 250gr gauges by a 5kg Scaime sensor, we chose du minimize the changes in the code, hence duplicating the value of the single force sensor on the 4 channels and only modifying the α F coefficient that converts its output tension to units of force. This also explains the "4" coefficient within the feedback loop diagram figure 2.3 page 141. It is encrusted within the lid, in a silicon bed, and it comes out a little from the bottom surface of the lid so we make sure it is in contact with the grains. z-axis of the accelerometer is oriented toward the top to the bottom of the cell. Note the polymer cubes that ensure the damping of between Linear stage and the cell lid, in order to prevent rigid contact and therefore an acoustic bridge with the medium.

DATA PROCESSING AND SELECTION

a damped mass-spring system, excited by the acoustic pressure of the grains. Note that it is not the acceleration of the grains. Indeed, with respect to the wave or the (low frequency) fluctuations of the effective g, the lid is supposed to be fixed (and no displacement implies no acceleration). The reception accelerometer is therefore a pressure sensor that we need to calibrate. We shock the accelerometer in its compartment in the lid. The system is free to oscillated and beside the first shock pulse, received no other excitation. Red is the time response and green is the theoretical prediction.

Probing the reception sensor A first test is to measure the free response of the system. We settle the lid horizontally, with the accelerometer hanging in place and we shock (hit) the sensor. Figure 2.11 shows the temporal response of the system.

For such system, Newton's second law writes:

d 2 Z dt 2 + 2ξω 0 dZ dt + ω 2 Z = F (t) m (2.4)
where F (t) = 0 in the present case as the is no driving force ; m being the mass of the studied system (the accelerometer) ; ω 0 the undamped angular frequency and ξ the damping ratio. The damped oscillations hence take the form:

Z(t) = Ae -ξω 0 t sin ω 0 1 -ξ 2 t + φ (2.5)
Amplitude A and phase φ being given by the initial conditions. Fitting with the data (see figure2 .11) gives ω 0 6740r ad.s -1 and ξ 0.1. We find ξ < 1/ √ 2 therefore the system is underdamped. It therefore has a resonant frequency which is given by ω r = 6663r ad.s -1 . ω 0 and ω r respectively correspond to f 0 = 1070Hz and a resonant frequency of f r 1060Hz. One also likes to mention the corresponding quality factor Q = 1 2ξ 4.8 which characterized the ratio of stored energy over the energy lost per unit cycle 4 .

Figure 2.12: Response curve in amplitude, for the accelerometer. Propagation for a sine wave through air (empty cell). The two colors correspond to two separated data sets (me moved the lid and the set-up in the meantime).

A second and more complete test is to measure a whole frequency response of he system, by exciting it with a sine signal spanning several frequencies, from the lowest to the highest allowed by the actuation devices (the PZT). To do so, me empty the cell and we propagate a wave from the source to the lid (within which the reception

DATA PROCESSING AND SELECTION

accelerometer is embedded). This allows to only consider the cell itself and not the granular medium. At constant source amplitude, we propagate a sine wave for which we change the frequency. Figure 2.12 presents the frequency response of the system, left completely free and without any contact with granular matter: the cell is empty (see inset figure 2.12) and the lid is facing the vibrating plate without being in direct contact with the cell (we mounted it over foam cubes, in yellow on the graphics).

We send sine waves of constant command amplitude and display here the restitution coefficient with respect to the sine frequency.

In this test, the accelerometer is not in contact with any grain. It is free to move in is slot. The contact with the granular medium might have an influence on the dynamical response of the moving mass of the sensor, and hence modify its vibrational response. However, in our tests in the cell filled with grains, we did not see any evidence of a drastic change in the resonant frequency of the set-up.

We can then assume reasonably that the momentum equation on this mass-spring system is dominated by the acceleration so that:

Γ = S lid P/m
where Γ is the acceleration of the moving mass m and P the acoustic pressure provided by the grains.

Propagation model in the acoustic cell (compression)

In this section, we choose, at first order, to neglect the gravity. Therefore the stress state is only characterized by σ = σ zz with z axis along the vertical direction and from the source (bottom plate) to the reception, at the lid of the box.

The momentum balance writes:

∂σ zz ∂z = ρ ∂ 2 u ∂t 2 (2.6) CHAPTER 2. PROPAGATION OF A PLANE WAVE 155
where u is the displacement of particles along z axis. The normal stress σ zz is linked to the longitudinal deformation ∂u ∂z by the constitutive relation:

σ zz = K ∂u ∂z (2.7)
Equation (2.6) thus writes:

∂ 2 u ∂t 2 -c 2 ∂ 2 u ∂z 2 = 0 (2.8)
where

c 2 = K ρ
is the celerity of sound in the effective continuous medium. K being the bulk modulus and ρ the material effective density.

As K and c are actually the quantities we want to measure in this study, it is preferable to express K with respect to ρ (that we can measure independantly) and c:

K = ρc 2
PDEs like (2.8) have solutions of form (plane progressive waves):

u(z, t) = f (t - z c ) + g(t + z c )
This is a second degree PDE therefore the problem has only two boundary conditions:

• u(0, t) = u 0 (t) • u(H, t) = 0
Where φ is the phase (of the carrier wave of pulse ω) at the center t m of the Gaussian wave packet (of displacement amplitude A and width σ).

The first BC (z = 0) of PDE (2.8) implies:

u(z, t) = u 0 (t - z c ) (2.9)
and therefore

u 0 (t) = -Ae -1 2 ( t-tm σ ) 2 ω sin (ωt + φ) + t -t m σ 2 cos (ωt + φ) (2.16
) which gives at the source (z = 0, at time t), the acceleration:

u 0 (t) = -A ω 2 cos (ωt + ψ) exp - 1 2 t -t m σ 2 (2.17) + 2A ω (t -t m ) σ 2 sin (ωt + ψ) exp - 1 2 t -t m σ 2 + A (t -t m ) 2 -σ 2 σ 4 cos (ωt + ψ) exp - 1 2 t -t m σ 2
Therefore the vertical stress at the wall is, in fine:

σ zz (H, t) = +2ρcAe -1 2 t-tm -H c σ 2 × ω sin ω(t - H c ) + φ + t -t m -H c σ 2 cos ω(t - H c ) + φ (2.18)
At the limit at which the lid (top wall) of the box is completely non-mobile, the momentum equation of the accelerometer settled in its compartment with a matrix of silicon gel writes:

m d 2 Z dt 2 = -kZ -α dZ dt -σ zz (H, t) S (2.19)
with S the surface of the accelerometer in contact with the grains (typically ∝ πD 2 /4

with D the dimension of the accelerometer) ; m, k and α respectively the moving mass, spring and damping coefficient of the system.

According to the fact that ω >> ω 0 (where ω 0 6700 s -1 is the natural of the system {accelerometer + silicon} in the lid, see page 152 paragraph 2.2.1 "probing the sensor"), we can make the assumption that the first two terms of the right hand side of 2.19 are negligible with respect to the acoustic pressure force P (H, t) S.

DATA PROCESSING AND SELECTION

Therefore equation 2.19 leads to:

d 2 Z dt 2 = - S m 2Aρc e -1 2 t-tm -H c σ 2 × ω sin ω(t - H c ) + φ + t -t m -H c σ 2 cos ω(t - H c ) + φ (2.20)
by which we should fit the wave packet signal coming from the lid accelerometer (reception), in the non-dispersive case.

In the dispersive case, the carrier wave do not necessarily travel at the same celerity as the modulation (i.e. here the Gaussian envelope). We thus write c φ the phase velocity, i.e. the celerity of the carrier wave and equation 2.20 becomes in the dispersive case:

d 2 Z dt 2 = - S m 2Aρc e -1 2 t-tm -H c σ 2 × ω sin ω(t - H c φ ) + φ + t -t m -H c σ 2 cos ω(t - H c φ ) + φ (2.21)

Raw data and data treatment

Synchronization of the data During the whole flight, we constantly record the time and any output from the force sensors and the accelerometers (along the 3 axis, at the position of the experiment which is slightly different from the acceleration recorded by the plane equipment and is also sampled at a higher frequency). We also record the displacement of each stage separately and of course the signals from the high-frequency accelerometers that measure the source and reception of the wave packet in the cell.

For each parabola, a folder is automatically created that contains signal.txt Selection of the time windows (refine crop of the parabola) For each parabola one needs to isolate the time window within which a criterion based on the "quality" of the "zero-G" is satisfied: remnant variations of the gravity must stay in a 0±0.05G

range.

Also, because it is hard to predict whether the calibers of the accelerometers will be properly set in advance, we allow the operator of the experiment to adjust manually their value during the parabola, in real-time. It is therefore important to be able to manually crop the previous time windows accordingly.

Filtering Beside the filtering role of the cubic polymers inserted between the lid of the cell and the moving plate hanging on the force gauges, we need to filter the signal in order to get rid of the remnant G-jitter and any other noise or parasite signal that would jeopardize the fitting of the expected signals. We therefore use a band-pass filter that only keeps frequencies within the range [f -1500Hz ; f + 1500Hz]. Then we use the Hilbert transform 5 to get the envelope of the signal. 5 The Hilbert transform is a linear operator that propagates a function (signal s(t)) to the complex plane so that it satisfies the Cauchy-Riemann differentiability necessary and sufficient condition in this space. It is often used to calculate instantaneous attributes of a time series, especially the amplitude and frequency. It is defined by the convolution:

H(s)(t) = 1 π +∞ -∞ s(x) t -x dx (2.22)
Here, we just use the Matlab function of the Hilbert transform, that provides the envelope for a given signal.

DATA PROCESSING AND SELECTION

Computing the physical parameters: fit of the data Recall from 2.23 that at the source (x = 0, at time t), the acceleration (and therefore the function by we fit the signal from the accelerometer) writes:

u 0 (t) = -A ω 2 cos (ωt + ψ) exp - 1 2 t -t m σ 2 + 2A ω (t -t m ) σ 2 sin (ωt + ψ) exp - 1 2 t -t m σ 2 + A (t -t m ) 2 -σ 2 σ 4 cos (ωt + ψ) exp - 1 2 t -t m σ 2
Similarly, we recall that at the reception (x = H, at time t, see equation 2.21), the signal has to be fit by:

d 2 Z dt 2 = - S m 2Aρc e -1 2 t-tm -H c σ 2 × ω sin ω(t - H c φ ) + φ + t -t m -H c σ 2 cos ω(t - H c φ ) + φ (2.23)
We bring to the reader's attention that in this formula, the phase of the carrier wave (the cosine/sine parts) does not necessarily travels at the same speed as the Gaussian envelope. When we derive the displacement field u(z, t) with respect to time t, we implicitly assumed it was the case, i. e. that the medium was not dispersive. In order to leave some flexibility on the fit, regarding this assumption, we unconstrain the model by allowing the carrier wave to travel at a different speed then the envelope. We therefore distinguish the group velocity c (at which the envelope travels) from the phase velocity c φ at which the carrier wave travels.

The fit procedure we adopted lies on these considerations. First we fit the source:

1. Gaussian fit of the envelope of the source (first term of formula 2.23): it provides a first estimate of the amplitude of the acceleration A s -Aω 2 (A is the amplitude of the displacement, for the source ; ω its pulse), of the width σ and center t m located at the maximum of the packet.

2. Fit of the source using formula 2. 23. This fit measures the pulse ω, the value provides a measure of the actually frequency and the phase at origin. These remain fixed in the fitting of the reception signal (bottom graph). Note the fit with the non-dispersive expression of the response (assume c φ = c).

These fit values are of course all saved and we compute the group 8 and the phase velocity 9 . respectively from H c and H c φ (respectively the flight time of the wave packet and of the carrier wave).

The phase of a sine wave being defined at modulo 2π, the phase shift corresponding to -ω H c φ is therefore defined modulo 2π. Therefore, in order to properly 8 Group velocity The group velocity of a wave corresponds to the velocity at which its overall shape (i.e. its modulation, its envelope) propagates. In our system, it is the velocity travel of the Gaussian-like packet. More generally, the group velocity is defined as:

c = dk dω
k and ω being respectively the wavenumber and the angular frequency. We have computed the group velocity using the measurement of the flight time (i.e. the temporal shift) between the generation of the wave and its arrival at the receptor. 9 Phase velocity The phase velocity is defined as:

c φ = k ω
Physically, it corresponds to the celerity at which the information transmitted by the signal (also called the carrier wave) propagates. It thus corresponds to the velocity at which a "phase event" travels through the medium.

DATA PROCESSING AND SELECTION

compute the phase velocity we need to unfold the phase. Once the phase has been unfold according to this continuity criterion, we can properly compute the phase velocity as:

c φ = H δφ (2.24)
This computation will be useful to check if c ∼ c φ , i.e. whether or not the non-dispersivity we assumed in the first place is a valid assumption and coherent with the theory we developed.

Selection of the wave packets

In this experiment, many packets are transmitted (during each campaign, over 40000 wave packets), but the sources of noise and other perturbations are numerous. It is therefore important to be able to evaluate the quality of the fits with respect to the data and then to check the validity of the packets we considered with respect to the theory we want to probe.

Success of the fit

We define a the following quantity χ 2 in order to measure the quality of the fit, on a purely mathematical point of view:

χ 2 = y signal -y fit max (y signal ) 2 (2.25)
For each wave packet, a value of χ 2 is therefore associated. The smaller, the more accurate the fit. The fit is computed within a time range [t 20% ; t 30% ] where t 20%

is the instant -before the max of the wave packet-after which the amplitude of the signal is larger than 20% of the max of the wave packet (similarly, t 30% is the instant -after the max of the wave packet-before which the amplitude of the signal is larger than 30% of the max of the wave packet). One can choose to evaluate this χ 2 over the same time range as the one over which we compute the fit. Thought some wave packets exhibit a "twin" (such as the one figure 2.15). They can be the consequence of some reflection of the wave packet against the lid and the bottom plate. These packets must be discriminated in order to remove them in the selection of the packet we keep for further treatment. Therefore we choose to compute the χ 2 value over an extended time range: [t 20% ; t 30% + 5σ] (σ being the width of the Gaussian envelope). Also, the ratio signal/noise is sometimes rather small, and therefore the fit does not capture the wave packet. Extending the time range over which χ 2 is computed enables to also remove the packets for which this is the case.

In the following, we will only mention χ 2 , and it will only refer to its definition over the extended time range. Figure 2.17 presents the scatter plot of the group velocity of the wave packets versus the confining pressure at the top of the box.

Maker color characterizes the value of χ 2 . Changing the colorbar range allows to distinguish the trend characterized by the more successful fits (in dark blue). We hence define a threshold for χ 2 and select the wave packets for which χ 2 < 0.04.

Defining an acoustic limit

Once we made sure we kept in our selection only the wave packets well described by formula 2.21, we can use the actual measurements we get to dig into the physics. The validity of the packets is then evaluated with respect to the hypothesis we make within the theoretical framework. First, we only consider the case on linear propagation.

DATA PROCESSING AND SELECTION

This implies that wave packets do not change shape during the propagation through the medium. Therefore, by computing the distortion, one can select and extract the very packets that are representative of a linear propagation (i.e. for which the distortion is close to zero).

The Distortion D σ . It characterizes the temporal stretching of the wave during the propagation in the medium. We define it by the ratio of the respective widths of the Gaussian envelopes:

D σ = σ Reception σ Source -1 (2.26)
D σ = 0 means the propagation is linear as the wave packet is transmitted without any deformation. We choose a threshold distortion of 20% as a selection criterion for the wave packets we keep.

Small amplitudes. All the theoretical development we make are in the framework of small perturbations. We therefore use the selected wave packets to extrapolate the value of the group velocity at source amplitudes A sr c → 0m/s 2 . To do so, we gather wave packets by typical values of P top within a range of ±5P a around each focus values of P. Then, for each pressure, we make a linear fit of the group velocity V g versus the source amplitude A sr c and we extract the value of V g(A sr c → 0) (see an example for wave packets at 4kHz, 180P a, figure 2.19).

We then have a dataset of group velocities for a range of pressure, at an acoustic limit defined by small (zero) source amplitudes of linearly propagated wave packets.

We use the same process to compute the phase velocity at the acoustic limit (results are displayed figure 2.21).

On figure 2.22 we superpose the group and phase velocity hence obtained. We see that the phase velocity is quite superior to the group velocity, especially at highest confining pressures. However, we need to consider this discrepancy with respect to errors we have on the velocity measurements, such as regarding the fit of the theoretical wavepacket formula, the precision of our measure of the cell height (which

Granular contacts in jammed weakly confined granular media

Mindlin's group has pioneered on the mechanics of granular system under compression. Using Hertz contact law, they have predicted a P 1/3 dependence of the elastic modulus. However many experimental data (among which [45]) seem to give a different power law: E ∝ P 1/2 . In 1996 De Gennes publishes a paper [38] in which he presents a model of grain contact involving a soft shell on the outer surface. Using the same arguments as before, using the contact area and the estimation of the pressure, he predicts the effective elastic modulus of such a sphere assembly to be proportional to P 1/2 at low confining pressure. The isotropic pressure in the effective medium is defined by:

P = - ∂F ∂V (2.27)
So here it satisfies:

P 4πR 2 dδ φ = Z dF c (2.28)
where F c is the free energy of one contact and φ is the volume fraction of grains.

Therefore the pressure is:

P = Zφ 4πR 2 dF c dδ (2.29)
Notice that the normal force on one grain is given by:

F N = dF c dδ (2.30)
Also, the compressibility is defined by:

K -1 = - 1 V dV dP (2.31)
Here, V = 4 3 π (Rδ) 3 for one grain squeezed by a displacement δ. So:

dV V = - 3 dδ R (2.32)
The sound celerity is given by:

c = R dP dδ 3ρ (2.33)
M where ρ is the density of the effective media thus ρ = ρ g φ. Therefore we have an implicite relation between the speed of sound and the pressure P. One can plot this relation and find an equivalent empiric formula that will be used to fit the data.

Indeed, we can express the pressure with respect to δ: The fit will therefore provide the independent measurement of R/δ c and Zφ.

At small δ/δ c , the expressions simplify into: (2.42)

P P = 6 δ δ c 2 (2.
Conversely, at large δ/δ c , the expressions simplify into: (2.44)

P P = 5 δ δ c 3
It is interesting to identify the elastic modulus E using the mean field calculation.

Starting from the normal force between grains:

f = 8 3 µ g 1 -ν g R 2 δ R 3/2
(2.45) one derives the mean field bulk modulus:

K = ρc 2 = 1 3π
µ g 1 -ν g Zφ (2. 46) where µ g is the shear modulus and ν g the Poisson ratio. From which we get: (2.50)

E = 16 √ 2 15 µ g 1 -ν g (2.
A Taylor expansion yields in the limit of P → 0: The experiment seems to show that, even down to very low confining pressures, the Hertz scaling does not match the data. We show that the P 1/4 scaling relation is very robust. We provide an interpretation consistent on data in higher pressures

c ∝ R 1/4 P 1/4 . ( 2

SOFT SUPERFICIAL LAYER

Mindlin's group has pioneered on the mechanics of granular system under compression. Using Hertz contact law, they have predicted a P 1/3 dependence of the elastic modulus. However many experimental data (among which [45]) seem to give a different power law: E ∝ P 1/2 . In 1996 De Gennes publishes a paper [38] in which he presents a model of grain contact involving a soft shell on the outer surface. Using the same coarse-grained approach as before, and a contact model in which only the shell of the grain deforms, he predicts the effective elastic modulus of such a sphere assembly to be proportional to P 1/2 at low confining pressure.

Soft superficial layer

In the following, we present a theoretical framework which aims at describing the mechanical properties of granular packing from high to weakly confining pressures.

We consider a grain in contact with an infinite plane by a normal force F N . It has a soft superficial layer of thickness e (see figure 3.1 page 181). The Hertz analysis gives the contact area A = πRδ where δ is the overlap between particles and R the radius. The bulk free energy takes dimensionally the form:

F b ∼ E(Rδ) 3/2 δ 2 b
Rδ while the superficial energy reads:

F s ∼ ηERδe (δ -δ b ) 2 e 2
η 1 is the elasticity contrast (i.e. the ratio between the shell and the bulk elastic moduli). Minimizing the free energy with respect to δ b , we get: The isotropic pressure in the effective medium is defined by:

δ b = δ 1 + δc δ
P = - ∂F ∂V (3.1)
So here it satisfies:

P 4πR 2 dδ φ = Z dF c (3.2)
where F c is the free energy of one contact, φ is the volume fraction of grains and Z the average number of contact per grain in the medium. Therefore the pressure is:

P = Zφ 4πR 2 dF c dδ (3.3)

SOFT SUPERFICIAL LAYER

Notice that the normal force on one grain is given by:

F N = dF c dδ (3.4)
Also, the compressibility is defined by:

K -1 = - 1 V dV dP (3.5)
Here, V = 4 3 π (Rδ) 3 for one grain squeezed by a displacement δ. So:

dV V = - 3 dδ R (3.6)
The bulk modulus is given by:

K = R dP dδ 3 (3.7)
and assume the proportionality between the speed of P-waves c ∝ K/ρ where ρ is the density of the effective media thus ρ = ρ g φ. Therefore we have an implicit relation between the speed of sound and the pressure P. As relations are defined within a multiplicative constant, we can absorb this constant into the elastic modulus E, which remains to be identified. We can express the pressure with respect to δ: The fit will therefore provide the independent measurement of R/δ c and Zφ.

At small δ/δ c , the expressions simplify into: We can now rewrite the expression of the celerity (see formula 3.19): (3.28)

c = P * ρ   

Fit of the model

In the latter expression of the velocity, P * depends on the grain material (µ g ; ν g ), configuration (coordination number Z) and on the grain dimension (radius R). Also, the parameter δ c reflects characteristics about the soft superficial layer (its thickness e and its elastic modulus ηK bulk ). It is therefore hardly possible to compare our experimental data with the one obtained by Makse and Domenico, as we have no information about the genuine nature of the grains (roughness of their surface 1 ...) and on all the experimental parameters that were not directly controlled but which can still have an influence on the nature of the contact (for instance the humidity rate in the room, hence the condensation at the surface of the grains). seem in rather good agreement between Makse's and Domenico's. The numerical data were obtained by 3D simulation of Hertzian frictional packings 2 . The simulation seems also in agreement with the experimental results, at least for higher pressures.

Below ∼ 10 7 P a, the numerics start to differ from the experiments. Authors invoke a loss of the signal and experimental difficulties at lower pressures.

We hence propose to use our data (which were made at much lower pressures, at least six orders of magnitude from Domenico's and Makse's) to compare with the prediction from the numerics. In order to fit our experimental results (which were made at confining pressures about six orders of magnitude smaller in our experiment), we proceed in two step:

in the first place, fitting of the parameter P * , at high pressures, and then lock this value of P * and fit the remaining parameter, δ c /R. While P * characterizes the bulk material properties of a Hertzian packing of beads, the parameter δ c /R is determined by the mechanical response of the packing at vanishing pressures, in the regime in which the deformation of the soft superficial layer is dominant.

Thanks to Domenico [46] (see left hand side of expression 3.24), we can use the numerical data from [48] to compute the celerity for the P-waves, in the case of Hertzian contacts (we recall here that Makse's simulations use a Hertzian potential of interaction between the particles). These values are plotted figure 3.3 page 188, with respect of the confining pressure. We then fit the first parameter of expression 3.28 with these numerical data for the lower pressures available and find P * = 41GP a.

This is the blue dashed line (it is a power law, following the Hertz prediction) on figure 3.3. The reader should note the importance on fitting the parameter P * only at lower pressures: P * depends of the average contact per grain in the packing Z, which is a function of the confining pressure P . Following numerical simulations by Makse [48], it seems that Z(P ) saturates at a minimum value for low confining pressures, for frictional 3D packings. We therefore consider P * to be constant with respect to the pressure range we consider. Then we lock this value into the fit formula and fit the parameter δ c /R (this is the continuous red line on figure 3.3).

We measure δ c /R = 3.10 -5 .

The reader should note that P * is not the value of the pressure at which crossover between the two regimes occurs. However thanks to the values we got from the fit of the model, we can estimate the cross-over at P 3500P a, which seem to be consistent with the plot figure 3.3.

Conclusion the experiment

In order to characterize the nature of the soft superficial layer we modeled, we made Atomic Force Microscopy (AFM) measurements of the surface of the glass beads we used in our main experiments. The beads we consider have a millimetric diameter, therefore, from the fitted In order to have a more quantitative understanding of the elasticity contrast with respect to the nature of the asperities of the grains, we need to look at the power density spectrum of these roughness measurements. Also, it would be interesting to simulate the effective contact area between the surfaces we measured by AFM. The work is currently in progress. The plate in contact with the granular media is fixed to a metallic frame (figure a) which motion is controlled by a piezo-electric actuator and chrysocal blades in order to create a restoring force to maintain the contact with the actuator that only has a pushing work power. A setting screw (not visible here) enables to properly set the pre-stress on this spring system. ing. Indeed, shear waves are expected to travel at lower velocities that compression waves. Chapter 5

Conclusion on the acoustic in weakly confined granular media

The present study aims at measuring the mechanical response of granular packings as close as possible from the jamming point. For the sake of experimental concerns, we chose to use the confining pressure as the main control parameter of our experiment.

The question of whether it is possible or not to reach the jamming point by only decreasing the pressure is controversial. It seems that the friction within the packing locks the configuration, even at vanishing pressures. The sample is very historydependent for hysteresis due precisely to friction in between the grains. We tried to initialize our preparation before each flight, passing a grid through the whole sample and by vibrating it for a finite time with a sine wave, but there is certainly aging in the system, especially because of the successive 0 to 1.8G phases inherent to the parabolic flights.

However, we have shown that going to vanishing pressures ensures to reach a regime at which the mechanical response of the packing mainly relies on friction, i.e. asperities at the surface of the grains. We proposed a theoretical model that enables to qualitatively understand the nature of the contacts at such pressures. We propose further developments that are currently on progress in order to have a more quantitative understanding of the model and on the fit parameters it provides. being the width of the cell (which we know because we choose it), the total surface occupied by the particles, with respect to y if given by the product of W and the integral of the volume fraction with respect to y . Then, by dividing this quantity by the average area for one grain (πd 2 /4), one gets the count of particle with respect to y .

This allows one to localize the physical position of the bottom of the numerical cell and therefore to properly measure the height of grains in the numerical cell.

0. 1 .

 1 STATES OF MATTER, CHANGE OF STATES 0.1 States of matter, change of states 0.1.

Figure 1 :

 1 Figure 1: Typical P -T phase diagram. The dashed line is the boundary between the solid and the liquid phase for water and bismuth. A is the triple point at which the three phases coexist ; C is the critical point of the system.

FIGURE 1 . 6 -

 16 FIGURE 1.6 -Représentation schématique de la fonction de corrélation de paire pour un solide cristallin idéal (a), et pour un fluide simple homogène (b).

Figure 2 :Figure 2

 22 Figure 2: Typical pair correlation function for (a) a solid crystal and (b) a homogeneous liquid. The solid exhibits a structural order at long range, whereas the liquid does not. For liquids, the pair correlation function tends to 1, uniform probability (Figure extracted from [1])

Figure 3 :

 3 Figure 3: Pair correlation function in respect of the temperature. When the temperature increases, the order disappears. Simulation of CCP packing of Nickel. Data from lecture of Matthieu Micoulaut (UPMC).

Figure 4 :

 4 Figure 4: Typical stress-strain curve.In green, the cycle is reversible. Stress does not rise over the elastic limit σ e . In red, the cycle is irreversible as the stress goes over σ e : the material is plastically deformed and the sample exhibits a remnant deformation.

Figure 5 :

 5 Figure 5: Typical fluid behaviors. Shear stress in respect to the shear rate.

Figure 5

 5 Figure 5 presents different common fluid behaviors. Diverging from the Newtonian fluid, shear thickening and shear thinning respectively describe a hardening or a easing of the shearing of the fluid for a given shear rate. Some fluids also present a yield under which the shear stress is not sufficient to deform the material. For all these material, the viscosity therefore varies with the shear rate. Another very interesting fluid behavior is the thixotropy, defined as the dependence of the fluid viscosity on the loading time. Example of thiroxotropic fluids are gels, concrete and other liquids containing solid inclusions, etc. Research in this field is

Fig. 1 .

 1 Fig. 1. Experimental and numerical results for soft and granular matter displaying non-locality. (a) Velocity profile of the flow of a jammed emulsion in a micro-channel with rough surfaces, data from [5]. Local rheology would predict a plug flow. (b) Velocity profile of a granular flow in a 2D shear cell. Local rheology would predict a linear velocity profile. The numerical data are from [37]. (c) Experimental (triangles) and numerical (squares) velocity profiles of a foam flow in a Couette cell, data from[38] and[72], respectively. Local theory would predict a localised failure at wall.

Fig. 3 .

 3 Fig. 3. Image sequence showing the temporal evolution of the local shear rate Γ for a system of frictional hard grains (kn/P = 2 • 10 4 ) in the quasi-static limit (I = 5 • 10 -4 ). Similar results are found for frictionless particles. Color code from blue ( Γ = -40 γ) to red ( Γ = 40 γ). Time lapse between two successive images: γ∆t = 10 -3 .

Fig. 4 .

 4 Fig. 4. Image sequence showing the temporal evolution of the local shear rate Γ for a system of frictional soft grains (kn/P = 10) in the quasi-static limit (I = 5 • 10 -4 ). Similar results are found for frictionless particles. Color code from blue ( Γ = -40 γ) to red ( Γ = 40 γ). Time lapse between two successive images: γ∆t = 10 -3 .

Fig. 5 .

 5 Fig. 5. Space-time diagrams showing the local contribution of Γ to the shear rate γ, measured on the central line of the cell, for a system of sheared hard (a) and soft (b) grains. For this example, the ratio of the grain contact stiffness k n to the overall pressure P is 2 • 10 4 and 10 respectively; the shear rate corresponds to an inertial number I = 5•10 -4 . Color code from blue ( Γ = -40 γ) to red ( Γ = 40 γ). (c) Probability distribution function (PDF) over time of the space average ⟨ Γ ⟩, for hard (black line) and soft (red line) grains. (d) PDF of the spatial standard deviation δ Γ . Red dotted line: δ Γ computed when |⟨ Γ ⟩| ≥ 5 γ (soft system).

Fig. 6 .

 6 Fig.6. Relaxation length L for frictionless (a) and frictional grains (c) below (red and yellow circles) and above (blue and green squares) yield conditions, data from[37]. Solid lines: fit of the data to eqs. (29) and (35), diverging as |Y -1| -1/2 when Y → 1. (b), (d) Log-log plot of the same quantities.

  Figure 2.1 presents different averaging of the same film in which we monitored the flow in the narrow channel we study in the following. Each image corresponds to a typical exposure time as it is computed from different averaging times from the same image acquisition.

  exp. 0.025 sec exp. 0.25 sec exp. 1.5 sec

Figure 2 . 1 :

 21 Figure2.1: By changing the exposure time of the photo shot, we can see that it is not obvious that a definite layer of grains flows. Indeed, there is a so-called "creep motion" in the depth, even below the critical depth under which the yield parameter Y < 1.

Figure 2 . 2 :

 22 Figure 2.2: Photo of the narrow channel experiment. The optical axis of the Fastcamera is perpendicular to the main axis of the channel. During the experiment, stabilized yellow LED spots lights the observation windows on the channel wall.Beads are all of the same kind (same reference from Silibeads) but in two different colors in order to increase black and white contrasts on the images.

1

  Silibeads references: 4504-288 LS for the blue and 4504-110 LS for the yellow beads.

Figure 2 . 3 :

 23 Figure 2.3: Flowing beads in the narrow channel. Side view. This picture has been shoot with an SLR camera. The conversion in gray scale (right) shows an equivalence of the contrast we get from the fast camera (Phantom M340).

Figure 2 . 4 :

 24 Figure 2.4: Reservoir and flow-rate control. The system is divided into three parts: the main reservoir feeds the intermediate cylindric bin within which the Beverloo law is verified. The outlet of the overall system is removable so we can adjust the tube diameter D.

Figure 2 . 5 :

 25 Figure 2.5: Measure of the flow rate (mass flow rate) with respect to the aperture diameter D. Fit with the Beverloo law [18]. On the graph, we fit with an expression a + b D and find a = -0.001864 ± 0.000163 and b = 0.85199 ± 0.0118

64

 64

1 .

 1 image acquisition (and measurement of the live mass flow rate) 2. compute an image of the average intensity of the whole film 3. free surface detection by a wavelet algorithm (convoluate each vertical line of the image, transverse to the direction of the flow, with a wavelet, and take the maximum). Fit this curve by a line of equation y = ax + b.

4 .

 4 rotates all the images by an angle θ r otate = tan -1 a b .5. compute the velocity profile (see explanation of the algorithm section 2.2.2)

Figure 2 . 7 :

 27 Figure 2.7: Localization of the free surface. Superposition of an image average (over the whole film) and the corresponding intensity profile. Images are on 12bits. Relative intensity is computed so that the maximum possible value for a pixel 2 12 ≡ 100%.

Figure 2 . 8 :

 28 Figure 2.8: Raw velocity profiles (length are in pixels). Semilog scale. Crosses show the velocity at the position of the free surface.

2. 4

 4 Continuous analytical model for dense granular avalanches2.4.

11

  Using σ xz = η(z) dux dz = µ(I) P (z) we can show that the effective viscosity of the medium η(z) ≡ µ(I) I d P (z) ρ p .

Figure 2 . 10 :

 210 Figure 2.10: tan θ vs dimensionless flow rate Q. Fit of relation 2.5.3. We find µ c = 0.41353 ± 0.00303 and B = 0.22687 ± 0.00982

Figure 2 .

 2 Figure 2.12 shows the linear fit of h (in SI base units) versus h. Using the relation h = L h we then have access to the characteristic length L = W µ c µ w , which, with the previous determination of µ c (and W = 31.3mm), gives an estimation of the friction coefficient of the glass walls: µ w = 0.174 ± 0.004.

Figure 2 . 11 :

 211 Figure 2.11: Velocity profiles -side view. Fit with the local rheology (dashed lines) for several flow rates. From purple to red, mass flow rates are: Q m = 130.2 -95.46 -74.50 -48.61 -29.07 -16.61 -9.435 -4.994g/s (±0.04%). We get V a L5/2 = 12102 ± 1 and respectively h = 0.0185116 -0.0160965 -0.0148754 -0.0127977 -0.0110859 -0.00911331 -0.00686528 -0.00640506 (in m).

Figure 2 . 12 :

 212 Figure 2.12: Dimensional flowing height h versus dimensionless height h. From definition 2.23 and relation 2.21, one can measure the characteristic length L.

Figure 2 . 13 :

 213 Figure 2.13: Relative error made with the plug flow assumption. Q represents the dimensionless flow rate measured from the weighting scale.

Figure 2 . 14 :

 214 Figure 2.14: Velocity profile measure from top, for several flow rates. From purple to red, mass flow rates are: Q m = 130.2 -95.46 -74.50 -48.61 -29.07 -16.61 -9.435 -4.994g/s (±0.04%).

Figure 2 . 15 :

 215 Figure 2.15: Velocity profiles in semilog scale-side view (see same graph in linear scale figure 2.11 page 78). The fit with the local rheology (dashed lines) predicts a zero velocity below at depth h(θ). It therefore cannot describe the flow below this limit. From purple to red, mass flow rates are: Q m = 130.2 -95.46 -74.50 -48.61 -29.07 -16.61 -9.435 -4.994g/s (±0.04%).

CHAPTER 2 .

 2 NON-LOCAL RHEOLOGY: EXPERIMENTAL APPROACH 81

2 ∇ 2 II

 2 . The constitutive relation can be expanded around the relation Y = µ(I), valid in the homogeneous case according to

Figure 2 . 16 :

 216 Figure 2.16: Asymptotic fit of the tails of the velocity profiles for aperture diameters ranging from 12 to 20mm, i.e for several flow rates. From purple to light green, mass flow rates are: Q m = 130.2 -95.46 -74.50 -48.61 -29.07g/s (±0.04%).

  figure2.17.

  figure2.17.

Figure 2 . 17 :

 217 Figure 2.17: Velocity profile measure from top, for several flow rates. Fit with a phenomenological cosh profile. From purple to red, mass flow rates are: Q m = 130.2 -95.46 -74.50 -48.61 -29.07 -16.61 -9.435 -4.994g/s (±0.04%).

Figure 3 . 1 Figure 3 . 1 :

 3131 Figure 3.1 shows experimental measures on the incline plane. We clearly see the

3. 1 .

 1 NUMERICAL SIMULATIONS OF GRANULAR MEDIA otherwise), and k t = k n /2 following [32].

4

 4 and a softness parameter S = 3.10 3 . These considerations are very important regarding the simulations we run, espe-CHAPTER 3. NON-LOCAL RHEOLOGY: NUMERICAL APPROACH 101

Figure 3 . 3 :

 33 Figure 3.3: Evolution of the critical friction coefficient µ c with respect to the softness parameter S = k n /P . Measures are made at fixed shear rate γ m P d and fixed pressure P = 1. Data from [1].

Figure 3 . 4 :

 34 Figure3.4: Evolution of the kinetic energy (a) and the total number of contact in the system (b) with respect to dimensionless time ( γ-1 being the characteristic time of the deformation), at the onset of the flow (example of the incline plane: from a preparation at vertical gravity, we incline the latter by an angle θ = 0.2. Frictionless case.)

Figure 3 . 6 :Figure 3 . 7 :Y < 1 .

 36371 Figure 3.6: Top wall position with respect to time. Fluctuations of the height of the cell is negligible with respect to the grain diameter d.

3. 2 .

 2 TEST CASE: NON-LOCAL RHEOLOGY IN A PLANE SHEAR CELL As seen in the first chapter, the non-local rheology predicts velocity profiles of form: u P/ρ g = I ∞ z + δI sinh(z/L) (3.12) The deviation of the velocity profile from the local contribution is characterized by the term δI sinh(z/L) that measure the non-locality of the flow. Fitting the profile by the non-local model hence means in particular to measure these relaxation length L with respect to the control parameter Y.

Figure 3 . 8 :

 38 Figure3.8: Setup for the plane shear cell. Pressure is imposed, so as the shear strain (wall velocity). This is a snapshot of a simulation of hard frictionless particles. Colors correspond to the instantaneous velocity along axis x.

Figure 3 .

 3 Figure3.9b. shows that the volume fraction is homogeneous in the bulk, so as the stresses P and τ (figure3.9c). and the stress parameter τ /P (figure3.9d).

Fig- ure 3 .Figure 3 . 9 :

 339 Figure 3.9: Mean profiles for the custom plane shear cell, for a run at Y > 1. (a) Velocity profile with respect to the height z/d. Black line is the result from simulation. Green is the prediction of the local rheology, red is the fit of the nonlocal rheology. (b) Corresponding profile for the volume fraction φ. (c) Stress state. Normal (pressure P ) and shear stress τ with respect to the height z/d. (d) Profile of the dynamical friction coefficient µ = τ /P with respect to the height z/d.

Figure 3 . 10 :

 310 Figure 3.10: Mean profiles for the plane shear cell, for a run at Y < 1. Buffer zones are close to the walls, in yellow. (a) Velocity profile with respect to the height z /d. (b) Corresponding profile for the volume fraction φ. (c) Stress state. Normal (pressure P ) and shear stress τ with respect to the height z/d. (d) Profile of the dynamical friction coefficient µ = τ /P with respect to the height z/d.

) 110 3 . 2 .

 32 TEST CASE: NON-LOCAL RHEOLOGY IN A PLANE SHEAR CELL with L the relaxation length. According to whether the system is above (Y > 1) or below (Y < 1) the flow threshold, we have:

Figure 3 . 11 :

 311 Figure3.11: Effective friction coefficient µ versus I ∞ . Data points are obtained from the global fitting on all our runs. Solid line is the fit of the phenomenological law µ(I), in which the critical friction coefficient is fixed at µ c = 0.0953 obtained from global fit. We get a = 5.25. Note that this result is slightly differs from the result from the global fit.

3. 2 .Figure 3 . 12 :

 2312 Figure 3.12: Relaxation lengths L relative to the yield parameter in the bulk. Values of L are obtained from the global fit.

Figure 3 . 13 :

 313 Figure3.13: Setup of the incline plane. Gravity is inclined by an angle θ. Control parameters of the numerical experiment are θ and the microscopic friction between particles (µ p = 0 or 0.5). On left and right sides of the cell, boundary conditions are periodic, and the bottom of the cell is made of fixed grains (in black) in order to create a rough/no slip boundary condition. On this graph, colors characterize the velocity of each particle along the main flow direction x.

20 )Figure 3 .Figure 3 . 14 :Figure 3 . 15 :

 203314315 Figure 3.14 presents the stress profiles and the effective friction coefficient µ = τ /P along depth z, for a run at θ = 0.2. As specified on the figure, z axis is downward and its origin z = 0 localizes the free surface. The linear fit of the stresses, from 3.20 also provides a value of ρg, where ρ = ρ g φ. In our code, g = 1 and the average mass of one grain of diameter d (taken at d = 1). The system is 2D therefore

Figure 3 . 16 :

 316 Figure 3.16: Incline plane. Velocity profiles for θ = 0.2 for several flowing heights h (all runs from h ≈ 30.3 and then -50 grains between each consecutive run). Frictionless case.

Figure 3 . 18 :

 318 Figure3.18: Velocity profile on the incline plane (solid line). Prediction of the local rheology (dashed line). Inset is the same graph but in semi-logarithmic scale. For this runs, θ = 0.2 and h ≈ 30.3.

Figure 3 .

 3 Figure 3.18 presents a typical velocity profile obtained for a simulation of a steady flow on the incline plane, for a height of flowing grains ≈ 30.3 and an inclination angle θ = 0.2. In dashed line is the prediction from the local rheology, using the

Figure 3 . 19 :

 319 Figure 3.19: Incline plane. Inertial number with respect to depth. All runs at same θ = 0.2.Dots : numerical data from simulations. Lines: cofit of all runs with the linear model. For this particular inclination (θ = 0.2), we find I b = 0.0455 and L = 1.54. Frictionless case (µ p = 0). Runs for heights H/d = 28.91 -22.01 -15.15 -10.95. Inset is the same plot but in semilog scale. Top graph is the residuals of the cofits.

. 31 )

 31 with L the same for all heights of grains H at fixed θ and therefore fixed Y b .This expression of I(z) is the solution of a linearized equation (equation3.28). It is thus expected to be valid around the base state I b .

3. 3 .Figure 3 . 20 :

 3320 Figure 3.20: Inertial number profiles for different θ at a given depth.

Figure 3 . 21 :

 321 Figure 3.21: Global fit results on inertial profiles at same yield parameter (θ = 0.2) and for several heights of grains (H = 28.91 -22.01 -15.15 -10.95). Inset is the same plot in semi-logarithmic scale. Top graph are the fit residuals.

Figure 3 . 22 :

 322 Figure3.22: We report the fitted relaxation length L measured from the fitting of the inertial profiles on the incline plane and report them on the calibration from the plane shear cell. Corresponding values of Y are computed from its definition (Y = τ /(µ c P ) using µ c = 0.0953 we obtained from the previous calibration, section3.2.3).

  is a function of both Y and , hence of the rheology. The value we obtain should therefore be consistent with the test case of the shear cell. On figure 3.22 we report the values of L we get from the global fit. From the consistency of the values of I b we found (as attested by figure 3.21) in respect to the yield stress parameter values we probe, the discrepancy of the fitted values of L particlularly raises the question of the validity of the adjustment of the non-local parameter . We recall that these cofit are done using the solution (3.31) of the linearized equation (3.28). Therefore this is not the exact solution of the non-local equation(3.27).

Figure 3 . 23 :

 323 Figure 3.23: Estimation of the value of the inertial number at the free surface with respect to the distance from the yield parameter Y. Colors are for three different height of grains: H/d = 28.91-22.01-15.15. Dashed line is a linear fit of expression y = a + bx. We find a = -0.007895 ± 0.00196 and b = 0.1623 ± 0.0023

126 3 . 4 . 1

 341 CONCLUSION AND PERSPECTIVES accurate numerical solution of the whole non-local model. For further investigation, we prescribe to use a shooting method in order to solve the non-local equation and hence to accurately measure the parameters of the model. Doing so would also provide a way to better define the boundary condition at the free surface and to relate it to the yield parameter. Elasticity of jammed granular matter . . . . . . . . . . . . . 1.1.1 Elastic moduli in the vicinity of the jamming transition . . 1.1.2 Mean field approach . . . . . . . . . . . . . . . . . . . . 129 130 1.1. ELASTICITY OF JAMMED GRANULAR MATTER 1.1 Elasticity of jammed granular matter 1.1.1 Elastic moduli in the vicinity of the jamming transition

Figure 12 .

 12 Figure12. Bulk (K ) and shear (G) modulus as a function of distance to jamming for two-dimensional bidisperse systems, with interaction potential V ∼ δ α (see equation (2)). The closed symbols denote moduli calculated by forcing the particles to move affinely and the open symbols correspond to the moduli calculated after the system has relaxed. Slopes are as indicated (adapted from[2] with permission-copyright by the American Physical Society).

Figure 1 . 2 :

 12 Figure 1.2: Compressive wave velocity c with respect to the confining pressure P . In dashed line are the prediction from EMT (blue tight dashed line is with Z = 5, = 1, green dashed line is with Z = 5, = 0); range squares are experimental results from Domenico and red circles are experimental results from Makse. Black solid line is a phenomenological curve in P 1/4 . Data extracted from [47].

Figure 2 . 1 :

 21 Figure 2.1: The two racks of the experiment. The additional screen and the PC desktop unit were used for another experiment we left aside for the moment.

2. 1 . 2

 12 The cell and its confining systemThe cell contains glass beads with a diameter range 1 < d < 1.3m from Silibeads. It is a parallelepiped box of height h = 11cm and of section dimensions W = 15.5cm, L = 16cm and fill with about ∼ 3, 3kg of grains (depending on the flight campaign).The inner walls and the top of the box are made of 3D-printed ABS. Their external surface is plane but the core is a honey nest: this 3D printed structure, in addition to the ABS material used, insures a phonic insulation from the exterior of the box (no acoustic by-pass) and also prevents reflexion of the sound wave studied within the box. The bottom of the box is a metallic plate mounted on three synchronized piezoelectric pistons (P 840.10 from PI) producing a tunable vertical oscillating motion at frequencies up to several kHz.

Figure 2 . 2 :

 22 Figure 2.2: CAD model of the test cell "BoxSon". One can distinguish the inner walls and the lid 3D printed in ABS. The remaining parts, among which the vibrating plate and the adjusted square around it are made of aluminum alloy.

  To measure the wave packets we transmit through the medium, we use two Brüel & Kjaer high frequency accelerometers (bandwidth 20Hz -100kHz). The source accelerometer, mounted underneath the vibrating plate is a 4705-B (unidirectional) ; the reception accelerometer is a 4735-B (3D). There sensibility is 1mV/G and their factory (and supposedly 2.1. EXPERIMENTAL SET-UP

Figure 2 . 4 :Figure 2 . 5 :

 2425 Figure 2.4: Jitter remnant acceleration during a parabola. Pressure is imposed on top of the cell lid. a) Schematic of the cell, with confining pressure from at the top σ 0 ; b) Effective gravity (on the three axis) and pressure at the top σ 0 and at the bottom σ B of the cell, with respect to time. The pressure at the top is controlled by a feed-back loop in order to maintain σ B at the desired value. c) Acceleration signal at the source and at the reception, for one transmitted wave packet. The measure of the flight time τ f provides a value for the group velocity at which the wave travels.

Figure 2 . 6 :

 26 Figure 2.6: Current g and adjustment of the confining pressure by the feedback loop during one parabola. Top graph: relative gravity (g 0 = 9.81m/s 2 ) felt in the aircraft, along a parabola. Center graph: vertical g (along axis x) in the plane and position of the linear fine stage with respect to it. Bottom graph: direct measure of the confining pressure σ 0 (at the top of the cell) and regulated pressure at the bottom, σ b .

Figure 2 . 7 :

 27 Figure 2.7: Experimental data: source and reception.

We developed a

  LabView program to control the experiment and compute the feedback loop to adjust the pressure in real time. The control parameters of the experiment are loaded from a text file previously prepared. On each parabole, i.e. each line of the ParabolaConfig.txt file, we choose the target pressure for the retroaction, the carrier wave frequency, the frequency at which we send the wave packets (which we maintain at 200Hz), and the amplitude of the generation of the wave packets. We probe 3 different source amplitudes for each parabola as the code was made to load and generate identical series of 3 wave packets. During each parabola we propagate about 300 wave packets, therefore ∼ 100 wave packet for each imposed amplitude.The ParabolaConfig.txt file is prepared prior to the flight as the ZeroG conditions are not viable for clear thinking and smart estimations of the setting for the control parameters of the experiment.

Figure 2 . 8 :

 28 Figure 2.8: User interface developed for the parabolic flight experiment.

Figure 2 . 8

 28 Figure 2.8 presents the user interface of the experiment. The experimentalist chooses the parabola number (on the left of the panel), loads the parabola parameters by clicking on the "select" button. During the loading of the parameters, the larger linear stage adjusts its position to an estimated value corresponding to the target pressure specified in the text file. If the retroaction is enabled (which we do in normal situation, no debbugging etc.), the finer linear stage will adjust its position during the zeroG phase. At the bottom right of the window we visulamize the estimated pressure at the top and at the bottom of the cell. Pressure at the top (P top )

Figure 2 . 9 :

 29 Figure 2.9: Setup of the source accelerometer. It is fixed to the vibrating plate, on the opposite side from the granular medium. z-axis of the accelerometer is oriented toward the top to the bottom of the cell.

Figure 2 . 10 :

 210 Figure 2.10: Setup of the reception accelerometer.It is encrusted within the lid, in a silicon bed, and it comes out a little from the bottom surface of the lid so we make sure it is in contact with the grains. z-axis of the accelerometer is oriented toward the top to the bottom of the cell. Note the polymer cubes that ensure the damping of between Linear stage and the cell lid, in order to prevent rigid contact and therefore an acoustic bridge with the medium.

Figure 2 .

 2 Figure2.11: We shock the accelerometer in its compartment in the lid. The system is free to oscillated and beside the first shock pulse, received no other excitation. Red is the time response and green is the theoretical prediction.

(

  time, vertical source acceleration and the 3 axis of the reception acceleration), the time origin for this parabola and the value of the control parameters in two other .txt files. The synchronization lies in dividing the permanent recordings into the series of parabola. It creates a data.mat file within each parabola folder that contains synchronized signals of force measurements and accelerations, with time. It also contains signals generated by/for the feed-back loop for the servo-mechanism (retroaction) that control the confinement: the force command, the fine stage displacement and the estimated pressure at the bottom of the cell.

2. 2 . 4

 24 Measure of the phase velocity: unfolding of the phaseThe computation of the phase velocity is made by first unfold the phase by the proper number of 2π shifts and by assuming a continuity of δφ with respect to the confining pressure. In the first place, by plotting δφ vs P top we use a phenomenological fit formula to shift by 2nπ (n being an relative integer) and collapse all the values of phi (one for each wave packet) on a main curve (more precisely in a range of ±π around this curve).

Figure 2 . 14 :

 214 Figure 2.14: Phase shift δφ = H c φ vs confining pressure P top .Blue circles: raw values of δφ ; black curve is the phenomenological fit (δφ = 4.44 tan -1 (1.5 -0.02 top) + 16.6); red dots: δφ after the unfolding process. Dashed lines define the limit around the phenomenological fit ±π.

Figure 2 . 15 :

 215 Figure 2.15: Twin wave packets. For this wave packet, χ 2 = 7.47 10 -5 whereas χ 2 extended = 0.29.

Figure 2 . 16 :

 216 Figure 2.16: Group velocity vs confining pressure at the top of the cell. Marker color characterizes the value of χ 2 . (a) χ 2 is computed only over the wave packet ; (b) χ 2 is computed over an extended time range. Red arrow points at the wave packet we studied figure 2.15.

Figure 2 . 17 :

 217 Figure 2.17: Group velocity vs confining pressure P . Marker color according to χ 2 which characterizes the success of the fit.

Figure 2 . 18 :

 218 Figure 2.18: Group velocity vs. confining pressure P top , for weakly distorted wave packets.

Figure 2 . 19 :

 219 Figure 2.19: Linear fit of V g vs A sr c in order to measure the velocity at acoustic limit (A sr c → 0). Example for P top = 180 ± 5P a.

2 )Figure 2 . 20 :Figure 2 . 21 :

 2220221 Figure2.20: Group velocity c φ vs P top at acoustic limit. Small squares are the wave packets retained from the previous data selection. Colors reflects the amplitude of the source. Blue squares are values at the acoustic limit (A sr c → 0). Error bars are of the order given by the size of the markers.

Figure 2 . 22 :

 222 Figure 2.22: Group and phase velocity with respect to confining pressure P .

Figure 2 . 23 :

 223 Figure 2.23: Soft Shell contact model: grain in contact with a plane surface. Both the shell and the bulk of the grain deform under the compression.

2. 3 . 1 2 b 2 1722. 3 .e 2 η 2 R

 312232 Building a mean field theory for granular contacts under very low confinementWe consider a grain in contact with an infinite plane by a normal force F N . It has a soft superficial layer of thickness e (see figure3.1 page 181). The Hertz analysis gives the contact area A = πRδ where δ is the overlap between particles and R the radius. The bulk free energy takes dimensionally the form:F b ∼ E(Rδ) 3/2 δRδ while the superficial energy reads:F s ∼ ηERδe (δδ b ) 2 e GRANULAR CONTACTS IN WEAKLY CONFINED GRANULAR MEDIA η 1is the elasticity contrast. Minimizing the free energy with respect to δ b , we get:and get the total free energy for one contact:

P 3 .

 3 GRANULAR CONTACTS IN WEAKLY CONFINED GRANULAR MEDIAand the rescaled pressure:

2 ( 2

 22

2 / 3 P 1

 231

  We estimate φ = 0.6, µ g = 30 GPa, ν g = 0.2, Z = 5.G = 2 √ 2Zφµ g 15π(1 -ν g ) = 11.25 GPaAn excellent approximation can be obtained from a phenomenological cross-over between these asymptotics as:

. 51 )

 51 If we fit by the relation c = A P 1/4 we find A Domenico = 20.734 and and on our dataA CNES =38.254. The ratio A Domenico /A CNES = 1.844 which is close to the ratio of the radius of beads radius in both sets of data Note P 1/4 which is consistent with all the experimental results (see figure3.3).

1762. 3 .Figure 2 . 24 : 1 3 . 2 3 . 3 3 . 4

 32241323334 Figure 2.24: Compression wave velocity with respect to confining pressure. Measurements on glass beads of diameter 45µm for experiments from Makse and Domenico and 1.15mm (±0.15mm) in our experiment.

Figure 3 . 1 :with δ c = e 2 η 2 R

 312 Figure 3.1: Soft Shell contact model: a grain of radius R and soft superficial layer of thickness e is pushed in contact with a plane surface by a normal force F N . Both the shell and the bulk of the grain deform under the compression. δ is the total deformation ; δ b is the deformation of the bulk.

P 3 ( 3 . 11 )

 3311 We therefore get a universal relation between the rescaled celerity:

  large δ/δ c , the expressions simplify into:

Figure 3 . 2

 32 Figure 3.2 reproduces experimental and numerical measures of bulk and shear moduli from Makse and Domenico (data extracted from [48]). All experimental data

10

 10

Figure 3 . 2 :

 32 Figure 3.2: Elastic moduli vs confining pressure. a) Bulk modulus K and b) shear modulus G. Experimental data from Makse (filled circles ; experiments on glass beads of diameter 45µm) and Domenico (opened squares ; experiments on sand, which average grain diameter is 81µm) ; numerical simulations by Makse (black linked circles). Extracted from [48].

Figure 3 . 4

 34 shows an AFM image of the surface of the glass beads we used in our experiment. The typical height of the asperities, which corresponds to the thickness e of the soft superficial layer in our theoretical model, is of order ∼ 100nm.

Figure 3 . 3 :

 33 Figure 3.3: Celerity of compressive waves (group velocity) with respect to confining pressure. Cyan circles are our experimental data ; black circles are computed data points obtained from simulations by Makse [48]. Blue dashed line is the fit of the Hertzian part of the model (we get P * = 41GP a) ; red continuous line is the complete fit of our theoretical model, after locking the value of the parameter P * to the previously fitted value. We get δ c /R = 3.10 -5 .

Figure 3 . 4 :

 34 Figure 3.4: AFM measurements of the surface of the beads. The typical height of the asperities is of order ∼ 100nm. Colors are function of the height.

Figure 4 . 1

 41 Figure 4.1 shows the schematic of the new system. The vibrating plate is fixed onto the aluminum frame. The latter is actuated by a single PZT (the same as the one used for the cell dedicated to compressive waves).

Figure 4 . 1 :

 41 Figure 4.1: Schematic of the shear cell. a) Schematic layout of the shear cell. b) CAD representation of the system without the vibrating plate. We can see the single piezo-actuator that vibrates horizontally the frame on which the plate is fixed (removed for better visualization).

Figure 4 . 2

 42 Figure 4.2 show a picture of the system without the walls of the cell (without and with the vibrating plate). The walls, made of an ABS chamber confined between aluminum plates is similar to the one used for the previous system, for similar acoustic

Figure 4 . 2 :

 42 Figure 4.2: Base of the acoustic shear cell. Without (a) or wih (b) the vibrating plate.The plate in contact with the granular media is fixed to a metallic frame (figure a) which motion is controlled by a piezo-electric actuator and chrysocal blades in order to create a restoring force to maintain the contact with the actuator that only has a pushing work power. A setting screw (not visible here) enables to properly set the pre-stress on this spring system.

Figure 4 . 3 :

 43 Figure 4.3: Preliminary results on the acoustic shear cell. Shear wave velocity c s with respect to confining pressure. Colors of the markers according to the value of χ 2 which characterizes the success of the fit of the wave packet.

Figure B. 2 :

 2 Figure B.2: Localization of the bottom wall. Numerical count of the number of free grains with respect to altitude y in the cell.

Table 2 .

 2 1: Fit values from the optimized fitting of datasets corresponding to outlet aperture 12 to 20mm.

	.1:

  . Data extracted from[47].

	138	2.1. EXPERIMENTAL SET-UP

  Fit of a typical wave packet. Successive fits of the signals according to the fit strategy previously described in table 2.1. Fits of the source (top graph)

							Source signal
	10						Enveloppe (Hilbert transform)
							Gaussian fit enveloppe
							Dispersive Acc. fit
	0
	-10
	11.173	11.174	11.175	11.176	11.177	11.178	11.179	11.18
							Reception signal
	2						Enveloppe (Hilbert transform)
							Gaussian fit enveloppe
							Dispersive Acc. fit
	0
	-2
	11.173	11.174	11.175	11.176	11.177	11.178	11.179	11.18
	Figure 2.13:

Merci à mes amis, et tout particulièrement à Benjamin pour être un ami si proche et depuis maintenant tant d'années. Merci à toi aussi Théau, pour ton amitié, les sessions de grimpe du samedi matin, les soirées et les discussions de rhéologues de bétons et granulaires qui se sont infiltrées dans tous ces moments de nos vies.Olympe, merci d'être là, depuis si longtemps, et de prouver en permanence que la distance ne compte pas.Je voudrais aussi remercier ma famille, qui m'a soutenu, même si le monde de la recherche scientifique lui était inconnu. Je voudrais remercier ma grand-mère, qui m'a toujours dit de faire "au mieux que je pouvais, et que c'était déjà très bien". Merci à mes grand-parents d'avoir été là. Merci à mes parents qui m'ont toujours donné tout l'amour que l'on peut espérer, et qui m'ont toujours soutenu avec tant de conviction. Merci à mon père de m'avoir poussé à toujours continuer malgré les reliefs de la vie. Merci aussi de m'avoir donné le goût de la science, des techniques, et de l'éclectisme. Merci à ma mère de m'avoir ouvert l'esprit, à l'art et à tant d'autres choses qui font la Vie. Merci à tous les deux de m'avoir rendu curieux et de m'avoir tant transmis de vous. Marion, merci d'être là. Merci pour tout ce que tu es, et ce que tu représentes pour moi.

Paul Ehrenfest (1880-1933)

Nowadays, this classification has been slightly modified in order to take into account cases in which derivatives of the free energy diverge (this is the case at the thermodynamic limit which is defined by a diverging volume and number of particles in the system, and a fixed density).

in the case of non-isotropic fluids, the viscosity can be different according to the direction of shear and can therefore be a tensor.

here we make the assumption that the volume fraction of the granular media, i.e. the ratio between the volume occupied by grains over the whole volume we consider, is about

60%. This is a typical value for granular packings of grains of relative similar size diameter and without external loading.

Actually, one can estimate the minimum particle size above which the system become athermal:

As this tube diameter is larger than any of which we use to control the flow rate at the outlet of the smaller reservoir, the flow rate is imposed by the latter. Thus, the flow rate in the channel is set by the diameter of the outlet tube of the intermediate bin.

The Beverloo scaling links the flow rate Q to en effective diameter equal to the aperture diameter minus a few grain sizes:Q ∝ (D -D m )

5/2

Phantom M340, 12 bits version. Sensor size: 2560x1600 ; full frame at maximum 800 frames per second.

. The "chi square test" (or χ 2 test) is a statistical test that probes the matching between two data series. It is based on finding the shift between the two data series for which the sum of the one to one squared difference is the smallest.

We can take θ ≈ θ c = tan -1 µ c as a good approximation.

Note that the value of a can not be provided by the fit of the asymptotic solution of the non-local rheology we do here. However, as the fit of the non-local rheology provides values of µ c and µ w that are consistent with what provides the local model, we can hope that the estimation of a is good.

 rheology . . . . . . . . . . . . . . . . . . . . . 1203.4 Conclusion on the numerical approach and perspectives .. . 125

According to equation

3.8 we have τ c ≈ 1/ √ k n . From[21] we know that for kn ≈ 10 4 , for slow deformations, we can take dt = τ c /100.

We chose to remove ∼ 50 grains, which corresponds to ≈ 1.5d on the flowing height.

As long as one has calibrated the rheology, i. e. as determined values of µ c and a (and , as the rheology is calibrated as a whole: the calibration lies on the use of the non-local model), the value of I b is predictable.

1.1. ELASTICITY OF JAMMED GRANULAR MATTER

Q provides a rough idea of the number of oscillations during the exponential decay 1/(ω 0 ξ).

The roughness of the grains can be a direct consequence of their fabrication process and also on their history (erosion).

Friction coefficient was set by the authors to a "large value" in order to avoid sliding at contact and therefore only probe infinitesimal strain perturbations, i.e. the linear elastic regime.

3.4. CONCLUSION THE EXPERIMENT

Acknowledgments

Je voudrais tout d'abord remercier chaleureusement les membres du jury. Re-

position for the confining pressure we wish to impose for the next parabola.

Feed-back loop

In this experiment, we probe the mechanical response of a granular packing at very low confinement pressures, down to 10P a. We therefore need to be very precise on the control of the pressure on (and within) the grains. The main cause of fluctuations of pressure within the medium is the remnant g-jitter caused by the human imprecision of the pilot and the turbulence in which the aircraft evolves. The feed-back loop uses the value of the current effective gravity g current measured by the 3-axis low frequency accelerometer installed on the experiment frame, and the force sensor that monitors the pushing force of the linear stage on the cell lid, and therefore on the sample. The actuators are the two linear stages previously mentioned.

The behavior of the feed-back loop is described in the bloc diagram figure 2. 3.

The feed-back loop aims at controlling the pressure at the bottom of the box P bot

DATA PROCESSING AND SELECTION

Which is the solution for the propagation of a plane wave in a semi infinite guide (no lid on the cell).

The second boundary condition characterizes the zero-displacement at the lid (at z = H), at any time t. We are therefore looking for solutions of form:

This illustrates the fact that the wall (at z = H) acts like a mirror. Close to the wall, the structure of the solution is finally:

which gives the vertical stress:

The vertical stress on the wall is: .14) where in our case

of t m , the width σ and the phase φ at the origin.

Then we fit the reception. Keep in mind that the phase φ does not change as the wave travels through the medium, so as the carrier wave pulse ω. These quantities are therefore constrained in the fit of the reception: An excellent approximation can be obtained from a phenomenological cross-over between these asymptotics as:

3 5

Identification of the elastic modulus E using the mean field theory

It is interesting to identify the elastic modulus E using the mean field calculation.

From the normal force between grains of shear modulus µ g and Poisson ratio ν g :

one derives the mean field bulk modulus:

and the mean field shear modulus:

where would be 0 in the no-friction case and is 1 in the infinite friction limit, which is the best approximation here.

From Domenico [46], we have the expression of the celerity for the compressional waves in terms of elastic moduli. Thus, using also relations 3.10 and 3.18: The two acoustic studies we presented are complementary: on the one had, the study of compressive waves provided new insights on the nature of contacts, while on the other hand, the shear wave study aims at providing a direct measurement of the shear modulus 1 . This experiment will also provide a better knowledge of the mechanical response of the soft granular contacts to shear excitations. 1 The celerity for shear wave is given by Domenico [46]:

In this thesis we presented several studies from both sides of the jamming transition. In the first part, we investigate the rheology of dense granular flows, especially in the framework of non-locality which we detailed in the review, first chapter. We have then presented an experiment of flow in a narrow channel, where non-local effects are exhibited. In this situation, the stress configuration is not homogeneous, neither along the depth, nor along the width of the channel. However, we chose to probe the non-local rheology model assuming a plug flow, i.e. neglecting the three dimensionality of the system. We showed that the non-local rheology is capable of describing the tails of the profiles, in the so called "creep regime", as well as the flowing region. These supposedly two regimes are in fact altogether described by a single rheological model: there is only one flowing regime. It exhibits non-local effects, which the non-local model we presented describes.

In a second chapter, we present a numerical study which aims to measure the boundary condition at the free surface of a dense granular flow. In order to validate our code, we first recover results on the measurement of relaxation lengths in a custom numerical set-up of a plane shear cell. We then focus on the study of the incline plane. We use the linearization of the non-local rheology around a base state in order to measure the parameters of the rheology. This linearization does not seem valid on the whole height of grains, as seen by the discrepancy of the measurements of the relaxation lengths previously obtained on the shear cell. However, we manage to provide estimates of the boundary condition at the free surface using this linear solution which appears to be a rather good phenomenological fit function of the inertial profile in the incline plane. The value of the inertial number at the free surface seems to be selected by the distance to the yield Y -1 following an affine relation.

In a second part of this thesis, we investigate the elastic properties of granular packing under very low confining pressure. In a first chapter, we presented the state of the art regarding the acoustic anomaly and the scaling expected for the elastic moduli, with respect to the average number of contacts in the packing, and the confining pressure. We then presented the experimental set-up we designed and used on board of the Airbus ZeroG during parabolic flight, which allowed us to reach vanishing pressures in our granular sample. We then detailed the post-treatment of the data and provided measurements of compression wave velocities with respect to the confinement. We then proposed a theoretical framework in order to explain the scaling of sound celerity with the pressure. In a last section, we presented preliminary results of shear wave propagation.

Appendices

Appendix A

Physical principles of parabolic flights

As Isaac Newton published in 1687 in his Principia, the gravitational interaction is unique and universal. As long as the studied body has a mass and is distant from another massive body, both attract each other with the famous law:

where F 12 is the gravitational force exerted by body 1 over body 2, G is the gravitational constant (G = 6.674 10 -11 N.m 2 /kg), m 1 and m 2 are respectively the mass of body 1 and body 2 separated by the distance d ; u 12 is the unit vector pointing from body 1 to body 2. The minus sign in the formula shows the two bodies are attracted.

Therefore it is impossible to avoid gravity. Even in space, any body feels the distant force of gravity from the other bodies. This is exactly the case for a satellite orbiting around Earth: it is only subject to gravity and thus constantly falls down on our planet. Because it has initially been launched with an no-zero velocity tangent to the direction of gravity, its trajectory is not a straight line pointing at the center of Earth but rather an ellipsoid: it constantly missed the ground and keep rotating around Earth . It is a constant free fall, with no friction.

As the state of weightlessness of a massive body is to be only submitted to its weight (and no contact force), the key idea of parabolic flights is to launch the body Again, the apparent gravity is about 1.8g.

Appendix B

Real position of the bottom in the incline plane simulations

In the numerical code, the bottom of the cell is characterized by y = 0 at a physically arbitrary position. However, in order to create a rough wall, we numerically glued the layers of grains at the bottom of the cell (in black on figure 3.13). The first step is to measure the real position of the physical bottom, in order to compute, in a second step, the height of flowing grains which is an important parameter in the incline plane setup.