Quiero Agradecer

Stéphane Morel

Moisés Hinojosa

Igualmente Agradezco

Gilles Pijaudier-Cabot

Daniel Bonamy

Sylvie Gentier

Antonio García Loera

Daniel Billaux Sus

Rogelio Garza

Dr Michel Dumon

Gracias A Hélène

Margaux Michel

Thomas Molina

Cristophe, Huyen, Phan, Myriam Frank

Jean-Marc Amadu

Corine Blain

Pierrette Wiss

Antoine Marache

Joëlle Riss

Shear stress Dilatancy Normal stress Shear stress Residual Peak Ør'

Keywords: joints rocheux, rugosité, fermeture, cisaillement, modèle élément discret (DEM). Mechanical behavior of rock joints: Influence of joint roughness on its closure and shear behavior rock joint, roughness, closure behavior, shear behavior, discrete element model (DEM). Primeramente

Comportement mécanique de joint rocheux : Influence de leur rugosité dans le comportement de fermeture et cisaillement

Résumé

Le comportement mécanique en cisaillement sous contrainte normale constante de joints rocheux est étudié en utilisant une approche numérique par éléments discrets (DEM Discrete Element Model). Les influences respectives de la rugosité des surfaces des joints, de l'élasticité des épontes, de la rupture des aspérités de surface et du niveau de contrainte de compression sur les comportements en fermeture et cisaillement des joints rocheux sont particulièrement analysées. Pour la première fois la rugosité des joints considérée comme auto-affine est utilisée avec DEM pour étudier le frottement des joints rocheux. Cette rugosité est décrite par l'intermédiaire de trois paramètres : exposant de rugosité auto-affine, longueur de corrélation auto-affine et variance des fluctuations de hauteur. Sur la base d'un algorithme fondé sur la méthode spectrale, huit surfaces auto-affines isotropes correspondant à différentes rugosités ont été générées. Ces surfaces numériques sont utilisées comme moules permettant de générer les surfaces composées d'éléments discrets utilisées dans la suite de l'étude. La modélisation par éléments discrets s'appuie sur une calibration des propriétés élastiques effectuée à partir d'un volume élémentaire représentatif suivie de l'implémentation d'un critère elliptique de contraintes de rupture (au niveau des lois d'union entre éléments) permettant de simuler les grandes lignes du comportement quasi-fragile d'un mortier (utilisé lors d'expérimentations antérieures). Sur cette base et une fois les surfaces rugueuses implémentées dans les modèles DEM, les essais de fermeture (test de compression) des huit joints sont effectués sous deux niveaux de contrainte de compression : 14 MPa et 21 MPa. Par la suite, les joints sont cisaillés selon deux directions perpendiculaires. Pour chaque direction de cisaillement et chaque niveau de contrainte de compression, les joints sont testés en utilisant trois modèles mécaniques différents : 1) modèle rigide dans lequel, à l'exception des surfaces de joint en contact, les épontes ne peuvent pas se déformer, 2) modèle élastique dans lequel les épontes peuvent se déformer dans leur volume et 3) modèle élastique-fracture dans lequel les épontes peuvent se déformer dans leur volume et les liens entre les particules peuvent rompre selon le critère elliptique de contrainte. L'utilisation de ces trois modèles mécaniques différents permet d'étudier de façon systématique l'influence de la rugosité seule (modèle rigide), l'influence de l'élasticité et de la rugosité (modèle élastique) et enfin, l'effet combiné de la rugosité, de l'élasticité et de la rupture (modèle élastique-fracture). L'étude des résultats obtenus lors des simulations DEM est accompagnée d'une analyse énergétique permettant d'estimer l'évolution de l'énergie élastique stockée dans le système, de l'énergie de friction, du travail associé à la dilatance du joint et de l'énergie dissipée au cours de l'essai de cisaillement.

Comportamiento mecánico de juntas rocosas: Influencia de la rugosidad en los fenómenos de cierre y cizalladura

Resumen

En esta tesis se estudia la fricción en juntas rocosas utilizando el Método de Elementos Discretos (DEM). En particular, se estudia la influencia de la rugosidad de las superficies de la junta, la elasticidad, la fractura, y el nivel de carga de compresión sobre el comportamiento de cierre y de cizalla de las juntas rocosas. Por primera vez la rugosidad de las juntas considerada como auto-afín es utilizada para estudiar la fricción de juntas rocosas, la rugosidad se describe mediante tres parámetros: el exponente de rugosidad, la longitud de correlación auto-afín y la varianza de alturas. Mediante un algoritmo de computadora basado en métodos espectrales, ocho superficies autoafines isotrópicas con diferente rugosidad fueron creadas. Posteriormente, las ocho superficies fueron utilizadas como moldes para generar las juntas utilizando elementos discretos. Antes de realizar las simulaciones de compresión y cizallaura, se calibraron las propiedades elásticas y de fractura (criterio de fractura elíptico basado en esfuerzos) de las juntas numéricas a los datos experimentales (obtenidos previamente) de unas muestras de mortero mediante la utilización de un volumen elemental representativo (REV). Una vez que las propiedades mecánicas de las juntas se obtuvieron mediante la calibración del REV, se realizaron las pruebas de cierre (prueba de compresión) de las ocho juntas DEM. Se utilizaron dos niveles de esfuerzo de compresión para las pruebas de cierre: 14 MPa y 21 MPa. Después, las ocho juntas DEM fueron cizalladas en dos direcciones mutuamente perpendiculares. Para cada dirección de cizalla y cada nivel de esfuerzo de compresión (14 y 21 MPa), las juntas fueron cizalladas usando uno de los tres modelos mecánicos siguientes: 1) un modelo rígido, en el que las juntas no se pueden deformar, excepto en su superficie, 2) un modelo puramente elástico, en el que las juntas se pueden deformar en todo su volumen y 3) un modelo elástico con fractura en el que las juntas se pueden deformar en su volumen y, si el esfuerzo sobre las uniones entre partículas excede cierto nivel de esfuerzo máximo, las uniones se rompen de una manera irreversible. El uso de estos tres modelos mecánicos nos permitirá estudiar de manera sistemática: la influencia de la rugosidad (modelo rígido), la influencia de la elasticidad y rugosidad (modelo puramente elástico) y, finalmente, el efecto combinado de la rugosidad de las juntas, la elasticidad y la fractura (modelo elástico con fractura). El estudio de los resultados obtenidos de las simulaciones DEM es seguido por una análisis energético el cual permite estudiar la evolución de los diferentes tipos de energía en función del desplazamiento de cizalla: energía elástica almacenada en el sistema, energía de fricción entre elementos discretos, el trabajo relacionado con la dilatación de la junta y la energía disipada por el amortiguamiento interno del DEM.

Palabras clave: juntas rocosas, rugosidad, comportamiento de cierre, cizalla, modelo de elementos discretos (DEM).

Introduction

This thesis proposes to study numerically the influence of the morphology of rock fracture surfaces on their shearing mechanical behaviors under normal stress. The present subject belongs to the rock mechanics discipline which started in the 1950's with the study of the mechanical response of a rock submitted to a disturbance. In general, some of the effects produced by these disturbances can be the translation of the rock mass, rotations, deformation and fracture. In order to understand how the rocks respond to the loads applied to it, is necessary to go deep about the nature of rocks, how they are formed and understand the several forms of loading of a rock mass. Before the study of the mechanical response of rock joints, that is at the heart of this thesis, a general introduction to rock mechanics and generalities of rock joints are going to be presented in this chapter.

Firstly, let us start with a brief resume of the structure and nature of rocks. There are three general categories of rocks according to how they are formed: igneous rocks, metamorphic rocks and sedimentary rocks. Igneous rocks are formed when the hot lava in the asthenosphere rises to the earth surface and cools forming plutonic rocks or volcanic rocks. Igneous rocks were the first formed since the earth formation. As time goes, the solid igneous rocks had been exposed to various chemical or mechanical watering mechanisms. These watering mechanisms decompose the rocks in small parts which are transported and deposited in specific regions (sedimentary basins). The process of watering is continuous and throughout the years, the sediment deposits can reach a thickness such that the material found at the base are exposed to huge pressures and temperatures that causes the particles of sediments to bond together forming rigid sedimentary rocks. In the other hand, metamorphic rocks form when the igneous or sedimentary rock changes its micro-structure and mineralogical composition to form new rocks. The modification of the micro-structure, and hence, its mineralogical and physical properties, can be caused by changes in pressure and temperature at which the rocks had been exposed. The above mentioned mechanisms concerning the rock formation cannot be understood without considering the earth's structure. The modern geophysical theories consider the earth as formed by layers of different composition and properties (Figure 0-1) where the core is the innermost of these layers. The core is a sphere with a radius of about 3470 kilometers and is composed largely of iron and nickel; it is separated in two parts named: inner core and outer core. Because of the high temperature in this region, the outer core is molten although the elevated pressure on it, which can be as much as 1 million times the atmospheric pressure at sea level.

The second layer of the earth is called mantle and lies directly below the crust. It is almost 2900 kilometers thick and makes up 80 percent of the Earth's volume. Although the chemical composition may be similar throughout the mantle, Earth temperature and pressure increase with depth. These changes cause the strength of mantle rock to vary with depth, and thus they create layering within the mantle. The upper part of the mantle consists of two layers:

The Asthenosphere extends from the base of the lithosphere to a depth of about 350 kilometers and is formed by weak, plastic rock. Because the temperature increases with depth, the mechanical properties of the rock in this layer change over a vertical distance of only a few kilometers. The gradual increase of temperature, causes that 1 to 2 percent of the asthenosphere stays in liquid state and flows slowly, perhaps at a rate of a few centimeters per year. At the base of the asthenosphere, increasing pressure causes the mantle to become mechanically stronger, and it remains so all the way down to the core.

The Lithosphere is relatively cool and hence is formed by solid rock whose mechanical behavior is similar to that of the crust. The lithosphere is considered the outer part of the earth and includes both the uppermost mantle and the crust. In most regions, the lithosphere varies from about 75 kilometers thick beneath ocean basins to about 125 kilometers under the continents. The lithosphere is less dense that the asthenosphere and consequently, floats on the last one. The lithosphere is broken into seven large tectonic plates (Figure 0-2) and several smaller ones. The plates move slowly, at rates ranging from less than 1 to about 16 centimeters per year. Because the plates move in different directions, they bump and grind against their neighbors at plate boundaries. The great forces generated at a plate boundary when they collide build mountain ranges and cause volcanic eruptions and earthquakes. These processes and events are called tectonic activity. In contrast to plate boundaries, the interior portion of a plate is usually tectonically quiet because it is far from the zones where two plates interact.

The crust is the outermost and thinnest layer. Because it is relatively cool, it consists of hard, strong rock. Oceanic crust is 5 to 10 kilometers thick and is composed mostly of basalt family rocks. In contrast, the average thickness of continental crust is about 20 to 40 kilometers, although under mountain ranges it can be as much as 70 kilometers thick. Continents are composed primarily of a granite rocks family. It is in this layer of the earth that the human activity develops: construction of dams, buildings, tunnels, mines (mineral deposits), nuclear waste deposits, gas and petroleum extraction, hydrological reservoirs and hot water to thermal energy exploitation, to name just some examples. Thus, the focus of the present thesis and the results derived from it, are expected to apply to the solid rocks forming the earth crust.

Rock mechanics and discontinuities: the influence of geological factors on rocks and rock masses

In the past decades, rocks mechanics has been broadly applied to civil and mining engineering projects. Such projects, as mentioned before, are developed in the earth crust, which can be considered as a stressed medium by natural process, i.e., tectonism, gravity, etc. The earth crust contains different kind of features like joints, faults, folds, etc., that make it non continuum. These natural discontinuities make rock mechanics a special subject and differentiate it from conventional engineering in which the working material (as steel, concrete or soil which is considered as a continuum material in soil mechanics) have not significant fractures/discontinuities and is not prestressed. Further, when surface excavations are made in rock, it is helpful to be able to evaluate the stability of the excavation. This highlights another crucial aspect which has only really been developed since the 1970s, and that is, understanding the full role of the rock structure, i.e. not only the intact rock but also the rock fractures and their three-dimensional configuration. In general, the stability of nearsurface excavations is governed by the rock structure, whereas deeper excavations can be more affected by the intact rock and pre-existing stresses. Thus, the rock structure is particularly important in civil engineering and open-pit mines and so it is necessary to be able to characterize and understand the mechanics of these discontinuities.

As mentioned before, the fractures in the rock generally govern the stability of near surface structures and the natural in situ stresses govern mainly the stability of deep structures. For example, the stability of a dam foundation will depend critically on the deformability and permeability of the underlying rocks, which in turn are dictated by the nature and geometrical configuration of the fractures in the rock mass. This is also true for the stability of the side slopes of surface excavations and the roof and sides of near surface underground excavations. However, at medium depths in weak rocks (for example the channel tunnel between England and France) and at considerable depths in strong rocks (for example South African gold mines) the natural stress, which is altered by the engineering, can be the dominant problem. Furthermore, stability can be influenced by other factors; e.g. whether the rock is wet or dry, cold or hot, rigid or soft. Typical circumstances where these factors are important are the degradation of chalk and mudstones on either exposure to water movement or desiccation. Normally, a rock mechanics project involves several variables: the kind of rock and its mechanical properties (intact rock), the discontinuities forming the rock structure, in-situ stress, boundary conditions (loads applied, stiffness of the surrounding media, i.e.), pore water pressure, and, in some instances, it is also important to consider the rheological behavior of the rock mass, that means, the deformation of the rock structure when the loads applied persisted during several years and giving as a result creep phenomena of the rock structure. In the following sections, a brief discussion of each of the above mentioned variables is going to be presented.

Intact rock

Intact rock is defined in engineering terms as rock containing no significant fractures. However, on the small scale it is composed of grains that form the microstructure, such microstructure is governed mainly by the basic rock forming processes. Subsequent geological events may affect its mechanical properties and its susceptibility to water penetration and weathering effects. The most useful single description of the mechanical behavior is the complete stress-strain curve in uni-axial compression. There are several features of interest obtained from this test: the elastic modulus, the Poisson ration  and the compressive strength (Figure 0-3a). Another feature is the steepness of the descending portion of the curve which is a measure of the brittleness, as illustrated in Figure 0-3b. The form of the complete stress-strain curve is dictated by the nature of the microstructure. For example, basalts with fine-high strength grains and low porosity exhibits a very brittle behavior and has a higher stiffness than a highly porous limestone rock with medium grain size ,furthermore, it presents a more gentle descending part of the curve caused by the gradual deterioration of the microstructure as it is progressively damaged. The intact rock can also have such characteristics as inhomogeneity and anisotropy.

In situ pre-existing rock stress state

May already be an in situ preexisting state of stress in the rock due to different causes. The first of these causes is gravity, which induces a vertical stress component caused by the weight of the rock of the overlying strata, this effect increases with depth. In the other hand, a high horizontal stress is necessarily to satisfy the stability conditions, moreover, horizontal stresses can be induced by tectonic activity. There are also factors such as surface topography and erosion which affect the vertical and horizontal stress state. In the other hand, there are cases when the load is applied as the project develops, i.e., a dam or nuclear power station foundation. In the excavation of a tunnel or mine stope, no new loads are applied in unsupported excavations: it is the pre-existing stresses that are redistributed by the engineering activity with decompression problems. In the examples given above, the final result is redistribution of the stresses as the project develops.

It is very important in every rock mechanics project to be aware of the stress state that can be present. In particular it should be remarked that, in the majority of the stress states measured throughout the world, one horizontal component of the stress field has greater magnitude than the vertical component. Local concentrations of stresses (in magnitude and direction) must also be taken into account in projects.

Pore fluids and water flow

Groundwater comprises almost 60% of the clean water in the continents. This water travels through the soil and rocks by the effect of gravity and the flux of water depends on the permeability and hydraulic conductivity of the soil and rocks through which it travels. It turns out that the permeability of rock masses is mainly dictated by the discontinuities that form the rock structure, when water travels through the discontinuous rocks, it can be under pressure and hence reduce the effect of the applied stresses. This leads to the concept of effective stresses which have proved to be so important in soil and rock mechanics, both from the theoretical and applied points of view.

Both the in-situ stress and the water flow are significantly affected by engineering activity, i.e. all excavation affects the stress state and all excavations act as sinks because the in situ hydraulic pressure is reduced to atmospheric pressure on the excavation boundary. There are some applications in which an understanding of fluid flow through intact rock can be critical, e.g. in reservoir engineering for the petroleum industry, in hot dry rock geothermal energy projects where the success depends on achieving the required borehole-to-borehole water flow, the case of radioactive waste disposal where the engineer can only state that the design of the waste repository is valid if the radionuclide dosage back to the biosphere can be calculated, and the storage of CO2, which is very common recently. This can only be done if the three-dimensional flow of water through fractured rock masses can be accurately modeled.

Influence of time

Another major factor of importance in rock mechanics is the influence of time. Nature have taken millions of years to create massive rocks on a meta-stable equilibrium, while engineering projects develops in the order or years or decades. This time contrast must be taken into account in rock mechanics. One of the basic tools in mechanics is the theory of elasticity which links stresses and strains by the instantaneous response of the rock. Because there is no time component in elasticity, this theory is unlikely to fully explain geological processes. The theory is, however, likely to be of considerable assistance in engineering when we are interested in the initial redistribution of the stress field upon excavation.

The influence of time is important because of factors such as the decrease in rock strength through time, and the effects of creep and relaxation. Creep is increasing strain at constant stress; relaxation is decreasing stress at constant strain. We might be considering processes which occur very rapidly, in particular, stress waves travelling through the rock. These could be caused by natural processes, as in earthquakes, or by artificial processes such as blasting or mechanical excavation using picks, discs or button cutters. Hence, throughout the time range from milliseconds to millions of years, the engineer should have some understanding of rate effects.

Discontinuities and rock structure

One of the principal factors that differentiate the subject of rock mechanics from other branches of mechanics of materials is the presence of significant discontinuities in the rocks. Rocks can be fractured when a sufficiently high stress is applied to them through natural or artificial processes.

Commonly, the rock may have already failed and formed faults and joints, these faults and joints may be the 'weak links' in the rock structure. During the lithification process and throughout geological history, there have been orogenic periods and other less severe loading processes applied to the rock. The result in terms of the rock fracturing is to produce a geometrical structure (often very complex) of fractures forming rock blocks (Figure 0-4). Because in the general uses of mechanics and stress analysis it is assumed that a material is continuous, the geological features such as: faults, joints, bedding planes, etc., are termed 'discontinuities'. These discontinuities have many geometrical and mechanical features which often govern the whole behavior of the rock mass (especially at shallow deeps in the earth crust). The discontinuities may have certain shapes, sizes and can be orientated in certain directions. The overall geometrical configuration of the discontinuities in the rock mass is termed rock structure. It is fundamental to understand this geometrical structure when developing engineering projects in/on rock.

From a fracture mechanics point of view, there are three ways in which a fracture can be formed: one by pulling apart and two by shearing (Figure 0-5), these modes are called Mode I, Mode II and Mode III. The mode of fracture can give a rough idea of the likely mechanical characteristics of the discontinuities. Given that such discontinuities exist in all rock masses at a variety of scales, it is not surprising that they significantly affect the deformability, strength and failure of rock masses. Moreover, other key characteristics such as the permeability, can be governed almost entirely by the rock structure configuration. Hence, part of this introductory chapter is dedicated to present the principal kinds of discontinuities, and, at the end of this chapter, focus is going to be make on one kind of discontinuity called joints. The mechanical behavior of such rock joints is at the heart of this thesis and detailed study of joints is going to be presented in Chapter 1.

Major types of structural features

Structural features and their origins are well described in several textbooks on general, structural and engineering geology. All that will be given here is a catalogue of the major types of structural feature and brief description of their key engineering properties.

Bedding planes

Bedding planes represent interruptions in the course of deposition of the rock mass Figure 0-6, they are generally highly persistent features, although sediments laid down rapidly from heavily laden wind or water currents. Bedding planes may contain material of different grain size from the sediments forming the rock mass, or may have been partially healed by low-order metamorphism. In either of these two cases, there would be some 'cohesion' between the beds; otherwise, shear resistance on bedding planes would be purely frictional. Arising from the depositional process, there may be a preferred orientation of particles in the rock, giving rise to planes of weakness parallel to the bedding.

Folds

Folds are structures in which the attributes of the beds are changed by flexure, i.e., resulting from the application of post-depositional tectonic forces (Figure 01234567). They may be major structures on the scale of a mine or mining district or they may be on a smaller local scale. Folds are classified according to their geometry and method of formation. The major effects of folds are that they alter the orientations of beds locally, and that certain other structural features are associated with them.

In particular, well defined sets of joints may be formed in the crest or trough and in the limbs of a fold. During the folding of sedimentary rocks, shear stresses are set up between the beds where slip may occur. Consequently, the bedding plane shear strength may approach, or be reduced to, the residual shear strength. Axial-plane or fracture cleavage may also develop as a series of closely spaced parallel fractures resulting from the shear stresses associated with folding.

Faults

Faults are defined as fractures on which identifiable shear displacement has taken place (Figure 0-8).

They may be recognized by the relative displacement of the rock on opposite sides of the fault plane. The sense of this displacement is often used to classify faults. Faults may be pervasive features which traverse a project's area or they may be of relatively limited local extent on the scale of meters; they often occur in echelon or in groups. Fault thickness may vary from meters in the case of major, regional structures to millimetres in the case of local faults. This fault thickness may contain weak materials such as fault gouge (clay), fault breccia (recemented), rock flour or angular fragments. The wall rock is frequently slickenside and may be coated with minerals such as graphite and chlorite which have low frictional strengths. The ground adjacent to the fault may be disturbed and weakened by associated structures such as drag folds or secondary faulting. These factors result in faults being zones of low shear strength on which slip had readily occurred.

Shear zones

Are bands of material, up to several meters thick, in which local shear failure of the rock has previously taken place (Figure 0-9). They represent zones of stress relief in an unaltered rock mass.

Fractured surfaces in the shear zone may be slickenside or coated with low-friction materials, produced by the stress relief process or weathering. Like faults, shear zones have low shear strengths but they may be much more difficult to identify visually.

Dykes

Dykes are long, narrow intrusions of generally fine-grained igneous rock with steep or vertical and approximately parallel sides (Figure 0-10). They may vary in width from a few centimeters to several meters and may appear as dyke swarms. Some dyke rocks are more resistant to weathering than the country rock, but the basic igneous dyke rocks such as dolerite can weather to montmorillonite clays which are noted for their swelling characteristics. The dyke margins are frequently fractured and altered during the intrusion. They form potential seepage paths and zones of low stiffness and shear strength in which movements will tend to be concentrated. Because of their high stiffnesses, unweathered dyke rocks can develop high stresses and so be susceptible to stress-induced failure or rock burst, as in the deep-level gold mines of South Africa.

Joints

Joints are the most common and generally the most geotechnically significant structural features in rocks. Joints are fractures of geological origin (Figure 0-11), and it is common that such discontinuities do not occur at completely random orientations: they occur with some degree of 'clustering' around preferred orientations associated with the formation mechanisms. Hence, it is sometimes convenient to consider the concept of a discontinuity set (which consists of parallel or sub-parallel discontinuities), and the number of such sets that characterize a particular rock mass geometry a group of parallel joints is called a joint set. Joints may be open, filled or healed. They frequently form parallel to bedding planes, foliations or slaty cleavage, when they may be termed bedding joints, foliation joints or cleavage joints respectively. Sedimentary rocks often contain two sets of joints approximately orthogonal to each other and to the bedding planes. These joints sometimes end at bedding planes, but others, called master joints, may cross several bedding planes. Veins, or cemented joints, are mineral infillings of joints or fissures. They may be sheet-like, tabular or irregular. They are generally of igneous origin but may also result from sedimentary processes. They are commonly associated with metalliferous orebodies and have been found to have major influences on orebody cavability and fragmentation. Veins may be weaker or stronger than the wall rock and should be taken into account in rock mass studies. Joints have a significant influence in the mechanical and hydraulic behavior of rock masses. However, nowadays it is not completely understood how the roughness, mechanical properties and boundary conditions of the joints influence their frictional behavior; one of the principal signs to arrive to such argument, is that there exists numerous empirical and theorical models which try to explain and predict the joint friction, [START_REF] Bandis | Experimental studies of scale effects on the shear behavior of rock joints[END_REF], [START_REF] Swam | Prediction of shear behavior of joints using profiles[END_REF], Maksinovic (1992), [START_REF] Leong | A model for rock interfacial behavior[END_REF], [START_REF] Haberfield | The application of energy principles to the determination of the sliding resistance of rock joints[END_REF], [START_REF] Maksimovic | The shear strength components of a rough rock joint[END_REF]. However, no one can give a clear explanation how the friction resistance in rock develops and which are the main variables involved.

The occurrence of joints

As mentioned before, all rock masses are fractured, and it is a very rare case where the spacing between discontinuities are appreciably greater than the dimensions of the rock engineering project. Very often major discontinuities delineate blocks within the rock mass, and, within these blocks, lead to a further suite of discontinuities. Such hierarchical systems of discontinuities may be much more complex, but this will not affect the idea that in general it is expected a relation of the form:

~1 ~1 0-1
Some of the most important aspects of discontinuities are their occurrence (mean value and distribution of spacing between discontinuities), orientation, persistence, discontinuity sets, aperture and roughness. Below it is given a brief description of each of these concepts.

Spacing and frequency: spacing is the distance x between adjacent discontinuity intersections with a measuring scanline. Frequency λ (i.e. the number per unit distance) is the reciprocal of spacing (i.e. the mean of these intersection distances). Figure 0-12 illustrates a sampling line through a rock mass, which intersects a number of discontinuities. The length of the sampling line is , the number of discontinuitites intersected is , and hence, the discontinuity frequency is:

= 0-2
and the average spacing is:

̅ = 0-3
It is also possible to consider the distribution of the individual spacings between fractures. When a sufficiently large sample of these individual spacing values is plotted in histogram form, a negative exponential distribution is often evident. It should be noted that the general trend of this histogram is to be many small spacing values and few very large spacing values in the distribution, which can be expressed as the probability density function: ( ) = . Note that the mean of this distribution is 1/ , and the standard deviation is also 1/ . It is a one-parameter distribution, with the mean and the standard deviation being equal. This negative exponential distribution is the spacing distribution associated with the Poisson process of random events. The negative exponential distribution is analogous to the normal distribution, except that it is the distribution to which the spacing values converge when successive spacing distributions of any type are superimposed on the sampling line. In other words, the occurrence of the negative exponential distribution is expected as the result of a suite of superimposed geological events, each of which produces fracturing of a given distribution. It should be remarked that the measured discontinuity frequency is expected to vary with the direction of the sampling line relative to the orientation of the discontinuities.

Orientation, azimut/dip angle: If we assume that a discontinuity is a planar feature then, the azimut (the compass bearing of the steepest line in the plane) and the dip angle (the angle that this steepest line makes to the horizontal plane) uniquely define the orientation of the discontinuity.

Persistence, size and shape: the persistence is defined as the extent of the discontinuity in its own plane, incorporating factors such as the shape of the bounded plane and the associated characteristic dimensions.

Block size:

In terms of excavation and support, it is helpful to have an estimate both of the average block size and the block size distribution, which is (in situ) analogous to of the particle size distribution used in soil mechanics. The block size is directly related to the spacing distribution of joints, their persistence and intersection between joint sets.

L x

Roughness: although discontinuities are assumed to be planar for the purposes of orientation and persistence analysis, the surface of the discontinuity itself is rough. Discontinuity roughness may be defined either by reference to standard charts or mathematically. The word 'roughness' is used to denote the deviation of a discontinuity surface from perfect planarity, which can rapidly become a complex mathematical procedure utilizing three-dimensional surface characterization techniques, whether these be by polynomials, Fourier series, noise waveforms, fractals or geostatistics. Roughness is related to various mechanical and hydraulic properties of discontinuities. On the purely geometrical side, it is possible to predict the amplitude of asperities from the roughness and profile length of the joints. On the mechanical side, shear strength is intimately related to the joint roughness. Moreover, there are obvious implications for aperture and variation in aperture as a function of discontinuity roughness.

Aperture: is defined as the perpendicular distance between the adjacent rock surfaces (wall joints) of the discontinuity. This will be a constant value for parallel and planar adjacent surfaces, a linearly varying value for non-parallel but planar adjacent surfaces, and completely variable for rough adjacent surfaces. This parameter has mechanical and hydraulic importance, and a distribution of apertures for any given discontinuity and for different discontinuities within the same rock mass is to be expected. Current research in the hydraulic context indicates that a discontinuity cannot be approximated as two parallel planes because of the phenomenon of channel flow, where the fluid mainly flows through certain channels within the discontinuity created by tracks of larger local apertures. We will see in future chapters that the formation of such channels depends of the relative displacements between both sides of the discontinuities.

As mentioned before, the present thesis studies the mechanical behavior of rock joints, in particular, the main objective is to develop a methodology which permits to study numerically (Discrete Element simulations) the influence of the roughness, elasticity and fracture of rock joints on their compression and shear behaviors. After this introduction to the importance of rock joints in the mechanical and hydraulical behavior of rock masses, Chapter 1 will resume some of the theories for the closure and shear behavior of rock joints, then, Chapter 2 will be dedicated to the study of the morphology of rock joints and a mathematical procedure will be proposed allowing the generation of self-affine surfaces that will be used to mimic the actual rough surfaces of rock joints. Then, Chapter 3 will present the Discret Element Method (DEM) that will be used to study the mechanical behavior of joints. After the calibration of elastic and fracture properties on a Representative Element Volume (REV), the approach used to implement rough surfaces on DEM is detailed and the approximation related to the discretization of the rough surfaces is evaluated. Finally, Chapter 4 presents the numerical simulations of the compression and direct shear tests and the results compared to experimental data previously obtained by Flamand (2000) (model validation). Parallel to the model validation, the results obtained from these numerical simulations are analyzed and discussed in light of the differences in behavior obtained from simulations performed for rigid, elastic and elasticfracture material behavior. An energetic analysis of the direct shear test is finally proposed to estimate the contribution of each behavioral components involved in the response (i.e., the energies linked to the friction and the dilatancy of the joint, the elastic energy stored in the system and the dissipated energy linked to the fractures). These are just some examples of the potentialities of the methodology proposed to study numerically the mechanical behavior of rock joints.

Mechanical behavior of rock joints

This chapter is dedicated to discuss the principal mechanical properties of the rock joints, i.e., closure behavior and shear behavior. Firstly, the closure behavior observed experimentally is exposed followed by the presentation of several models. Then, the experimental shear behavior of rock joints is presented including the set of parameters that influence it. Finally, a review of several models for the shear stress is presented.

The load distribution inside a discontinuous rock mass can be very complex, in order to simplify the study of stability of rock joints it is essential to obtain a set of two resultant forces applied to the discontinuity, one perpendicular to the mean plane of the discontinuity (normal load) and the other one acting parallel to the mean plane of the discontinuity (shear load). From this viewpoint, the discontinuity can be considered as submitted to compression loads and shear loads. On this basis, let us analyze the response of a joint when compression and shear loads are applied to it.

Stiffness

When rock joints are submitted to a compression load, joint's surfaces are gradually pushed together, with an obvious limit when the two surfaces are closed (because of the joint roughness, there is not a perfect closure). The stiffness associated with this compression process gradually increases with applied stress or displacement, again reaching a limit associated with the strength of the intact rock, as indicated in Figure 1-1a. In tension, no tensile stress can be sustained because discontinuities are regarded as having no tensile strength by definition. In the other hand, when a discontinuity is subjected to shear stress or shear displacement, the shear load vs. shear displacement curve is rather like the complete stress-strain curve for compression of intact rock, except of course that all failure is localized along the discontinuity. There is an initial shear stiffness, a shear strength and a post-peak failure locus (Figure 1-1b). [START_REF] Barton | Strength, deformation and conductivity coupling of rock joints[END_REF].

However, for the analysis it is convenient to assume either one global linear stiffness value or a composite piecewise linear approximation. As a first approximation we can consider the linear stiffness representations as for the normal case and for the shear case. We can also consider the possibility that a normal stress will cause a shear displacement (effect that is usually negligible), using a constant that a shear stress will cause a normal displacement, using a constant These stiffnesses have the dimensions of, for example, MPa/mm, because they relate stress to displacement. With these linear approximations for the stiffnesses:

a) b) K = + 1-1 = + 1-2
or, in matrix notation:

= 1-3
This final expression also permits evaluation of the displacements when the stresses are known through use of the inverse of the matrix . It should also be noted that this matrix, containing offdiagonal terms, provides a first approximation to the coupling of normal and shear stresses and displacements.

For the displacements of in situ discontinuities, the stiffness of the surrounding rock system will also need to be taken into account. There are many other practical aspects, such as the fact that the shear stiffness may be anisotropic in the plane of the discontinuity due to surface anisotropy itself (existence of striations or lineations, etc.).

Description of the experimental compression test

The deformation test of a rock joint under normal load, also called as joint's closure test, is normally carried out in an uni-axial press. At laboratory scale, samples of cylindrical shape are used (Figure 1-2a). Several displacement sensors are positioned around the sample in order to measure the closure of the joint when a normal load is applied. In order to correctly estimate the joint's closure, it is necessary to substract the intact rock deformation, [START_REF] Gentier | Morphologie et comportement hydromécanique d'une fracture naturelle dans le granite sous contrainte normale -etude expérimentale et théorique[END_REF]. Eq. 1-4 is used to compute the joints closure.

∆ = ∆ -. ℎ 1-4
where:

∆ is the joint's closure (mm). ∆ is the vertical displacement measured by the sensor (mm). σ is the normal load applied (MPa). h is the distance between both points of the displacement sensor (mm). E is the Young's modulus of the intact rock (MPa). [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF] after [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF].

It is important to distinguish two types of joint closure tests a) when the upper and lower surfaces match and b) when they do not match. An unmatched joint can deform more than a matched one. In the present thesis, only matched joints are studied.

Models of joint closure

As mentioned before, the present thesis is focused on the shear behavior of rock joints, however, in the next section, a brief review of several models for the joint closure is presented. For a more detailed discussion of the various closure models for rock joints, it is recommended to consult the thesis of [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF].

Shehata (1971)

Shehata proposes the following equation to estimate the closure of a rock joint:

∆ = ( ) - ( ) 1-5
where and are a reference normal stress and actual normal stress respectively, and is a constant that depends on the elastic behavior of the joint.

Goodman (1974)

Goodman introduces the idea that the joint closure tends to a maximum value at high normal stresses. He proposes the following function to estimate the closure of the rock joints.

= ∆ ∆ -∆ + 1-6
where and are a reference normal stress and the actual normal stress, and ∆ and ∆ are the actual joint closure and the maximum joint closure respectively. Then, [START_REF] Goodman | Methods of geological engineering in discontinuous rocks[END_REF] modifies Eq. 1-6 introducing the mechanical properties of the rock and a new equation to estimate the normal stiffness of the joint.

= ∆ ∆ -∆ + 1-7 = ( -) 1 - ∆ ∆ ∆ 1-8
where A and t are constants reflecting the mechanical properties of the rock, and is the normal stiffness of the joint.

Bandis (1980)

Bandis also consider a maximum closure of the joint and proposes the following equations to estimate the joint normal stiffness and joint closure:

= 1 - ∆ ∆ ∆ 1-9 and = ∆ -∆ 1-10
where a and b are constants that depends on the mechanical properties of the rock.

Bandis (1983)

In an effort to consider the influence of the joint roughness on its closure behavior, Bandis proposed the following relations for the normal stiffness of the joint and its closure as:

∆ = + ( ) + 1-11
Where A, B, C and D are constants and is the maximum opening of the joint. JRC is the Joint Roughness Coefficient and JCS is the Joint Compresive Strength. The normal stiffness of the joint is estimated by:

≈ 0.02 + 2 -10 1-12

Brown and Scholz (1986)

Brown and Scholz propose a logarithmic function of the normal stress to estimate the joint closure. They were the first to consider the influence of roughness in the joint closure:

∆ = + ( ) 1-13
where C and B are constants determined from the roughness of the joint. The normal stiffness is estimated from: = 1-14

Marache (2002)

Marache improve a closure model that was initially developed by [START_REF] Hopkins | The effect of surface roughness on joint stiffness, aperture, and acoustic wave propagation[END_REF][START_REF] Hopkins | The implications of joint deformation in analysing the properties and behavior of fractured rock masses, underground excavations, and faults[END_REF]. Improvements with respect to the initial model include the possibility of modelling greater fracture sizes, and the inclusion of total force boundary conditions. The fractures are modelled as the interface between two elastic half spaces. The surface roughness of the interface walls is modelled by discretizing the surface topography maps into a grid of elements. Each element thus represents an asperity on the elastic half space. In order to correctly evaluate the normal stress levels, asperities are represented by cylinders of radius R, with the height H of each cylinder being determined from the void space map.

In addition, the material's mechanical properties, i.e. Young's modulus (E) and Poisson's ratio (ν), are required. Different material properties can be assigned to each of the two half spaces and to each of the asperity subsets. The model is based on the computation of forces acting between both fracture walls, at each contact zone, for a specified applied normal load. It is then possible to calculate the total displacement at each contact asperity, as the sum of displacements arising from:

 deformation of the asperity itself due to the forces acting directly on it.  deformation of the bulk material surrounding the asperity due to the forces acting on the asperity and,  deformation of the bulk material surrounding the asperity due to the forces acting on all other asperities in contact (deformation of the bulk material results in mechanical interaction between contact zones).

A system of equations is established to describe the total displacement at each contacting asperity, in terms of the normal stress acting at each point. The system is solved by imposing either a far-field normal displacement, or a total force boundary condition, and by searching for the distribution of stresses which minimizes the strain energy in the system.

Xia et al. (2003)

Xia et al. propose a model considering the waviness and unevenness of the joint surfaces. To develop their model, they assume that: a) the asperity peaks of the unevenness component near the contact can be approximated as a sphere which has a root mean square radius and, the asperity peaks of the waviness component at the contact region can also be approximated as a sphere, which also has a curvature radius; b) the distribution of the asperities is random on the surface. c) The asperity deformations of the unevenness and waviness components obey the Hertz theory. d) The deformations are not affected by the adjacent asperities and waviness. e) The correction for tangential stress in oblique contact between asperities is neglected when the contact angle is small. The joint closure deformation is a function of the asperity deformation, the elastic deformation of the waviness and the distance between the center of the waviness peak and the contact boundary on the unevenness surface.

Shear strength of rock joints

Now let us consider the strength of discontinuities submitted to shear loads, but firstly a discussion of the typical shear tests is going to be presented.

Shear testing

In rock mechanics problems other than those involving only fracture of intact rock, the shear behavior of discontinuities can be important. Conditions for slip on major pervasive features such as faults or for the sliding of individual blocks from the boundaries of excavations are governed by the shear strengths that can be developed by the discontinuities concerned. In addition, the shear and normal stiffnesses of discontinuities can exert a controlling influence on the distribution of stresses and displacements within a discontinuous rock mass. The shear strength and shear stiffness can be measured in one single shear test. There exist various types of shear test, in which the main difference lies on the boundary conditions of the joint sample.

The most commonly used method for the shear testing of discontinuities is the direct shear test. As shown in Figure 1-3, the mean plane of discontinuity surface is aligned parallel to the direction of the applied shear force. The two walls of the specimen are fixed inside the shear box using a suitable encapsulating material, generally an epoxy resin or plaster. Direct shear tests in the configuration of Figure 1-3a are usually carried out at constant normal force or constant normal stress. Tests are most frequently carried out on dry specimens, but many shear boxes permit specimens to be submerged and drained and hence the shear tests are carried out with excess joint water pressures. During the test, the tangential force applied is measured by a load-cell and the vertical displacement (also called dilatancy) of the upper sample is registered by a displacement sensor. During the test, the upper and lower samples do not rotate. [START_REF] Jaeger | Friction and sliding of joints[END_REF].

The typical results of the direct shear test are presented in Figure 1-4 in a plot of tangential stress (τ) vs shear displacement (ΔU) and another plot of the vertical displacement (ΔV) of the upper sample vs shear displacement (ΔU), [START_REF] Archambault | The evolution of void spaces (permeability) in relation with rock joint shear behavior[END_REF][START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF].

a) b)

Figure 1-4 Typical Experimental results of a direct shear test [START_REF] Archambault | The evolution of void spaces (permeability) in relation with rock joint shear behavior[END_REF][START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF].

According to [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF], in the plot of τ vs ΔU, the curve A corresponds to a classical result for a rough rock joint, while curve B correspond to a planar joint. C and D show two examples of dilatancy curves. [START_REF] Archambault | The evolution of void spaces (permeability) in relation with rock joint shear behavior[END_REF] suggest to separate the curves τ vs ΔU and ΔV vs ΔU in five phases:

Phase I: the joint deformation is elastic and pseudo-linear. The tangential stiffness of the joint Ks is estimated during this phase. The upper and lower samples forming the joint tend to match perfectly often implying a contractancy of the joint (negative ΔV). The contact area increases progressively until the sliding of the joint starts, leading to a dilatancy behavior.

Phase II: A non linear relation between τ and ΔU is observed. The dilatancy increases implying a lower contact area than in phase 1.

Phase III: This phase corresponds to the behavior observed in the vicinity of the shear peak load which coincides with the higher angle of dilatancy. The degradation of the surfaces starts. For a planar joint (curve B) there is not such peak of the shear load.

Phase IV: The shear load and the dilatancy angle decrease. The degradation of the surface continues and the particles detached by the damage process pulverize.

Phase V: The shear load tends to a constant (residual shear load). This residual load/stress is reached quickly for planar joints compared to rough joints. In this phase, the slide of the surface depends only on the topographic structure of large wavelength. The orientation of such structures with respect to shear displacement direction makes that the dilatancy can be positive, null or negative. The contact area increases due to the contact of the joint walls with the detached particles.

Note that during a shear test, the zones in contact tend to be heated. The temperature of such heated zones depends on the normal load, shear speed and the thermal properties of the rock.

The triaxial cell is sometimes used to investigate the shear behavior of discontinuities. Specimens are prepared from cores containing discontinuities inclined at 25-40• to the specimen axis. A specimen is set up in the triaxial cell as shown in Figure 1-3b and then the cell pressure and the axial load are applied successively. The triaxial cell is well suited to testing discontinuities in the presence of water. Tests may be either drained or undrained, preferably with a known level of joint water pressure being imposed and maintained throughout the test. In an attempt to overcome the need to obtain, prepare and set up several specimens containing similar discontinuities, a stage testing procedure is sometimes used. A specimen is tested at a low confining pressure as outlined above. When appears slip on the discontinuity, loading is stopped, the cell pressure is increased to a new value, and loading is recommenced. By repeating this process several times, a number of points on the peak strength envelope of the discontinuity can be obtained from the one specimen. However, this approach exacerbates the major difficulty involved in using the triaxial test to determine discontinuity shear strengths, namely the progressive change in the geometry of the cell-specimen system that accompanies shear displacement on the discontinuity.

Factors that affect the shear strength of rock joints

Factors such as the degree of matching between the two walls of the joint, type of infilling, joint roughness, joint size, etc., can have an effect on the strength of discontinuities. Therefore, let us study in a general way the kind of outcome that the factors mentioned before can induce on the joint friction. This thesis being focused on the effect of the roughness on the shear response, the following of the chapter is dedicated to the current theories of the roughness effect on shear strength.

Experimental factors

Influence of scale

Another factor that can influence the mechanical behavior of rock joints is the size of the sample. One of the more remarkable studies on size effects of rock joints was made by [START_REF] Bandis | Experimental studies of scale effects on the shear behavior of rock joints[END_REF]. In their experiments, they show the influence of the specimen size on the shear behavior of rock joints, especially on the shear peak strength and on the shear stiffness. As the size of the specimen decreases the shear peak and the shear stiffness increase (Figure 12345). [START_REF] Bandis | Experimental studies of scale effects on the shear behavior of rock joints[END_REF].

There are two main phenomena involved in the size effect on rock friction. The first one deals with the size effect related to the joint roughness. In general terms, as the size of the specimen increases the roughness decreases. The second contribution to size effect of rock joints is linked to the fracture of the sample. There are numerous experimental data that show an inverse relation between the specimen size and its fracture strength.

Influence of boundary conditions

All of the previous phenomena are discussed considering direct shear test experiments at constant normal load. From this kind of test, the shear of the specimen leads to a dilatancy of the rock joint. Although this test may reproduce discontinuity behavior adequately in the case of sliding of an unconstrained block of rock from a slope, it may not be suited to the determination of the stressdisplacement behavior of discontinuities isolating a block that may potentially slide or fall from the periphery of an underground excavation at considerably depths. In the former case, dilation is permitted to occur freely, but in the latter case, dilation may be inhibited by the surrounding rock and the normal stress may increase with shear displacement (Figure 123456). This helps explain why limiting dilation on discontinuities by rock bolt, dowel and cable reinforcement, can stabilize excavations in discontinuous rock.

Influence of matching (roughness differences between joint surfaces)

Rock joints are formed at various stages and times. Through various geological processes, the initial morphology of joints can be distorted from weathering, shearing, loading and thermal cycles, and hence joints could present different surface profiles on each side and different degrees of matching. In general, joints that have been weathered, sheared or are under loading and thermal cycles are poorer in matching compared to fresh tensile fractures. It is not possible to characterize the degree of matching from the roughness alone. However, roughness and matching may have significant influences on the mechanical, hydraulic and thermal properties of the joints, e.g. as the degree of matching of the joint reduces, its shear peak strength approaches its residual shear strength and its hydraulic conductivity increases, Zhao (1997).

Natural factors

Infilling materials

Discontinuities may contain infilling materials such as gouge in faults, silt in bedding planes, lowfriction materials such as chlorite, graphite and serpentine in joints, and stronger materials such as quartz or calcite in veins or healed joints. The presence of these materials can influence the shear behavior of discontinuities. The presence of gouge or clay seems to decrease both stiffness and shear strength. Low-friction materials such as chlorite, graphite and serpentine can markedly decrease friction angles, while vein materials such as quartz can serve to increase shear strengths. Thus, if a discontinuity contains soft and weak infilling materials, this discontinuity can exhibit mechanical properties similar to those observed for clays and sills. The shear strength of these materials is usually described by an effective stress in the form of Coulomb law of friction. From a laboratory study of such filled discontinuities, [START_REF] Ladanyi | Shear strength and deformability of filled indented joints[END_REF] observed that:

 For most filled discontinuities, the peak strength envelope is ranged between the one corresponding to the infilling material and the one obtained for a similar clean discontinuity.  The stiffnesses and shear strength of a filled discontinuity decrease with increasing filling thickness, but always remain higher than those of the filling alone.  The shear stress-displacement curves of filled discontinuities often have two portions: the first reflecting the deformability of the infilling materials before rock to rock contact is active, and the second reflecting the deformability and shear failure of rock asperities in contact.  The shear strength of a filled discontinuity does not always depend on the thickness of the filling. If the discontinuity walls are flat and covered with a low-friction material, the shear surface will be located at the filling-rock contact. Swelling clay is a dangerous filling material because it loses strength on swelling and can develop high swelling pressures if swelling is inhibited.

Influence of surface roughness on shear strength

Shear tests carried out on smooth and clean discontinuity surfaces under constant normal stress generally give shear stress-normal stress curves as shown in Figure 1234567. When a number of such tests are carried out at a range of effective normal stresses, a linear shear strength envelope is obtained. Thus the shear strength of smooth, clean discontinuities can be described by the simple Coulomb law:

′ = ′ tan ′ 1-15
where ′ is the effective angle of friction of the discontinuity surfaces. However, discontinuity surfaces occurring naturally are never as smooth as the artificially prepared surfaces which gave the results exhibited Figure 1234567.

The shear force-shear displacement curve shown in Figure 1-8a is typical of the results obtained for clean and rough discontinuities. Under constant normal stress, the peak of the shear stress is reached after a small shear displacement and the post-peak regime can exhibit a residual strength for sufficiently large shear displacements. Moreover, if the shear tests are performed for various normal stresses, the shear stress peak and the residual strength exhibits usually the behaviors shown in Figure 1-8b. This behavior can be explained in terms of surface roughness using a simple model introduced by [START_REF] Patton | Multiple modes of shear failure in rock[END_REF] as shown in Figure 123456789. A smooth, clean, dry discontinuity surface has a friction angle , so that a limiting equilibrium for the direct shear test configuration gives:

= tan 1-16
If the discontinuity surface is inclined at an angle (Figure 1-9b) to the direction of the shear force, , then slip will occur when the shear and normal forces on the discontinuity, S* and N*, are related by: * = cos -N sin 1-18 and * = cos + N sin

1-19

Substitution of these values in Eq. 1-16 and simplification gives the condition for slip as:

= tan( + ) 1-20
Thus inclined discontinuity surface has an apparent friction angle of ( + ).

Patton extended this result to surfaces with a series of asperities and with a variety of surface profiles, he found that, at low values of N, sliding on the inclined surfaces occurred according to Eq. 1-20. Dilation of the specimens necessarily accompanied this mechanism. As the value of was increased above some critical value, sliding on the inclined asperity surfaces was inhibited, and a value of was eventually reached at which shear failure through the asperities occurred. Natural discontinuities rarely behave in the same way as this idealized model. However, the same two mechanisms -sliding on inclined surfaces at low normal loads and the suppression of dilation and shearing through asperities at higher normal loads are found to dominate natural discontinuity behavior. Generally, the two mechanisms are combined in varying proportions with the result that peak shear strength envelopes take the curved form of Figure 1-8b. Roughness can cause shear strength to be a directional property. Figure 1-10 illustrates a case in which rough discontinuity surfaces are prepared in slate specimens by fracturing them at a constant angle to the cleavage [START_REF] Brown | Shear strength characteristics of Delabole slates[END_REF]. When the specimens are tested in direct shear along directions of the ridges on the surfaces parallel to the direction of sliding (test A), the resulting shear strength envelope gives an effective friction angle of 22• which is comparable to the value of 19.5• obtained for clean and polished surfaces. However, when the shearing direction is normal to the ridges (test B), sliding up the ridges occurred with attendant dilation. A curved shear strength envelope is obtained with a friction angle of 67.5• at near zero effective normal stress and a friction angle of 46• at higher values of effective normal stress. [START_REF] Brown | Shear strength characteristics of Delabole slates[END_REF].

As it is possible to see from the last discussion, rock friction is a very complex phenomenon, it depends on many variables and involves different behaviors as the stick-slip and stable sliding.

Factors that can influence shear behavior of rock joints can be: matching of surfaces, gouge material, velocity of sliding, level of compression stress, elastic and fracture properties of the materials, asperities damage and sliding between asperities, time of contact between asperities of the joints surfaces, scale-effect, and surface roughness of joints, etc.

Models of joint friction

When studying the stability of rock masses, it is of capital importance to know the friction behavior of rock joints. As shown before, friction is a very complex phenomenon and is controlled by several parameters such as roughness, wear, scale effect, gouge material, boundary conditions (dilation allowed or not), speed of sliding, degree of matching, etc. To our knowledge, there is currently no model able to predict the shear behavior of rock joints and including all the latter variables involved in friction, only approximative solutions exists. To conclude with this chapter, some typical analytical models for joint's friction are going to be presented.

Newland and Alleyly (1957)

For rough surfaces, and assuming that the sliding surface consists of a series of sawtooth irregularities with an average angle , Coulomb Law is replaced by the following equation:

= tan(∅ + )
where ∅ is the friction angle (i.e., obtained for i=0).

1-21

Patton (1966)

Patton proposes a bilinear model for rough surfaces that account for two phenomena that have been observed experimentally: (i) overriding of asperities at low normal stress levels and (ii) shearing through asperities at higher normal stress levels:

= tan(∅ + ) when < 1-22 = + tan(∅ ) when > 1-23
where ∅ , ∅ , and are defined in Figure 1-11. For most practical purposes, ∅ = ∅ .

Barton-Bandis model (1974)

The effects of surface roughness on discontinuity deformation and strength have been taken into account by [START_REF] Bandis | Fundamentals of rock joint deformation[END_REF][START_REF] Barton | Strength, deformation and conductivity coupling of rock joints[END_REF] and [START_REF] Barton | Strength, deformation and conductivity coupling of rock joints[END_REF] in terms of a set of empirical relations between stress and deformation components and the parameters joint roughness coefficient, , and joint wall compressive strength, :

= tan + 1-24
where σ is the effective normal stress, is the joint roughness coefficient on a scale of 1 for the smoothest to 20 for the roughest surfaces, is the joint compressive strength and is the drained, residual friction angle. Eq. 1-24 suggests that there are three components of shear strength: a basic frictional component given by , a geometrical component controlled by surface roughness (

) and an asperity failure component controlled by the ratio ( / ). As shown in Figure 1-12, the asperity failure and geometrical components are combined to give the net roughness component, °. The total frictional resistance is then given by ( + °). Eq. 1-24 emphasizes that the shear strength of a rough joint is both scale dependent and stress dependent. Indeed, as increases, the term log 10 (JCS/ ) decreases, and so the net apparent friction angle decreases. As the scale increases, the steeper asperities shears off and the inclination of the controlling roughness decreases. Similarly, the asperity failure component of roughness decreases with increasing scale because the material compressive strength, JCS, decreases with increasing size. [START_REF] Bandis | Experimental studies of scale effects on the shear behavior of rock joints[END_REF].

The Barton-Bandis model takes explicit account of more features of discontinuity strength and deformation behavior than both the elementary Coulomb model and the Patton model discussed previously. However, its practical application may present some difficulties. In particular, the derivation of relations for the mobilization and degradation of surface roughness from a piecewise linear graphical format rather than from a well-behaved formal expression may lead to some irregularities in the numerical simulation of the stress-displacement behavior. Moreover, though three decades of experience has been gained in assigning values, the exercise remains a subjective one.

Archambault and Ladanyin's model (1972)

Archambaults and Ladanyi propose a model in which the shear strength yields:

= ̅ + tan ∅ (1 -) + 1 -(1 -) ̅ tan ∅ 1-25
where and ̅ are the proportion of total joint area sheared through the asperities and the rate of dilatancy at the peak shear strength, respectively. Both quantities are normal stress dependent and are such that:

= 1 -1 - 1-26 ̅ = 1 - tan 1-27
where 1 and 2 are empirical constants with suggested values of 1.5 and 4, respectively and is a transitional stress. The uniaxial compressive strength of the intact rock can be taken as an estimate of . is the shear strength of the rock compressing asperities and it can be described by any intact rock strength criteria.

Modified Ladanyi and Archambault Criterion (Saeb, 1990)

Saeb proposes a modification of Archambault and Ladanyi criterion such as:

= tan ∅ + (1 -) + 1-28
This criterion emphasizes the simultaneous contribution of shearing and sliding to the shear strength of a rock joint. In equation 1-28, is the proportion of joint surface area sheared through the asperities and (1 -) is the proportion on which sliding occurs. The relative contribution of shearing and sliding depends on the level of normal stress. At low normal stresses, sliding is dominant, at high normal stresses asperity shearing is dominant.

Swam and Zongqi Tribological model (1985)

Swam and Zongqi propose a tribological model derived from the extreme value statistics of a rough profile whose asperity peak height distribution ∅( ) is known and is the summed height = + of asperities of the upper and lower joints surfaces. Using a Hertzian contact formulation, normal load ( ) and shear force ( ) as a function of aperture is given by:

( ) = 4.739n b / AE (Z -e) / ∅(Z) dZ 1-29 ( ) = 9.870(nb) AS (Z -e) ∅(Z) dZ 1-30
where n is the peak density, b = mean peak radius, A = average area of contact, E'= Young modulus in plain strain and S = shear strength of the unit contact. Taking into account that the normal force ( ) and the shear force ( ) acts at a given inclination angle , the coefficient of friction of the joint at the peak is given by: = tan( + )

1-31 where is the average angle of friction between asperities and is the mean contact angle between asperities (i.e., or the average angle of the slopes of the asperities in the direction of shear).

Continuous-yielding joint model (1990)

The continuous-yielding joint model was designed to provide a coherent and unified discontinuity deformation and strength model, taking account of non-linear compression, non-linearity and dilation in shear, and a non-linear limiting shear strength criterion. Details of the formulation of the model are given by [START_REF] Cundall | Numerical simulation of fault instabilities with a continuouslyyielding joint model[END_REF].

The key elements of the model are that any shear displacement at a discontinuity has a component of plastic (irreversible) displacement, and any plastic displacement results in progressive reduction in the mobilized friction angle. The displacement relation is ∆ = (1 -)∆ where ∆ is an increment of shear displacement, ∆ is the irreversible component of shear displacement and is the fraction that the current shear stress constitutes of the limiting shear stress at the prevailing normal stress. The progressive reduction in shear stress is represented by:

∆∅ = - 1 (∅ -∅ )∆ 1-32
where ∅ is the prevailing mobilized friction angle, ∅ is the basic friction angle, and is a parameter with the dimension of length, related to joint roughness. The response to normal loading is expressed incrementally as:

∆ = ∆ 1-33
where the normal stiffness is given by:

= 1-34
in which and are further model parameters. The shear stress and shear displacement increments are related by:

∆ = ∆ 1-35
where the shear stiffness may also be taken as a function of normal stress, e.g.

= 1-36
in which and are further model parameters.

The continuously-yielding joint model has been shown to have the capability to represent satisfactorily single episodes of shear loading and the effects of cyclic loading in a manner consistent with that reported by [START_REF] Brown | Fatigue failure characteristics of some models of jointed rock[END_REF].

Leong and Randolph wear model (1992).

In Leong and Randolph model, the frictional resistance is attributed to three components. In addition to basic sliding and dilational component, they introduced a third component due to ploughing of the surface asperities and wear particles. The dilation and plough components are functions of the surface roughness, relative hardness of the two surfaces, normal stress and sliding distance. Degradation of surface roughness, and hence reduction in dilation and plough resistance, is formulated using wear theory. The model can reproduce not only the shear peak but also the post peak response as expressed by:

= tan(∅ + ) + 1 -tan ∅ tan 1-37
where ∅ is the sliding resistance between mineral grains at initial slip and when large shear displacements take place, ∅ is the residual angle of friction, is the angle of dilation and is the friction coefficient due to ploughing action on the surfaces by their asperities as well as for the wear debris.

Maksimovic (1996)

Maksimovic proposes the following model:

= tan ∅ + ∆∅ 1 + 1-38
where ∅ is the basic friction angle, i.e., the angle of the shearing resistance mobilized at high normal stress levels at which all dilatancy effects are suppressed, as all the asperities are sheared off forming the smooth shearing plane. This angle could be approximately equal to the angle of the physical friction between mineral grains. In Eq. 1-38 ∆∅ corresponds to the joint roughness angle, reflecting the surface roughness of the discontinuity and the associated dilatancy effects at zero stress level. ∆∅ can be described as the angle of maximum dilatancy which occurs on an undamaged rugged surface, and "the median angle pressure" is equal to the value of normal stress at which the contribution of dilation and the breakage of asperities is equal to one half of the angle of dilatancy for zero normal stress. It reflects mainly the strength and rigidity of asperities that make the discontinuity surface rough.

Kwafniewski and Wang fractal model (1997)

Kwafniewski and Wang proposed the following model for the shear peak resistance of a joint based on the fractal geometry that characterizes the joint surface roughness:

= √ + 1-39
where D is the fractal dimension of the joint surface (self-similar fractal surface), is the tensile strength of the rock material, is the topothesy of the joint surface, and is the height variance of the joint surface defined by:

= ( ) 1-40
where r c is the lag at the correlation length, and c 1 and c 2 are constants.

2 Morphology of rock joint surfaces

Overview

The rocks forming the earth crust are submitted to mechanical stresses due to natural process such as the movement of tectonic plates or the effect of gravity. Because of these stresses, rocks can fracture and produce sets of rock joints or discontinuities as shown in Figure 2-1, that illustrates an example of discontinuous rock masses. It is nowadays well known that geometrical properties of joints like their spacing, orientation, persistence and roughness, combined with elasto-plastic and fracture behavior of the intact rock, control the hydro-mechanical behavior of the rock-masses, [START_REF] Hudson | Engineering rock mechanics: an introduction to the principles[END_REF]. In practice, it is common to study the mechanical behavior of one single joint in order to understand and predict the behavior of rock masses, which usually contains various sets of joints.

This chapter is dedicated to the characterization of the morphology of rock joints. The tools used regularly in quantitative fractography are presented and the statistical parameters commonly used to characterize the surface roughness are particularly detailed. At the end of the chapter, a numerical algorithm is proposed to generate self-affine surfaces by means of the control of different statistical parameters such as the variance of heights, the self-affine correlation length and the roughness exponent (Hurst's exponent H).

Roughness description of rock joints (Fracture surfaces)

The surface roughness is a measure of the texture of a surface. It is quantified by the vertical deviations of a real surface from its ideal form. If these deviations are large, the surface is rough; if they are small, the surface is smooth. Roughness plays an important role in determining how a real object will interact with its environment. For instance, rough surfaces have higher friction coefficients than smooth surfaces and so exhibit a quicker wear compared to smooth surfaces.

Roughness can be quantified either from a profile (line) or from a surface (area). There are many different roughness parameters in use. Nevertheless, most parameters reduce all of the information to a single value and, hence, great care must be taken in applying and interpreting them. For instance, small changes in how the raw profile data are filtered or how the average line (2D)/average plane (3D) is calculated can greatly affect the calculated parameter.

There are different ways to record the topography of a surface (laser profilometer, atomic force microscopy, confocal microscopy, stereovison from scanning electron microscopy, etc.) usually, the height distribution of the surface is measured from a regular grid (in the average plane of the surface) and this leads to a set of point coordinates xyz of the surface.

Average plane

One of the most common parameters of a rough surface is the average mean plane, obtained by adjusting the equation:

= + + 2-1
to the xyz coordinates of the surface. Usually, the average plane is computed and subtracted to the raw data before the characterization itself of the joint morphology. For instance, the least square method can be used to estimate the coefficients a, b and c of the equation defining the average plane (Eq. 2-1).

Amplitude parameters

The amplitude parameters characterize the surface roughness on the basis of the vertical deviations of the profile from the average plane or profile. They represent the most common surface roughness parameters found on mechanical engineering and in the technical literature.

The first parameter is the arithmetic average of the absolute values of the heights ( , ). The average roughness is expressed in units of height:

= 1 | ( , )| 2-2
Unfortunately does not discriminate very well between profiles of quite different morphology. Another amplitude parameter is the root mean square average of profile heights ( , ) defined as:

= 1 ( , ) 2-3 2.2.2.1. Height distribution
Generally, a rough surface can be seen as the result of a stochastic process where ( , ) is a random variable representing the height of the surface at the point ( , ). The random variable ( , ) is completely specified by the range of values and the probability ( ) for each value. Thus, if one notes ( ) the probability of a height ranged from to + , then the cumulative probability that a height is lower than some level ℎ is:

(ℎ) = ( ) 2-4
On the other hand, the probability distribution ( ), whose average is at the height of the mean plane of the profile, may be characterized by its central moments:

= ( ) 2-5
The second moment , usually noted as , is the height variance and describes the scattering of the distribution from its average. It is usual to use the square root of the variance in order to obtain a quantity expressed in units of height. The square root of the variance corresponds to the standard deviation of the distribution and is formally identical to the roughness.

Skewness

The third moment characterizes the skewness and is usually noted as . The skewness corresponding to an odd moment, its measure is sensitive to the degree of asymmetry of the height distribution.

Kurtosis

The fourth moment characterizes the kurtosis and describes the flattening or more exactly the overall shape of the distribution (for instance, a Gaussian distribution has a kurtosis of 3).

Texture parameters and anisotropy

The average parameters such as or contain no information about the spatial or the textural variation of the profile, i.e., about the height variation from point to point of the profile. The description of a surface restricted to its variations in amplitude leads only to partial information of its morphology. The change in surface height as a function of position or the length scale of observation is also very important in the description of roughness. Indeed, as it will be shown in the subsequent chapters, the texture of the surface can influence strongly the mechanical behavior of rock joints. One tool to characterize the texture of a surface it is the so called correlogram. The basis of this method consists in investigating the correlation between pairs of points on a surface as a function of the distance between them.

In this respect, the power spectral density function (PSDF) that corresponds to the Fourier transform of the autocovariance function (ACVF) is commonly used to study the height correlation between the points of a surface. For a profile of length the autocovariance function ( ) is defined as:

( ) = 1 - ( ) ( + ) , 2-6
where ( ) and ( + ) are pairs of heights separated by a distance (Figure 2-2). The ACVF is very sensitive to periodic components of the surface and will detect these even if they are obscured by random components. If the points of the profile are entirely uncorrelated, the ACVF will decay asymptotically at a rate that depends on height variances of the surfaces. For natural rough surfaces (fracture surfaces for example), there exist a maximum separation distance = , above which there is no more statistical correlation between the heights.

The corresponding function of the ACVF in the frequency domain is the power spectral density function (PSDF). The Wiener-Khinchine relation defines the PSDF as the Fourier transform of the ACVF:

( ) = ( ) , 2-7
where is an angular frequency.

A single value of the PSDF has dimensions of height squared per unit frequency, and as for the height distribution, it is possible to express moments as:

= ( ) 2-8
The zero moment corresponds to the variance of the height distribution. Thus, if the power spectrum is integrated over a pass-band of surface frequencies ranged from 1 to 2, the area under the function's curve corresponds to the height variance of the surface. The second and r z z(r) z(r + τ)

fourth moments are the variances of the distributions of the slopes and curvatures of the surfaces respectively.

Commonly, instead of the ACVF ( ), the structure function ( ) is used:

( ) = 1 - [ ( + ) -( )] 2-9
or in discrete form:

( ) = 1 ( -) 2-10a
The structure function is related to the ACVF by:

( ) = 2 [1 -( )] ( ) = ( ) (0)

2-11

In geostatistics, instead of the structure function, a function called the variogram ( ) is commonly used. The structure function ( ) and the variogram ( ) are very similar, in fact, they are related by a factor of 2 ( ( ) = ( )/2 ).

Both the structure function ( ) and the variogram ( ) contain the same information as the ACVF and the PSDF, but offer some practical advantages: they are stable and easy to compute, they do not impose a periodogram model of the surface, and they do not require prior-high-pass filtering.

Another useful property is that for a random stationary profile ( ) → 2 (or ( ) → ) as → ∞.

The variogram or structure function can exhibit functional changes more clearly than the ACVF.

Hurst exponent

For some surfaces, such as most of the surfaces produced by fracture, their PSDF can be represented by a relation of the type:

( ) = 2-12
where is a dimensionless constant and is a constant with dimensions of length called the topothesy. These constant appear to be an intrinsic property of the surface; when a surface is described by Eq. 2-12, it is possible to express most of the statistical parameters discussed before as function of and . Mandelbrot [START_REF] Mandelbrot | How long is the coast of Britain?. Statistical self-similarity and fractional dimension[END_REF] pointed out that such surfaces are examples of fractals.

One characteristic of fractals is that they are continuous functions but not differentiable. They possess the property of self-similarity, that is, they appear the same at any scale of magnification. Self-similar fractals can be characterized by a single parameter, the Hurst exponent . Nevertheless, for some fractal objects, when the scale of magnification is changed, a scaling factor with dimensions of length must be introduced to restore the appearance of the object at lower magnifications. This scaling factor turns out to be the topothesy, and fractals of this kind are described as self-affine, [START_REF] Bouchaud | Scaling properties of cracks[END_REF]. Examples of self-affine fractals in nature include fracture surfaces [Metals; Pande (1987), rocks;[START_REF] Odling | Natural fracture profiles, fractal dimension and joint roughness coefficients[END_REF], [START_REF] Schmittbuhl | Scaling invariance of crack surfaces[END_REF], [START_REF] Power | Topography of natural and artificial fractures in granitic rocks: implications for studies of rock friction and fluid migration[END_REF], [START_REF] Develi | Quantification of natural fracture surfaces using fractal geometry[END_REF], [START_REF] Mourot | Anomalous scaling of mortar fracture surfaces[END_REF], [START_REF] Yang | On the applicability of self-affinity concept in scale of three-dimensional rock joints[END_REF][START_REF] Sakaguchi | Asperity height aperture of an artificial tensile fracture of metric size[END_REF], [START_REF] Dougan | Estimating the cut-off in the fractal scaling of fractured concrete[END_REF], [START_REF] Carpinteri | Three-dimensional fractal analysis of concrete fracture at the meso-level[END_REF]] and natural terrain [START_REF] Mark | Scale-dependent fractal dimensions of topographic surfaces: an empirical investigation, with applications in geomorphology and computer Mapping[END_REF], [START_REF] Hirata | Fractal dimension of fault systems in Japan: fractal structure in rock fracture geometry at various scales[END_REF], [START_REF] Murata | The variogram method for a fractal model of a rock joint surface[END_REF], [START_REF] Ehlen | Fractal analysis of joint patterns in granite[END_REF], [START_REF] Renard | Three-dimensional roughness of stylolites in limestones[END_REF], [START_REF] Ebner | Anisotropic scaling of tectonic stylolites: A fossilized signature of the stress field?[END_REF]].

An effective way to estimate the fractal parameters of a surface is to compute the structure function (or the variogram). For a fractal profile, it can be shown that the structure function obeys the following relation, [START_REF] Costa | Fractal description of rough surfaces for haptic display[END_REF]:

( ) = 2-13
In other words, the structure function of a fractal profile obeys a power law, so it plots as a straight line on a log-log scale . This is an easy way of establishing the fractal behavior, and from the slope and intercept of this straight line both, the Hurst exponent and the topothesy can easily be calculated.

Anisotropy

Some surfaces present different roughness parameters when measured in different directions [START_REF] Power | Topography of natural and artificial fractures in granitic rocks: implications for studies of rock friction and fluid migration[END_REF], [START_REF] Ponson | Ph.D Thesis, Crack propagation in disordered materials: How to decipher fracture surfaces[END_REF], [START_REF] Lopez | Anomalous scaling of fracture surfaces[END_REF], [START_REF] Renard | High resolution 3D laser scanner measurements of a strike-slip fault quantify its morphological anisotropy at all scales[END_REF], [START_REF] Mourot | Scaling properties of mortar fracture surfaces[END_REF], [START_REF] Candela | Characterization of fault roughness at various scales: implications of threedimensional high resolution topography measurements[END_REF]]. This anisotropy of the roughness development can be put in evidence when the Hurst exponent , the height variance or the correlation length exhibit change as a function of the considered direction in the surface. For instance, there are experimental evidences that the Hurst exponent (also called as roughness exponent) computed for profiles parallel to the direction of the crack propagation takes a different value from profiles measured perpendicular to the crack propagation direction (Figure 2-4) , [START_REF] Ponson | Ph.D Thesis, Crack propagation in disordered materials: How to decipher fracture surfaces[END_REF]. Values of the roughness exponents of 0.8 (measured for profiles parallel to the crack front direction) and 0.6 (perpendicular to the crack front direction) have been computed for fracture surfaces of different materials and for various modes of fracture; the existence of such values of the roughness exponent for several materials and modes of fracture gave the idea to some authors to say that such values are universal [START_REF] Bouchaud | Scaling properties of cracks[END_REF].

2σ² 2H Log(τ) Log [S(τ)] σ² 2H τ γ(τ)
On the other hand, it has been shown recently, [START_REF] Morel | Transient damage spreading and anomalous scaling in mortar crack surfaces[END_REF], from profiles recorded in the direction parallel to the crack front in a mortar, that another roughness exponent value H=0.4 could be observed at large length scales. This roughness regime could be characteristic of the length scales for which the material can be considered as linear elastic, [START_REF] Morel | Transient damage spreading and anomalous scaling in mortar crack surfaces[END_REF]. Thus, if any consideration of the profile direction is made, the roughness exponent H value is ranged between 0.4 and 0.8 In addition, fracture surfaces of some kind of rocks such as shale or slate, among others, exhibit anisotropy in the height variance or correlation length for the reason that such materials breaks in some preferential planes (Figure 2345). Moreover, as the rock joints normally are exposed to many years of physical or chemical weathering, anisotropic rough joint surfaces can be produced by preferential weathering caused by groundwater, for example.

Bearing area (Abbott's curve)

The Abbott-Firestone curve was first described by E.J. Abbott and F.A. Firestone in 1933 and describes the surface texture of an object. The curve can be found from a profile trace by drawing lines parallel to the base line and measuring the fraction of the line which lies within the profile (Figure 23456). Mathematically, it is the cumulative probability density function of the surface profile's height and can be calculated by integrating the profile trace. It is useful for understanding the properties of sealing and bearing surfaces. The Abbott's curve concept will be used later in this work to identify the self-affine surface generated at the inner of the discrete element model.

Directional Parameters

Z2, Z3

Other common parameters used to describe the morphology of rough surfaces are the directional parameters and which quantifies the mean slope and mean curvature of a profile. Normally, these parameters are computed for various directions in the surface in order to put in evidence any anisotropy.

These parameters are classically defined as:

( ) = 1 - 2-14 x z z z x z z b) Probability density function z c) Bearing area (Abbot's curve) a) Profile trace ( ) = 1 -2 + 2-15
where N is the number of discrete points of the rough profile and is the separation between points, i.e., the sampling step of the surface or the calculation step (which not necessarily coincides with the sampling step). Note that for a self-affine profile (i.e., rock fracture surfaces), ( ) and ( ) are related with ( ) in the following way:

( )² = ( ) ² = ² 2-16a ( ) = 2- 17
( )² = ( ) = 2-18a ( ) = 2-19
Note from Eq. 2-15 and Eq. 2-16 that for a self-affine surface, as increases, ( ) and ( ) tends to zero.

Colatitudes 2D

The parameters and synthesize the slope and curvature of the surface in a single average value. It is possible to go a step further in the roughness description and compute the slopes between the discrete adjacent points of the surface and obtain their distribution. For this, a new parameter called colatitude 2D is defined:

= - ∆ 2-20
where ∆ is the step of sampling or more generally the distance between both points and gives the slopes between neighboring points. Usually is computed for different directions in the surfaces and separated in positive and negative values, [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF].

As explained before in this chapter, the roughness of rock joints is fractal (more specifically, a fracture surface which has been freshly created in a rock is expected to be self-affine) and because of different reasons like weathering, rock microstructure or by the simple fact of fracturing, the roughness of the rock joint can be isotropic or anisotropic. So in order to study the influence of the joint roughness in the mechanical response of rock joints it is necessary to be able of generating isotropic or anisotropic self-affine surfaces controlling the roughness exponent, correlation length and height variance. The rest of this chapter is dedicated to present a methodology to numerically generate this kind of surfaces.

Numerical method to generate fractal surfaces

Let us consider the rock joint surface as a stochastic process ( ). Additionally, the surface is second-order stationary with zero-mean and a well defined variance. In this case, the stochastic process ( ) can be transformed in an equivalent process ( ) in the frequency domain by:

( ) = ( ) 2-21
In Eq. 2-21, ( ) is the Fourier transform of ( ) and both ( ) and ( ) contain the same information but expressed in different ways. ( ) decomposes the content of ( ) into different frequencies present in that process and helps to identify periodicities. As we have seen before in this chapter, the autocovariance function is usually used to study the correlation between points in the surface and has an equivalent function in the frequency domain called the Power Spectral Density Function. The so-called Wiener-Khintchine theorem states that any second order stationary process has an autocovariance function ( ) of the form:

( ) = ( ) 2-22
where ( ) is the power spectral density function defined by:

( ) = | ( )| 2-23
Conversely, if one considers a function represented by Eq. 2-22, there exists a stationary process ( ) with ( ) as autocovariance. Both functions, covariance and spectral density, contain the same information (Figure 234567). Thus, in order to generate numerically a self-affine surface, we propose to use the scaling property of the PSDF. For this we consider firstly a discrete random field that is a linear combination of sines and cosines with random phases and random amplitudes. In our case, the amplitudes are important and the phases are of no interest for the moment. As it is well known by Fourier analysis, a sequence ( ) of N points can be transformed into a finite set of Fourier coefficients:

( ) = 2 + 2 2-24 where, = 1 ( ) cos 2 2-25 = 1 ( ) sin 2 2-26
and, where N corresponds to the number of points of the sequence (number of points of the profile) and {j = 0 . . . . , N -1},.

Another way of representing the finite discrete sequence is as a complex exponential Fourier series:

( ) = ( ) / 2-27
where

( ) = ( ) ( )

2-28

corresponds to the discrete amplitude and, i = (-1), and {r = 0,…….., N-1}.

Coefficients ( ) are related to and coefficients by the equations:

( ) = + 2-29 ( ) = - 2-30
The discrete amplitude spectrum is the representation of amplitude ( ) vs frequency 2πj/N and the phase spectrum is the plot of phase ( ) vs frequency 2πj/N. Moreover, the Parseval-Rayleigh theorem states that the sum of squared amplitudes is equal to the total power of {z(r), r = 0 . . . . . N -1 } that can be identified with the height variance of the surface :

( ) = 1 | ( )| = 2-31
and, the discrete amplitude spectrum can be related directly to the discrete spectral amplitude:

( ) = ( ) 2-32
As shown previously, this result can be also obtained from the Fourier transform of the autocovariance function (Eq. 2-7).

On the other hand, the spectral density function of a real stochastic process is an even function:

( ) = (-) 2-33
and, in this case, the discrete amplitude spectrum can be written as:

( ) = (-) 2-34
Then, the amplitude spectrum has to be related to the spectral density that depends on the covariance model that we wish to impose on the stochastic process. The phase spectrum does not affect the covariance structures, then it can be taken at random from a uniform distribution between 0 and 2π:

( ) ~ (0, 2 ) 2-35
In this way, we construct the complex Fourier coefficients:

( ) = ( ) ( ) 2-36 ( ) = ( ) ( ) - ( ) ( ) 2-37 : ( ) = ( ) ( ) 2-38 : ( ) = ( ) ( ) 2-39 ( ) = ( ) - ( ) 2-40
The Fourier transform of a real stochastic process is an Hermitian function, then coefficients ( ) must be hermitian, i.e., with an even real part:

( ) = (-) 2-41
and an odd imaginary part:

( ) = - (-) 2-42
By calculating the inverse Fourier transform of the complex coefficients ( ), the discrete finite realization {z(r), r = 0 . . . . . N -1} is obtained with the specified covariance model:

( ) = ( ) / 2-43
If the number of points N is a power of two, the inverse discrete Fourier transform can be rapidly and efficiently computed with the fast Fourier transform (FFT).

For a real process, the Fourier coefficients must have some symmetries which are specified as follows:

(i) for a 1D real process, the complex Fourier coefficients must exhibit a symmetry as shown in Figure 2345678.

1 N/2 +1 N η χ ϑ χ* Figure 2-8 Complex Fourier coefficient representation for 1D real process.
In Figure 2-8, η and ϑ correspond to real coefficients (the imaginary parts are equal to 0) while χ are complex coefficients and χ * correspond to the complex conjugate of χ.

(ii) for a 2D real process, the symmetries of the complex Fourier coefficients are shown in Figure 23456789.

The coefficients of sectors corresponding to the different letters are independent, as well as the coefficients of a single sector are independent between themselves. Sectors A, N, S and X correspond to real coefficients (the imaginary parts are equal to zero) and the rest are complex sectors.

1 N 1 /2 + 1 N 1 1 A B X B* D E P G N 2 /2 + 1 N V S V* N 2 D* G* P* E* Figure 2-9 Complex Fourier symmetries in 2D u v 2.3.1. Procedure to generate self-affine surfaces a)
Power spectrum density of a self-affine surface.

First, it is necessary to generate a unitary matrix U m of size L x L. Here, L 2 will correspond to the number of discrete points of the resulting surface. By convenience, L must be an odd number. U m represents the power spectral density ( , ) of the surface and therefore must obey a power law behavior in the form:

( , ) ∝ 1 √ + 2-44
where = 2 + 2, and and are the frequencies. Note that the center of the ( , ) matrix is given by + 1, + 1 and must have a value equal to 0, i.e., ( = 0, = 0) = 0, in order that the average height of the surface is equal to zero. Moreover, the Fourier transform of a discrete surface ℎ( , ) is defined as:

( , ) = 1 ℎ( , ) ( / / )

2-45

The power spectral density of the surface is related to the Fourier coefficients ( , ) through the expression:

( , ) = ( , ) = 1 ℎ( , ) ( / / )

2-46

Note that each value of the matrix ( , ) represents the power spectral density of the surface in the direction = ( / ) and at the frequency = √ + as shown in Figure 2-10.

Making ( = 0, = 0) = 0 in Eq. 2-46 we get:

( = 0, = 0) = 1 ℎ( , ) 2-47 1 ℎ( , ) = 〈ℎ( , )〉 = 0 2-48 ( = 0, = 0) = |〈ℎ( , )〉| = 0 2-49
〈ℎ( , )〉 stands for the average value of ℎ. Assigning ( = 0, = 0) = 0 makes easier to adjust the height variance of the surface.

b) Imposition of the power law behavior

A self-affine surface must have a power spectral density ( , ) that follows a power law behavior, [START_REF] Makse | Method for generating long-range correlations for large systems[END_REF]:

( ⃗) ∝ 1 | ⃗| ( ) 2-50
where | | = √ + , ( ) = 2 ( ) + 2 where ( ) corresponds to the Hurst exponent (or the roughness exponent) in the direction = (see Figure 2-10). If ( ) = whatever the surface is isotropic in roughness exponent while, if ( ) evolves as a function of the direction , the surface is anisotropic. In this work, the function ( ) is described from the following expression:

( ) = 1 cos + sin H 2-51
where and correspond to the roughness exponents in the directions x and y respectively. Both directions are perpendicular between themselves and, the x axis corresponds to = 0 while the y axis is related to = π/2.

c) Fixing the correlation length

As previously mentioned, the structure function ( ) of fracture surfaces (Figure 2-3) exhibit usually an horizontal asymptotic behavior when the separation distance between two points of the surface reaches a critical value . The critical distance is called self-affine correlation length and, for distances > there are no more statistical correlation between the points of the surfaces while for distances < the height fluctuations of the surfaces exhibit a self-affine scaling invariance. The self-affine correlation has an equivalence in the frequency domain noted by such as :

= 1 2 2-52 P(u,v) u v θ
Fracture surfaces can have a self-affine correlation length (or in the frequency domain) that depends on the considered direction in the surface. This anisotropy in terms of self-affine correlation length or in terms of the corresponding frequency can be generated from an analogous approach to the one considered from Eq.2-51 for the roughness exponent H, which consists in considering an evolution of the frequency as a function of the angle such as :

( ) = 1 cos + sin 2-53
where and are the self-affine correlation frequencies in the directions x and y respectively.

Thus, using the Eq. 2-53 the power spectral density of the fracture surfaces is determined as follows:

( ⃗) ∝ | ⃗| ≤ | ( )| 1 | ⃗| ( ) | ⃗| > | ( )| 2-54 d)
Fixing the variance of the surface Finally, fracture surfaces can exhibit also a variance depending of the considered direction in the surface. This anisotropy in terms of variance can be generated from the approach previously used for the Hurst exponent (Eq.2-51) and the self-affine correlation frequency (Eq. 2-53). Indeed, if one considers an arbitrary value of variance noted as , its corresponding PSD can be expressed from the relationship = 1 ∑ . Let us now consider a particular value of variance noted as and its corresponding PSD , the relationship between the previous arbitrary values of the variance and the PSD and these particular ones can be expressed as :

= 2-55
On this basis, according to Eq.2-51 or Eq. 2-53, one can consider an evolution of the variance as a function of the angle and especially its correspondence expressed in terms of power spectrum density such as :

( ) = cos + sin 2-56
where and correspond to the expected values of the variance in the directions x and y respectively.

e) Self-affine surface in the real space (Inverse Fourier Transform) ( , ) is used to obtain the complex Fourier coefficients ( , ). The Fourier coefficients are defined as follows:

( , ) = ( , ) + ( , ) 2-57 ( , ) = ( , )cos [ ( , )] 2-58 ( , ) = ( , )sin [ ( , )] 2-59
where ( , ) correspond to the random phase and considered in a first step as uniformly distributed between 0 and 2π Finally, the inverse discrete Fourier transform of the matrix of complex coefficients ( , ) allows to obtain the rough surface with the desired power spectral density:

ℎ( , ) = ( , ) ( / / )

2-60

This is done rapidly and efficiently by the fast Fourier transform. Examples of self-affine surfaces generated with the procedure mentioned before are shown in Appendix 1.

Discrete elements model of an isotropic rock

Overview

The discrete element method (DEM) is any of a family of numerical methods for computing the motion and effect of a large number of particles. Even if DEM is very closely related to molecular dynamics, the method is generally distinguished by its inclusion of rotational degrees-of-freedom and often complicated particle geometries (including polyhedra). The general method was originally developed by Cundall in 1971 to solve problems in rock mechanics where the fundamental assumption is that the material consists of separate, discrete particles that interact with each other when they are in contact. These particles are considered as rigid bodies and may have different shapes, properties and obey Newton's laws of motion. Examples of materials modeled by DEM are: liquids and solutions, bulk materials in storage silos, such as cereals, granular matter and blocky or jointed rock masses.

Today DEM is becoming widely accepted as an effective method for addressing engineering problems in granular and discontinuous materials, especially in granular flows [START_REF] Fillot | Modelling third body flows with a discrete element method-a tool for understanding wear with adhesive particles[END_REF], powder mechanics Jauffrès ( 2012), and rock mechanics [START_REF] Fakhimi | Simulation of failure around a circular opening in rock[END_REF], [START_REF] Diederichs | Damage initiation and propagation in hard rock during tunneling and the influence of near-face stress rotation[END_REF], [START_REF] Wang | Modeling Lac du Bonnet granite using a discrete element model[END_REF], [START_REF] Shimizu | Distinct element analysis for class II behavior of rocks under uniaxial compression[END_REF].

Recently, the method was expanded to take into account thermal phenomena [START_REF] Wanne | Bonded-particle modeling of thermally fractured granite[END_REF] and it is possible to make hybrid models by coupling DEM-FEM methodologies [START_REF] Cai | FLAC/PFC coupled numerical simulation of AE in large-scale underground excavations[END_REF].

Discrete element methods are relatively computationally intensive but, take advantage of parallel processing capabilities to scale up the number of particles or length of the simulation. Also, the calculation method is a time-stepping explicit scheme. The solution scheme is identical to that used by the explicit finite-difference method for continuum analysis. The DEM is based upon the idea that the time-step chosen may be so small that, during a single time-step, disturbances cannot propagate further from any particle than its immediate neighbors. Then, at all times, the forces acting on any particle are determined exclusively by its interaction with the particles with which it is in contact.

Since the speed at which a disturbance can propagate is a function of the physical properties of the discrete system, the time-step can be chosen to satisfy the above constraint. There are several advantages to such a scheme, it makes possible to simulate the nonlinear interaction of a large number of particles without excessive memory requirements or the need for an iterative procedure, since matrices are never stored. DEM may be used to model static or dynamic problems, but the full dynamic equations of motion are solved even when static solutions are required. This is done in order to follow such phenomena as failure and "flow" of material in a realistic manner; it is not necessary to invoke some nonphysical algorithm, as it is done in some implicit methods.

A DEM-simulation is started by first generating a model, which results in spatially orienting all particles and assigning an initial velocity. Then the model executes many thousands of time-steps and at each time-step, Newton's second law (force = mass × acceleration) is integrated twice for each particle to provide updated velocities and new positions. Based on these new particle positions, contact forces are derived from the relative displacements between particles and are used to compute the new positions and velocities of all particles for the next time-step, this loop is repeated until the simulation ends.

When modeling materials at the meso-macro scale the typical forces involved during simulation are frictional forces, contact forces between particles and the effect of gravity.

In the present work a soft-contact approach is used in which the particles are treated as rigid bodies allowed to overlap at the contact points. The magnitude of the inter-particle forces depends on the degree of overlap between the particles. In addition to spherical particles, various of the DEMs codes also includes "walls." Walls allow one to apply velocity boundary conditions to assemblies of balls for purposes of compaction and confinement. The balls and walls interact with one another via the forces that arise at contacts. The equations of motion are satisfied for each ball. However, the equations of motion are not satisfied for each wall (i.e., forces acting on a wall do not influence its motion). Instead, its motion is specified by the user and remains constant regardless of the contact forces acting on it. Also, contacts may not form between two walls; thus, contacts are either ball-ball or ball-wall.

Despite the fact that the general DEM method allows to simulate particles of complicated shapes, in the present work spherical particles are used since contact detection between spherical objects is much simpler than contact detection between angular objects, moreover, it is also possible to create particles of arbitrary shape by attaching two or more spherical particles together such that each group of particles acts as an autonomous object, if the bonded particles in a group are not allowed to break out, then it will be called a clump. The spheres within a clump may overlap to any extent, and each clump behaves as a rigid body with deformable boundaries, thus, can serve as a super-particle of general shape. It is also possible to model a brittle solid by bonding every particle to its neighbors; the resulting assembly can be regarded as a "solid" that has elastic properties and is capable of "fracturing" when bonds break in a progressive manner.

When a material can be considered as a continuum, a common technique to model it is to create a grid and install the initial stress and boundaries. In a discrete element code, a compacted state cannot be pre-specified in general, since there is no unique way to pack a number of particles within a given volume. A process analogous to physical compaction must be followed until the required porosity is obtained. The initial stress state cannot be specified as independent of the initial packing since contact forces arise from the relative positions of particles. As we will see at the end of this chapter, these internal forces can influence the fracture behavior of the model so it becomes important to know and control the internal stresses in the model as a function of the porosity.

This chapter resumes the principal characteristics of the Discrete Element Method (DEM) as implemented in the software PFC-3D (Itasca C. G. Inc.). The main elements of the method will be explained and a DEM model of a mortar specimen is presented. The elastic properties of the DEM model and its fracture behavior will be obtained by the calibration of a representative elementary volume (REV) composed of discrete particles. The software PFC-3D developed by Itasca, Inc. was used in this thesis to perform all the simulations on a PC with 2 Processor Intel Xeon CPU E5-2640 0 @ 2.50 GHz with 64.0 Go of RAM and Windows Server 2012, 64 bits. To see more details about the DEM method used in this thesis, the author suggests to consult the User Manuals of PFC-3D, Itasca, C.G. Inc.

DEM general formulation

Notion of a single particle

The motion of a single rigid particle is determined by the resultant force and moment of force vectors acting upon it, and can be described in terms of the translational motion of a point in the particle and the rotational motion of the particle. The translational motion of the center of mass is described in terms of its position, , velocity, ̇ , and acceleration, ̈ ; the rotational motion of the particle is described in terms of its angular velocity, , and angular acceleration, ̇ . The equations of motion can be expressed as two vector equations: one relates the resultant force to the translational motion; and the other relates the resultant moment to the rotational motion. The equation for translational motion can be written in the vector form:

= ( ̈-) 3-1
where is the resultant force, i.e, the sum of all externally applied forces acting on the particle, m is the total mass of the particle; and is the body force acceleration vector (e.g., gravity loading).

The equation for rotational motion can be written in the vector form:

= ̇ 3-2
where is the resultant moment acting on the particle, and is the angular momentum of the particle. This relation is referred to a local coordinate system that is attached to the particle at its center of mass. If this local system is oriented such that it lies along the principal axes of inertia of the particle, then Eq. 3-2 reduces to Euler's equation of motion:

= ̇ + ( -) 3-3 = ̇ + ( -) = ̇ + ( -)
where , , and are the principal moments of inertia of the particle, ̇ , ̇ , and ̇ are the angular accelerations about the principal axes, and , , and are the components of the resultant moment referred to the principal axes. For a spherical particle of radius R, whose mass is distributed uniformly throughout its volume, the center of mass coincides with the sphere center. Any local-axis system attached to the center of mass is a principal-axis system, and the three principal moments of inertia are equal to one another. Thus, for a spherical particle, Eq. 3-3 can be simplified and referred to the global-axis system as:

= ̇ = 2 5 ̇ 3-4
The equations of motion, given by Eq. 3-1 and Eq. 3-4, are integrated using a centered finite difference procedure involving a time-step ∆ . The quantities ̇ and are computed at the midintervals of ± ∆ /2, while the quantities , ̈ , ̇ , , and are computed at the primary intervals of ± ∆ . The following expressions describe the translational and rotational accelerations at time t in terms of the velocity values at mid-intervals. The accelerations are calculated as:

̈ ( ) = 1 ∆ ̇ ( ∆ / ) -̇ ( ∆ / ) 3-5 ̈ ( ) = 1 ∆ ̇ ( ∆ / ) -̇ ( ∆ / )
Inserting these expressions into Eq. 3-1 and Eq. 3-4, and solving for the velocities at time ( + ∆ /2), results in:

̇ ( ∆ / ) = ̇ ( ∆ / ) + ( ) + ∆ 3-6 ( ∆ / ) = ( ∆ / ) + ( )

∆

Finally, the velocities in Eq. 3-6 are used to update the position of the particle center as:

( ∆ ) = ( ) + ̇ ( ∆ / ) ∆ 3-7
The calculation cycle for the law of motion can be summarized as follows. Given the values of ̇ ( ∆ / ) , ( ∆ / ) , ( ) , ( ) and ( ) , Eq. 3-6 is used to obtain ̇ ( ∆ / ) and ( ∆ / ) . Then, Eq.

3-7 is used to obtain ( ∆ ) . The values of ( ∆ ) and ( ∆ ) , to be used in the next cycle, are obtained by application of the force displacement law.

Contact model

Force-Displacement Law

The force-displacement law relates the relative displacement between two entities (ball-ball and ballwall) at a contact to the force acting on the entities. For ball-ball contact, an additional force and moment arising from the deformation of a "cement material" represented by a parallel bond can also act on each particle. Only the calculation of the forces arising from the contact at a point is described in this section. The force-displacement law can be described in terms of a contact point, [ ] , lying on a contact plane that is defined by a unit normal vector, . The contact point is within the interpenetration volume of the two entities. For ball-ball contact, the normal vector is directed along the line between ball centers; for ball-wall contact, the normal vector is directed along the line defining the shortest distance between the ball center and the wall. For ball-ball contact, the unit normal vector, , that defines the contact plane is given by:

= 3-8
where and are the position vectors of the centers of balls A and B, and d is the distance between the ball centers:

= - = - - 3-9 Contact plane A B x z y d R A R B P A P B Wall
For ball-wall contact, is directed along the line defining the shortest distance, d, between the ball center and the wall. This direction is found by mapping the ball center into a relevant portion of space defined by the wall. The idea is illustrated in Figure 3-2 for a two-dimensional wall composed of two line segments, AB and BC. All space on the active side of this wall can be decomposed into five regions by extending a line normal to each wall segment at its end points. If the ball center lies in region 2 or 4, it will contact the wall along its length, and will be normal to the corresponding wall segment. However, if the ball center lies in region 1, 3, or 5, it will contact the wall at one of its end points, and will lie along the line joining the end point and the ball center. For a three-dimensional convex polygonal walls, the above idea is extended such that the ball may contact the wall at a vertex, along an edge joining two wall segments, or on a face. The overlap , defined to be the relative contact displacement in the normal direction, is given by:

= + - ( - ) - ( - ) 3-10
where R is the radius of the ball. The location of the contact point is given by:

= + - 1 2 ( - ) + - 1 2 ( - )

3-11

The contact force vector (which represents the action of ball A on ball B for ball-ball contact, and the action of the ball on the wall for ball-wall contact) can be resolved into normal and shear components with respect to the contact plane as

= + 3-12
where and denote the normal and shear component vectors, respectively. The normal contact force vector is calculated by

= 3-13
where is the normal stiffness [force/displacement] at the contact. The value of is determined by the current contact-stiffness model.

1 2 3 4 5 6 A B C [ ] [ ] [ ] [ ] [ ]
The shear contact force is computed in an incremental fashion. When the contact is formed, the total shear contact force is initialized to zero. Each subsequent relative shear-displacement increment results in an increment of elastic shear force that is added to the current value. The motion of the contact must be considered during this procedure. Motion of the contact is accounted for by updating and [ ] every time-step. is updated by calculating two rotations: the first is about the line common to the old and new contact planes; and the second is about the new normal direction. The equations assume that the rotations are small. The first rotation is expressed by

{ } . = - 3-14
where [ ] is the old unit normal to the contact plane, is the Kronecker delta and (or ) is the permutation symbol that takes the value 0 if two indices coincide, = +1 if i,j and k permutate like 1, 2 and 3 and = -1 otherwise. The second rotation is expressed by

{ } . = . - 〈 〉∆ 3-15
where 〈 〉 is the average angular velocity of the two contacting entities about the new normal direction

〈 〉 = 1 2 + 3-16
where , are the rotational velocities of the entities A and B. The relative motion at the contact, or the contact velocity (which is defined as the veloicity of ball B relative to ball A at the contact point) is given by:

= ( ̇ ) -( ̇ ) 3-17 = ̇ | | + | | - | | -̇ | | + | | - | |
where ̇ | | and ̇ | | are the translational velocities of the entities A and B. The contact velocity can be resolved into normal and shear components with respect to the contact plane. Denoting these components by and for the normal and shear components, respectively, the shear component of the contact velocity can be written as

= - 3-18 = -
The shear component of the contact displacement-increment vector occurring over a time-step of ∆ is calculated by

∆ = ∆ 3-19
and is used to calculate the shear elastic force-increment vector,

∆ = -∆ 3-20
where is the shear stiffness at the contact. The new shear contact force is found by summing the old shear force vector at the start of the time-step (after it has been rotated to account for the motion of the contact plane) with the shear elastic force-increment vector:

= { } . + ∆ 3-21
The values of the normal and shear contact forces determined are adjusted to satisfy the contact constitutive relation. The constitutive behavior of a material is simulated by associating a contact model with each contact. In addition to the contact model, there may also be a bond and a dashpot. These three components, define the contact force-displacement behavior.

Slip behavior is provided by enforcing a relation between shear and normal force, such that the two contacting entities may slip relative to one another. The relation provides no normal strength in tension, and allows slip to occur by limiting the shear force. Slip behavior is always active, unless a contact bond is present, in which case the contact-bond behavior supersedes the slip behavior.

The slip behavior is defined by the friction coefficient at the contact μ [dimensionless], where μ is taken to be the minimum friction coefficient of the two contacting entities. The criterion of nonormal strength is enforced by checking whether the overlap is less than or equal to zero. If it is, then both the normal and shear contact forces are set to zero. The contact is checked for slip conditions by calculating the maximum allowable shear contact force:

= | | 3-22 If >
, then slip is allowed to occur (during the next calculation cycle) by setting the magnitude of equal to = .

Bond model and fracture bond model

It is possible to bond particles together at contacts. The two standard bonding behaviors are embodied in contact bonds and parallel bonds. The contact-bond glue is of a vanishingly small size that acts only at the contact point, while the parallel-bond glue is of a finite size that acts over a circular cross-section lying between the particles. The contact bond can transmit only a force, while the parallel bond can transmit both a force and a moment. Both types of bonds may be active at the same time. However, the presence of a contact bond preclude the slip behavior, not so the parallel bond. Once a bond is formed at a contact between two particles, the contact continues to exist until the bond is broken. Particles may be bonded only to one another; a particle may not be bonded to a wall. Both types of bonds are also created between all proximate particles within the given range, regardless of whether a normal force exists between them. Both types of bonds are deleted if their strengths are exceeded.

Contact Bond

A contact bond can be considered as a pair of elastic springs (or a point of glue) with constant normal and shear stiffnesses acting at the contact point. These two springs have specified shear and tensile normal strengths. The existence of a contact bond precludes the possibility of slip (i.e., the magnitude of the shear contact force is not adjusted to the allowable maximum). Instead, the magnitude of the shear contact force is limited by the shear contact bond strength. Contact bonds also allow tensile forces to develop at a contact. These forces arise when there is no overlap between balls. In this case, the contact bond acts to bind the balls together.

The magnitude of the tensile normal contact force is limited by the normal contact bond strength. A contact bond is defined by the following two parameters: normal contact bond strength, [force]; and shear contact bond strength,

[force]. If the magnitude of the tensile normal contact force equals or exceeds the normal contact bond strength, the bond breaks, and both the normal and shear contact forces are set to zero. If the magnitude of the shear contact force equals or exceeds the shear contact bond strength, the bond breaks, but the contact forces are not altered, provided that the shear force does not exceed the friction limit, and provided that the normal force is compressive. The force-displacement behavior relating the normal and shear components of contact force and relative displacement for particle contact occurring at a point is shown in Figure 3-3. At any given time, either the contact-bond or the slip behavior is active. In this figure ,  is the normal contact force, where: > 0 indicates tension; is the relative normal displacement, > 0 indicates overlap;

is the magnitude of the total shear contact force; and is the magnitude of the total shear displacement measured relative to the location of the contact point when the contact bond was formed.

Parallel Bond

A parallel bond provides the force-displacement behavior of a finite-sized piece of cementatious material deposited between two balls that are treated as spheres. These bonds establish an elastic interaction between particles that acts in parallel with the slip or contact-bond behaviors described above. Thus, the existence of a parallel bond does not preclude the possibility of slip. Parallel bonds can transmit both forces and moments between particles, while contact bonds can only transmit forces acting at the contact point. Thus, parallel bonds may contribute to the resultant force and moment acting on the two bonded particles. A parallel bond can be considered as an elastic beam of circular cross section. This elastic beam acts in parallel with the point-contact springs that are used to model particle stiffness at a point, and whose force-displacement behavior is shown in Figure 3-3. Relative motion at the contact (occurring after the parallel bond has been created) causes a force and a moment to develop within the bond material as a result of the parallel-bond stiffness. These force and moment act on the two bonded particles, and can be related to maximum normal and shear stresses acting within the bond material at the bond periphery. If either of these maximum stresses exceeds its corresponding bond strength, the parallel bond breaks. , with the convention that this force and moment represent the action of the bond on sphere B of Figure 34. Each of these vectors can be resolved into normal and shear components with respect to the contact plane as:

= + 3-23 = +
where , , and are plotted in Figure 3-4. These vectors are shown in Figure 34, where the parallel bond is depicted as a beam of elastic material. When the bond is formed, and are initialized to zero. Each subsequent relative displacement and rotation-increment at the contact results in an increment of elastic force and moment that is added to the current values. The elastic force-increments occurring over a time-step Δt are calculated by

∆ = - 3-24 ∆ = - ∆ with ∆ = ∆
and the elastic moment-increments by

∆ = ∆ 3-25 ∆ = -∆ with ∆ = - ∆
where the contact velocity is given by Eq. 3-18, A is the area of the circular cross section of the bond beam, J is the polar moment of inertia of the circular cross-section, and I is the moment of inertia of the circular cross-section about an axis through the contact point. These quantities are given by: A

B R A R B = 3-26 = 1 2 = 1 4
where is the radius of the parallel bond. The new force and moment vectors associated with the parallel bond are found by summing the old values existing at the start of the time-step (after the shear component vectors have been rotated to account for the motion of the contact plane) with the elastic force-and moment-increment vectors. The maximum tensile and shear stresses acting on the bond periphery are calculated (via beam theory) to be

= - + | | 3-27 = | | + | |
where the values of A, I and J are given by Eq. 3-26. If the maximum stresses exceed the normal and shear strengths respectively or a combination of these maximum stresses exceeds an equivalent stress, then the parallel bond breaks.

Energy dissipation

Energy supplied to the particle system is dissipated through frictional sliding. However, frictional sliding may not be active in a given model or, even if active, may not be sufficient to arrive at a steady state solution in a reasonable number of cycles. So a so-called local damping and viscous damping are available to dissipate kinetic energy. Local damping acts on each ball, while viscous damping acts at each contact. Local damping applies a damping force, with magnitude proportional to the unbalanced force, to each ball. Viscous damping adds normal and shear dashpots at each contact. These dashpots act in parallel with the existing contact model, and provide forces that are proportional to the relative velocity difference between the two contacting entities (ball-ball or ballwall). For compact assemblies, local damping is the most appropriate form to establish equilibrium and to conduct quasi-static deformation simulations. When a dynamic simulation of compact assemblies is required, the local damping coefficient should be set to a low value appropriate to energy dissipation of dynamic waves. Alternatively, viscous contact damping should be used. For problems involving free flight of particles and/or impacts between particles, local damping is inappropriate, and viscous contact damping should be used. In the present work, due to the fact that simulations will be performed for compact assemblies and for quasi-static configuration set up, only local damping will be used.

Local Damping

The local damping used is similar to that described by [START_REF] Potyondy | A bonded-particle model for rock[END_REF]. A damping-force term is added to the equations of motion, given by Eq. 3-1 and Eq. 3-4, such that the damped equations of motion can be written as:

+ = ; = 1 … 6 3-28 = ̈ ; = 1 … 3 ̇( ) ; = 4 … .6
where , , and are the generalized force, mass, and acceleration components, respectively; includes the contribution from the gravity force; and is the damping force:

= -| | ( ); = 1 … 6 3-29
The damping force is controlled by the damping constant α and can be specified separately for each particle. In the DEM literature the value of α = 0.7 is commonly used and in the present work is going to be implemented for all the particles. This form of damping has the following advantages:  Only accelerating motion is damped; therefore, no erroneous damping forces arise from steady-state motion.  The damping constant, α, is non-dimensional.  Since damping is frequency-independent, regions of the assembly with different natural periods are damped equally, using the same damping constant.

Critical time step, stiffness of the model and simulation time

The equations of motion are integrated using a centered finite-difference scheme expressed by Eq. 3-6 and Eq. 3-7. The computed solution produced by these equations will remain stable only if the time-step does not exceed a critical time-step that is related to the minimum eigen-period of the total system (PFC-3D, Itasca). However, global eigen-value analyses are impractical to apply to the large and constantly changing systems typically encountered in a DEM simulation. Therefore, a simplified procedure is implemented to estimate the critical time-step at the start of each cycle. The actual time-step used in any cycle is taken as a fraction of this estimated critical value.

Estimation procedure for solution stability

First, consider a one-dimensional mass-spring system described by a point mass, m, and spring stiffness, k, with the coordinate system shown in Figure 3-5a. The motion of the point mass is governed by the differential equation -= ̈ . The critical time-step corresponding to a secondorder finite difference scheme for this equation is given by:

= ; = 2 3-30
where T is the period of the system. Next, consider the infinite series of point masses and springs in Figure 3-5b. The smallest period of this system will occur when the masses are moving in synchronized opposing motion such that there is no motion at the center of each spring. The motion of a single point mass can be described by the two equivalent systems shown in Figure 3-5(c) and (d). The critical time-step for this system is found, using Eq. 3-30, to be:

= 2 (4 ) ⁄ = ⁄ 3-31
where k is the stiffness of each of the springs in Figure 3-5b. The preceding two systems characterize translational motion. Rotational motion is characterized by the same two systems, in which mass, m, is replaced by moment of inertia I , of a finite-sized particle, and the stiffness is replaced by the rotational stiffness. Thus, the critical time-step for the generalized multiple mass-spring system can be expressed as

= ⁄ , ( )
⁄ , ( )
3-32
where and are the translational and rotational stiffnesses, respectively, and is the moment of inertia of the particle. The system modeled in DEM is a three-dimensional collection of discrete bodies (spherical particles and clumps) and springs. Each body may have a different mass, and each spring may have a different stiffness. A critical time-step is found for each body by applying Eq. 3-32 separately to each degree of freedom, and assuming that the degrees of freedom are uncoupled. The stiffnesses are estimated by summing the contribution from all contacts, using only the diagonal terms of the contact stiffness matrix. The final critical time-step is taken to be the minimum of all critical time-steps computed for all degrees of freedom of all bodies.

Porosity and inter-particle forces

The porosity is defined as the ratio between the volume of voids divided by the total volume occupied by the material :

= = 1 - 3-33
stands for the volume occupied by the solid material. In the present case focused on natural rock, porosity can vary from 0 to as much as 90% depending on the kind of rock and its deep in the earth crust [START_REF] Goodman | Introduction to rock mechanics[END_REF]. Porosity is known to have a great influence on mechanical and hydraulical behavior. When modeling rock-like materials with DEM, the usual objective is to create an irregular or regular packing of spherical elements to fill some given space bounded by walls at a given porosity, and to ensure that the assembly is at equilibrium. It is important to note that if the term porosity is used in the following, this porosity must be considered as a numerical one, i.e., a porosity linked to the discrete character of the model which is not necessarily comparable to the experimental porosity of the modeled material. There are limits to the porosity that can be achieved; clearly it is impossible to pack spherical particles to an arbitrarily low porosity. At some values of porosity, the particles can be arranged without touching; at other values, there will be overlap between particles. In extreme cases of DEM models with very low porosity, the inter-particle forces can be so strong that could influence the fracture behavior of the model in the sense that, when the particles are bonded together and when a bond breaks, the elastic energy stored in the bond and then transformed in kinetic energy of the particles, can be high enough to produce a large avalanche of fractured bonds, an so, fracturing all the model in unrealistic manner.

Representative Elementary Volume (REV) in DEM

In order to obtain reliable results in a shear joint simulation, the DEM model must describe accurately the elastic behavior of rock (Young modulus and Poisson's ratio) as well as the fracture behavior and especially, the model must exhibit a ratio approximately equal to 10 between the compression and tension strength as experimentally has been observed in a quasibrittle material. In the present section, the procedure used to obtain the expected elastic properties of rock is described and a fracture criterion of bond model allowing to translate the fracture behavior of rock with a reasonable accuracy is proposed.

One of the major issues in DEM is to determine the minimum number of particles needed to correctly simulate the mechanical behavior of an isotropic material. The estimation of the minimum number of particles leads to the determination of the representative elementary volume (REV) of the material formed by discrete particles. [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF] have studied, in a systematic way, this problem. The main idea argued by the author is that the geometrical isotropy is a necessary condition to ensure the mechanical isotropy of the simulated material. According to [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF] an isotropic DEM model, is obtained if the radius dispersion of the arrangement of spherical particles is higher than 15% and the number of discrete elements is around 10 000. The radius dispersion can be expressed as:

= - 〈 〉 * 100% 3-34
, , 〈 〉 stands for minimum, maximum and average radius.

As previously mentioned, the choice is made to simulate the behavior of mortar replicas of granite samples created by [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF] and studied by [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF]. In this thesis, the spherical particles are chosen in order to mimic the granulometry of the mortar replicas used in [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF]. The radii of the particles are taken from a Gaussian distribution Figure 3-6 with mean 〈 〉 of 0.5 mm and standard deviation of 0.1 mm giving a radius dispersion > 15 % (the actual values of the minimum and maximum particle radius in the REV are = 0.19 , = 1.00

). The set of spherical particles are generated inside a cubic volume of 1 cm by side confined by boundary rigid walls as shown in Figure 34567. Initially the population of particles is created with artificially small radii at random coordinates and then expanded until the final desired radius is achieved and the mean porosity 〈 〉 is obtained. Note that, only the granular skeleton of mortar (sand) is modeled by the discrete elements, the cement matrix is not physically represented but modeled by the parallel bonds. Thus the REV (DEM model) presents a numerical porosity which does not have to be directly compared to the experimental porosity of mortar. As explained before (section 3.2.6), the porosity can influence the fracture behavior of the DEM. To avoid this kind of problems and with the aim of reproducing a coherent fracture behavior, a mean porosity 〈 〉 of 0.51 is used in this study. This value was chosen after trying various levels of the porosity and regarding the computation time and the dimensional stability of the final DEM model. Indeed, the numerical porosity (0.51) selected in this study results from 1) a compromise between the number of particles in the REV (and so, in the final joint model) conditioning the computation time and, 2) the internal stress level which controls the dimensional stability and fracture behavior of the DEM model (rough surface of the DEM joint). The mean porosity is defined as:

〈 〉 = 1 - 〈 〉 3-35
where 〈 〉 is the average volume of a spherical particle, , is the number of discrete elements in the REV and is the volume of the REV corresponding to 1 cm 3 cube bounded by walls. Substituting the values of the mean porosity and the average volume of the particles in Eq. 3-35, the number of discrete elements needed to generate the REV of 1 cm by side is ≈ 935. Note that this number is almost ten times lower than that proposed by [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF], nevertheless, as will see later in this chapter, this number of particles is sufficient enough to reproduce the mechanical and fracture behavior of a typical rock with reasonable accuracy. In summary, in the discrete element model, the spherical particles represent the granular skeleton (sand) of mortar while the cohesion linked to the cement matrix is performed from a lattice of beams (parallel bonds) and as a consequence, the void space of the DEM have not be compared to the void space of the actual material (mortar) as previously mentioned.

Once the system of particles in the REV is in mechanical equilibrium, they are bonded together and the boundary walls are removed, then, the system is cycled until the new mechanical equilibrium is reached. The contact bond and the parallel bond (see section 2.2.3 of this chapter) were used simultaneously in order to provide a symmetric behavior both in tension and compression.

Let us now detail the procedure used to obtain the elastic properties expected for mortar reported in [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF] after [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF]. A rock like granite is formed by grains with complex shapes that are bonded together. The geometry of the grains and the interlocking between them make their rotations difficult and, when a rock sample is loaded in quasi-static condition, the angular velocity of the grains could be very small. As noted by [START_REF] Cho | A clumped particle model for rock[END_REF], when spherical particles are used to model the microstructure of natural rocks, the interlocking of the natural grains cannot be captured correctly and the rotation of the discrete spherical elements become an important factor controlling the strength of the DEM model, particularly the unconfined compression strength. In the present work, the interlocking between grains in a natural rock is modeled by restricting the particle rotation when they are bonded together and allowing them to rotate freely when they are not bonded. However, a Coulomb's friction coefficient between particles of 0.5 (typical value used in the DEM literature) was used to represent the micro-roughness of the natural rocks grains.

It is important to remark that, in this study, there are no intentions to create a DEM model with the stiffness of a natural rock, which are in the order of tens GPa for the elastic modulus E. Trying to make a DEM with such elastic modulus, implies assigning large values for the contact stiffness between discrete elements. As previously mentioned in section 1.2.5 of this chapter, the critical time step during the simulations is inversely proportional to the contact stiffness, and assigning large values for the stiffness can produce very small and, hence, an excesive computation time.

The solution proposed in this study consists in considering a model sufficiently soft in such a way that the corresponding critical time step permits to carry out shear joint simulations in a couple of days even for very large systems of particles (the numerical rock joints generated and presented after in this thesis, are formed by around 100 000 discrete elements). A value of 0.1 N/m is assigned to the normal contact stiffness between particles and the shear contact stiffness is fixed ten times lower than in order to get an acceptable Poisson's ratio of the REV. To compare the results of the simulations with experimental data (mortar replicas of a natural rock joint created and studied by [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF], the external forces/stresses applied to the model will be of such magnitude that they will produce the same deformation in the REV as the experimental deformation. To clarify this point let us consider the example of a natural rock with elastic modulus submitted to an uniaxial compression stress . The normal stress that must by applied to the REV, of elastic modulus in order to produce a strain equivalent to the experimental one is:

= 3-37 where = =
The Young's modulus , shear modulus and Poisson's ratio of the REV, are estimated from the simulation of an uni-axial compression and a direct shear test in quasi-static conditions. A local damping coefficient = 0.7 was used to dissipate enough kinetic energy and ensure stable conditions during the simulations, (see section 2.2.4). To apply the normal and shear displacements in the simulations, an ensemble of discrete elements were used to form the clumps shown in Figure 3-8, these clumps acts as rigid bodies with deformable boundaries. During the compression-tension tests, the upper clump of particles was displaced vertically by moving it at constant velocity , = 1 μ / , while the lower clump was fixed. Both clumps are not allowed to rotate. Every timestep, the vertical force applied by the particles in contact with the clumps is estimated. The axial and lateral deformation of the REV are computed by applying equations 3-38 to 3-41 taking into account the particles located in the middle of the REV and labeled by a set number as shown in Figure 3

-8. = 〈 〉 -〈 〉 2 ⁄ 3-38 = 〈 〉 -〈 〉 3-39 = 〈 〉 -〈 〉 3-40 = + 2 3-41
〈 〉, 〈 〉, 〈 〉 stand for the average value of the x,y,z coordinates of the particles at the time t and belonging to the set i, L is the size of the REV, and , , are the strains in the x, y and z directions.

In the direct shear test simulation, the upper clump was displaced horizontally at constant velocity ( ) = 1μ / in the x direction, every time step the unbalanced shear force in the clumps were measured and the shear deformation computed by using Eq. 3-42: which is in agreement with the one of the shear modulus obtained from the shear simulation (the discrepancy between both values is less than 4%).

= 〈 〉 -〈 〉 2 ⁄ 3-42

Fracture behavior of REV

The rock strength is an important factor to consider in many engineering projects like the construction of dams, mines, radioactive wastes storages, building construction, and gas or petroleum extraction. The usual procedures to estimate the compression and tensile strength of the rock at the laboratory scale consist of performing uni-axial compression test and brazilian test respectively.

The aim of this section is to propose a fracture criterion for the bonds in the DEM-REV model capable of reproducing the typical fracture behavior of rocks under uni-axial compression and tension tests. The elastic properties of the REV were tested in the previous section, so the next step consists in choosing a fracture criterion for the bonds and calibrate it with respect to the experimental data shown in Table 3-1. This data were extracted from the thesis of [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF] who studied mortar replicas of granite created by [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF]. The fracture criterion used for the bonds is given by Eq. 3-46, which represents an ellipse centered at (ℎ,0) and with major and minor axis a and b respectively.

-ℎ + ≥ 1 3-46 and are the maximum normal stress and the maximum shear stress in the bonds. These stresses are computed using Eq. 3-27 based on beam theory as explained previously in this chapter. If the normal and shear stresses in a bond make the left hand side of Eq. 3-46 be greater or equal than 1, the bond will break irreversibly. The parameters a, b and h given by equations 3-47 to 3-49 are chosen in such a way that the shear strength of the bonds in pure shear is five times bigger (parameter = 5 in Eq. 3-49) than the tension strength of the bond in pure tension =

5

, and the ratio between compression strength and tension strength of the bond is 10.

ℎ = - 9 2 3-48 = * * ( -ℎ ) 3-49
Note that the fracture criteria for the bonds given by Eq. 3-46 has only one independent parameter, the tension strength of the bond = 0.025 , this parameter is used as a calibration variable and its value is chosen such that the tension strength of the REV is equivalent to the experimental tension strength given in Table 3-1. However, cannot be chosen in an arbitrary manner, will depend on the distribution of the internal stresses in the REV. As mentioned before, once the boundary walls used to create the REV are suppressed, the bonds between the particles will deform and develop stresses until the mechanical equilibrium of the system of particles is reached. Figure 3-10 a and b shows an example of such bond stresses when the REV reaches the equilibrium just after the boundary walls were removed. It is intuitive that cannot take any arbitrary value. Indeed, if is less that the maximum tensile stress in the REV, some of the bonds will break instantaneously before carrying the compression and tension tests simulations. So in this work, take values lying at the right side of the normal stresses distribution shown in Figure 3-10a such that the compression stress as well as = 11 2 3-47 the tensile strength obtained from the REV correspond approximately to the strength expected for mortar. Note that the internal stresses (in the parallel bonds) give an internal disorder to the REV. Note that the values of the compression and the tension strength obtained from the REV are in agreement with the experimental ones reported in Table 3-1 (74.6 MPa and 6.6 Mpa, respectively).

Also note that the ratio between compression strength and the tensile strength of the REV is around 10 as expected in natural rocks. On the other hand, Figure 3-12 shows the broken bonds (reflecting the micro-cracks expected in mortar) at the end of the tension (Figure 3-12a) and compression (Figure 3-12b) simulations. Note that, in the tension case the broken bonds are located within a more or less narrow region oriented perpendicular to the tensile loading axis. This result is in agreement with the preferential orientation of the micro-cracking in mortar well known to be perpendicular to the axis of loading during a tension test. In compression (Figure 3-12b), the microcracks spread in the whole volume without apparent preferential orientation. Note that the diffuse location of the broken bonds in this case involves some difficulties to observe accurately the eventual bands location (oriented perpendicular to the maximum direction of the extension of the material).

One can note, despite that the relative simplicity of the fracture criterion defined by Eq. 3-46, the fracture behavior obtained in tension as well as in compression exhibits a reasonable agreement with the quasi-brittle fracture behavior expected for mortar.

(a) (b)

Figure 3-12 Cracking of the REV at the end of (a) tension and (b) compression simulations.

Finally, Figure 3-13 shows the intact bonds before and after tension/ compression simulations. Note that the density of intact bonds is higher at the end of the tensions test than the compression test.

Generation of self-affine rock joints with DEM

Once the elastic and fracture properties of the REV are calibrated and compared to experimental data (mortar replicas of a rock joint created [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF]), the roughness of the same mortar replicas [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF] after [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF]) is used to generate the joints using DEM, particularly, the roughness of profiles obtained in the b direction of the mortar replicas, Table 3-2a. (see [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF] for more details). This section presents a methodology to generate such rough joints which will be used to carry out the compression and shear test simulations in the next chapter. Thus, to produce the rough joints with DEM the numerical generator of self-affine surfaces is used to create a set of square self-affine surfaces of 10 cm by side. A total of eight different surfaces were generated

X Y Z X Y Z X Y Z X Y Z X Y Z X Y Z X Y Z
for the same phase distribution of the surface S1, and corresponding to two values of the roughness exponent, the correlation length and the height variance as detailed in Table 3-2. The synthetic surfaces are isotropic in roughness exponent H, correlation length L c and height variance σ². Moreover, in order to compare the results of compression and shear tests simulations with experimental data, the morphology of the surface S5 was generated with the same height variance, correlation length and roughness exponent of the mortar replicas created by [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF] and studied by [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF]. The other values for H, Lc and σ² of the synthetic surfaces in Table 3-2 where chosen arbitrarily, however, these values lies on the range of values founded experimentally (0.4 to 0.8 for H, and 1.5 cm to 20 cm for L c ). The value of σ² L =0.328 mm² was select also arbitrarily to facilitate the creation of the numerical DEM joints. Figure 3 Note that, as shown in Figure 3-15, the fact of using the same phase distribution for the surface generation allows to obtain fixed locations of the hills and holes in the average surface plane. Also note that the synthetic surfaces S4 and S5 are the smoothest (High H, High L c and Low σ²: H H -L c-H -σ² L ) and roughest surfaces (Low H, Low L c and High σ²: H L -L c-L -σ² H ). Surfaces S1 to S8 are used as a "numerical mold" which serves to form the rough joints with discrete spherical elements in DEM. The upper and lower part of the joint are created independently and put it together a posteriori. The procedure to generate the rough surfaces in DEM consists in the following steps:

1. Creation of a mesh composed of triangular walls using the coordinates x,y,z of the synthetic surfaces S1-S8 (Figure 3-17a). 2. A set of spherical particles is generated in the region enclosed by the boundary walls and the numerical mold of the self-affine joint (Figure 3-17b). The procedure to generate the particles was described before in this chapter and the micro-mechanical properties of the particles are the same of those obtained for the calibrated REV. The walls of the numerical

Morphology of the rough surface formed from discrete elements

It is expected that the surface formed from the discrete elements will not match exactly the synthetic surfaces used as molds because the discrete character of the DEM and its intrinsic porosity.

Characterization of the surfaces generated with DEM.

In order to compare the surfaces formed from the discrete elements to the synthetic ones (molds) it is necessary to extract the heights (z coordinate) of the DEM surface at discrete points (x,y) that coincides with the (x,y) coordinates of the meshes of the synthetic surfaces used as molds. Thus, a "numerical profilometer", similar to the one used in The mesh used to obtain the DEM-joint roughness is formed by 365x365 points, the distance between points is Δx = Δy = 0.308 mm, 3.24 times larger than the mean diameter of the particles.

Comparison of Figure 3-18e and Figure 3-18f allows to visualize the effect of the discretization of the surface linked to the DEM. Several points of the Figure 3-18f exhibit height values greater than the heights of the corresponding points of the synthetic surface: these points correspond to spherical elements having penetrated the walls of the mold. Moreover, numerous points exhibit height values lower than the values of the corresponding points of the synthetic surface. These points exhibit the intrinsic porosity of the DEM and increase the tortuosity of the discrete surface. Thus, the problem linked to the discretization boils down to estimating the points of the discrete surface in agreement with the synthetic one and, as a consequence, those which exhibit height values out of the synthetic one and especially those linked to the porosity.

Modified Abbot's curve to find the mean plane of the DEM surface.

Once extracted the rough discrete surface from the DEM, the next step consists in estimating the average mean planes of DEM-joints, in order to compare it a posteriori with their corresponding ideal morphologies S1-S8. For this, a modified Abbot's curve procedure is implemented. The first step of the method consists in aligning the corresponding DEM and synthetic surfaces in the x-y plane. Then, in a second step, the synthetic surfaces is offset by a vertical distance h from its original mean plane (z-coordinate) while the DEM surfaces remains at a fixed position. Then, the third step consists in displacing vertically the synthetic surface towards the DEM surface and counting the number of contact points between both surfaces as a function of the relative vertical distance dv (Figure 3-19a).

In the ideal case (that means the DEM and synthetic surfaces are identical), when the mean plane of both surface coincides (dv = 0), all the points are in contact and a plot of number of contact vs dv must be a vertical line as the red-dotted-line in Figure 3-19c, note that the slope of this curve at dv = 0 corresponds to the highest value of 90°. When the DEM and synthetic surface are not identical, the plot of number of contacts vs dv exhibits the shape of the black-solid-line of Figure 3-19c. In the present case, the average mean plane of the DEM surface (dv mean-plane ) is considered to be as the one for which the slope of the number of contact points vs dv curve reaches the maximum value (green line in Figure 3 Thus, due to the discrete character of the DEM and its high porosity, only 40% of the points of the DEM joint-surface coincide with the original synthetic surfaces S1-S8. So we can consider that the DEM-joints are formed by the synthetic surfaces S1-S8 at which a noise in the z-coordinate is added. The main question is then: what is the maximum noise level for which the morphological characteristics of the surface (roughness exponent, correlation length and height variance) are retained?. To answer this question, let us analyse how the addition of an uncorrelated random elevation (i.e. noise) affects the scaling properties of the synthetic surface. For this, additional uncorrelated elevation maps are taken as a fraction of the height of the synthetic surfaces ranging from 10 to 100 % with a step of 10% and added to the z coordinates of the corresponding synthetic surfaces Figure 3-21a show how the addition of a noise affects the scaling properties of the synthetic surfaces. From Figure 3-21a it is possible to note that the noise affects essentially the roughness exponent and the height variance. Nevertheless, as long as the magnitude of the noise is lower than 30%, the roughness exponent and the height variance are little changed. As a consequence, the choice is made to consider as representative of the expected rough surface, only the points of the discrete surface for which the elevation is ranged between the elevation of the corresponding point of the synthetic surface considered at the vertical position dv mean-plane plus or minus three times the height variance of the synthetic surface, i.e.,. Note that in Figure 3-21a there is a cutoff of the self-affine regime around 1000 µm. This lower limit of the scaling behavior coincides with the average diameter of the discrete elements <D> = 1000 µm. Note also that, from the previous procedure used to identify the spheres forming the self-affine surface mentioned above, the DEM joint and the synthetic joint has the same variance and correlation length. However, the roughness exponent of the DEM joint is lower than the roughness exponent of the original surfaces used as a mold, but coincides with the roughness exponent of the original synthetic surfaces when a noise of 20% is added. Figure 3-23 shows of the variograms computed for all the DEM joint surfaces generated using the synthetic surfaces S1-S8 as molds and Figure 3-22 shows their corresponding elevations maps. Table 3-2b shows the roughness exponent H, correlation length Lc, and height variance σ² of the self-affine joints created with DEM. 3-2b that the roughness exponent of the DEM surfaces have a lower roughness exponents H that the corresponding surface "molds". However, the results of the compression and shear simulations that will be obtained for the joint DEM S1 could be compared with the

Note from Table

S1 S2 S3 S4 S5 S6 S7 S8 L C-H L C-L 2σ 2 L 2σ 2 H H H H L
experimental ones obtained for the shear direction b, [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF], because their roughness parameters are similar (Marache, 2000), Table 3-2c. Remark that the particles identified as being out of the expected rough surface, i.e., the yellow particles in Figure 3-21c, correspond mainly to particles located in the vicinity of the sharp holes linked to the intrinsic porosity of the DEM.

Thus, during the simulation of the mechanical tests, if those particles located out of the expected rough surface are implicated in the friction phenomenon, the contribution of these particles in the whole friction forces will be estimated (these contributions will be termed "artifact shear stress"). If this contribution remains weak, the simulation of the mechanical test will be considered as representative. Otherwise, the result will not be retained.

DEM simulations of direct shear test of rock joints under constant normal load (CNL)

Overview

This chapter proposes DEM simulations of direct shear test of rock joints under constant normal load.

The rough joints generated in Chapter 3 are successively compressed then sheared and this considering the different material behaviors: rigid, elastic and fracture-elastic.

Once calibrated the REV model (in Chapter 3) and in conjunction of the numerical generator of selfaffine surfaces explained in Chapter 2 of this work, let us use these results to generate several DEM rough joints surfaces and then compress them and sheared them to study the influence of the joint roughness, elasticity and fracture on their mechanical behavior.

Several authors have used numerical models to study the mechanical, thermal, and chemical behavior of intact rock and rock joints. For example, [START_REF] Karami | Asperity degradation and damage in the direct shear test: a hybrid FEM=DEM approach[END_REF] studied joint surface damage using a 2D hybrid FEM/DEM code. They use various joint profiles with different geometrical configuration to simulate quasi-static direct shear test under various normal loads. The roughness of profiles was selected using Barton's classification (JRC). One of the main results of their numerical simulations shows that dilation is controlled mainly by the joint surface roughness and the applied normal load. When the applied normal load is low, significant dilation occurs. In contrast, when the applied normal load is high, dilation is low and damage is composed of both surface wear and asperity breakage. The same year, [START_REF] Cho | Development of a shear one in brittle rock subjected to direct shear[END_REF] carried out experimental and DEM 2D-simulation of direct shear test at various constant normal loads for an intact brittle synthetic rock. Results of their DEM simulations show that a macro shear-band is formed as a consequence of progressive failure of tension-induced micro-cracks. The next year, Park and Song (2009) made 3D-DEM simulations of direct shear test for rock joints. They pay special attention to the influence of particle size, ball friction coefficient, joint roughness (based in the Joint Roughness Coefficient JRC defined by Barton) and joint contact bond strength in the shear response of the model. They showed that ballfriction coefficient was the most important factor governing the shear strength and dilation angle of the joints, on the other hand, they found that variation in joint roughness and contact bond strength had a larger effect on the cohesion and peak friction angle. Later, [START_REF] Lambert | Influence of calcium leaching on the mechanical behavior of a rock-mortar interface: A DEM analysis[END_REF] studied the influence of leaching on the mechanical behavior of rock (granite)-mortar interfaces by means of DEM simulations. They carried out direct shear stress simulations under constant normal load and several leaching times were used to simulate the increase of the macro-porosity resulting from the leaching. They found the existence of a transitional leaching depth, at which the joint exhibits neither dilation nor contraction. Moreover, numerical results exhibited a decrease of Young's modulus and a decrease of UCS with leaching time. [START_REF] Duriez | A discrete modeling-based constitutive relation for infilled rock joints[END_REF] used a 3D-DEM model to calibrate an incrementally nonlinear constitutive relation for in-filled rock joints that can reproduce rock joint features, i.e. the dilatancy. The same year, [START_REF] Renouf | Multiphysical modeling of third-body rheology[END_REF], used DEM models to simulate, mechanical, thermal and chemical behavior of third-body systems in quasi-stable and dynamic conditions. Finally, more recently, [START_REF] Bahaaddini | Numerical direct shear tests to model the shear behavior of rock joints[END_REF], made direct shear test simulations of 2D saw-tooth triangular joints. They showed that the smooth joint model used commonly in DEM, exhibits particle interlocking which takes place at shear displacements greater than the minimum diameter of the particles and that, this particle interlocking affects the mechanical behavior of the joints. In order to solve this problem, they proposed a new methodology consisting of generating separately the upper and lower parts of the rock joint and then applying the smooth joint model to the contacts formed between the two blocks.

Nevertheless, in spite of various numerical studies focused on shear test of rock joints, the knowledge of the influence of the roughness on the mechanical behavior of a joint is still rather imprecise. Indeed, to our knowledge, a 3D DEM simulation of a shear test using 3D realistic joint roughness has never been performed. This chapter will attempt to provide contributions to this scenario.

In the domain of rock mechanics, the loads applied the rocks are essentially promoted by gravity. In cases where the rocks are found deep in the earth crust, the loads applied by the upper rocks can be so high that they inhibit the dilation of the rock joints when they displace relative to each other. It is common to model the shear behavior of rock joints in such circumstances as corresponding to a constant normal stiffness. When the studied rock joints are at shallow depths up to 3 to 4 km in the earth crust, dilation of the rock joints can occur. Then, a common model for rock joins (located at moderate depth) consists of shearing the joint under a constant normal load. In the present study, the second case is analyzed: the shear behavior of rock joints is studied under constant normal load.

In this section, the compression (closure) and direct shear test simulations at constant normal load are presented.

Uniaxial Compression Test (UCT) simulation of rock joints

Joint preparation

In Chapter 3 it was mentioned that the upper and lower walls forming the numerical joints were generated separately. Thus, before the realization of the compression test itself, the upper and lower part of the joint are assembled by moving vertically each other until they have only one contact point between them (Figure 4-1a).

Boundary conditions in compression

The application of forces and displacements to the joint was made from clumps. For each joint, two clumps are used, one for the lower sample and one for the upper sample. In Figure 4-2 are Upper sample

Lower sample

Contact point represented these clumps in blue. To carry out the compression or closure test, the lower clump is considered as fixed in position and its translational (V x , V y , V z ) and angular velocities (ω x , ω y , ω z ) were set to zero, while the upper clump can move along the z axis according to a constant compression velocity V z = 6 x 10-5 m/s while its translational velocities Vx and Vy as well as it angular velocities ω x , ω y , ω z . were set to zero. During the closure test, the unbalance forces in the clumps were recorded, and the z-position of the upper clump was recorded every 500 cycles.

Results of the joint's closure simulations

The eight joints (DEM S1 -DEM S8) generated in the previous section, where subjected to the joint closure test. For each joint, three different mechanical behaviors were simulated: a) Rigid model, b) Elastic model and c) Elastic-fracture model.

In the rigid model, the particles belonging to the same sample (upper sample or lower sample) keep constant its relative positions during the simulation, this means that the upper and lower sample keep their initial shape during the test. Nevertheless, the surface particles of the upper sample can interact mechanically with the surface particles of the lower sample. So the rigid model is deformed only at the interface between the lower and upper samples.

In the elastic model, all the particles forming the joint sample can interact with its neighbors and change its relative positions when a load is acting upon them, but the bonds between them cannot break. As a consequence, the elastic model mimics a perfectly elastic material.

Finally, in the elastic-fracture model, as in the elastic model, each particle interacts with its neighbors at their contact points, but if tensile or shear stresses in their bonds exceed certain limits, the bonds can break irreversibly. The fracture thresholds for the normal and shear stresses are the same than the ones used in the REV.

Figure 4-3 shows the results obtained for the simulations of the joint's closure test in the case of the numerical samples DEM S1 to DEM S8. Each plot shows the results obtained for the rigid, elastic and elastic-fracture model. For each joint two different compression levels were applied, 14 MPa and 21 MPa. As expected intuitively, Figure 4-3 shows that the elastic-fracture models deform (or close) more than the elastic models which, in turn, deform (or close), more than the rigid model.

Moreover, for each joint, the joint stiffness Kn (black dotted straight line in Figure 4-3) is higher for the rigid model than for the elastic and elastic-fracture models, for which the joint stiffnesses are almost the same.

Finally, Figure 4-4a exhibits typical results for the contact points between upper and lower surfaces of the joints. Note that the contact points (black points) are uniformly distributed in the whole joint surfaces. This behavior is explained because upper and lower joint surfaces are very similar (i.e., they match very well).

In Figure 4-4b, the evolution of the contact area as function of the joint closure is shown. The contact area between upper and lower surfaces of the joint is estimated according the following equation:

% ≈ ∑ × 100 4-1
where is the contact radius of the two spherical particles estimated by:

= - 4 4-2
where is the mean radius of the two spherical particles in contact, U is the overlap between particles, L correspond to the projected surface area of the joint in the plane (x, y), and the sum is made over all the contact points "cp" between upper and lower samples of the joint. The percentage are of contact estimated in the simulations is very low compared to experimental results and this is mainly because the DEM models are highly porous (i.e., porosity of 51%). Figure 4-4c shows the evolution of broken bonds during the closure test for the joint DEM S2, this result is representative for the other joints DEM S1 to DEM S8 which behave similarly.

The fracture of bonds in elastic-fracture models allows a rearrangement of particles (relative displacements between neighbors particles) which leads to an increase of the joint closure compared to the elastic model.

Direct shear test simulation under CNL of rock joints

Once the closure tests for the joints DEM S1 -DEM S8 were made, the direct shear test simulations are carried out for all joints. As in the previous section, the results of the shear test simulations are shown according to the nature of the mechanical model (rigid model, elastic model and elasticfracture model). In this section, the boundary conditions of the direct shear test at constant normal load is recalled followed by the results of the shear test simulations for all the rigid DEM joints, then the results of the elastic DEM joints are presented and discussed. Finally, the results of the shear test simulation performed according to the elastic-fracture behavior are presented. This methodology allows understanding the influence of the joint morphology, friction phenomenon, elasticity and fracturing on the shear behavior of rock joints.

Boundary conditions in shear

During the direct shear test simulations, the lower clump was kept fixed in its initial position and it is not allowed to rotate (angular velocities in x, y and z directions equal to zero). The upper clump was displaced in the x or y direction at constant velocity while a constant normal load is applied to it. During the simulation, the unbalance shear force in the upper clump was recorded as well as the x, y No. of fractures/No. surface bonds and z coordinates of its center of mass. During the test, the upper clump is free to move in the z direction. Figure 4-5 illustrates the boundary conditions applied during the direct shear test simulations.

Figure 4-5 Boundary conditions applied during the direct shear tests simulations.

Influence of joint morphology on the shear behavior

For the shear simulations, two levels of normal load were used σ 14MPa and σ 21MPa ; for each normal load and each behavior model (rigid, elastic, elastic-fracture), the joints DEM S1 to DEM S8 were sheared in the x and y directions. In total, 96 shear simulations were performed; Figure 4-6 shows the chart of shear simulations: In the next sections, the results of the direct shear test simulations are presented. Firstly, the results obtained for the rigid model are discussed followed by the results of the elastic model. Finally the results obtained for elastic-fracture model are presented and discussed.

Rigid model

Figure 4-7 and Figure 4-8 shows the results obtained for the direct shear test simulations under the compression stress levels of 14 MPa and 21 MPa respectively. The DEM S1 to S8 (see Table 3-2 in Chapter 3) are sheared along the direction x (solid line) and along the y direction (dashed line). The dilatancy responses (red curves) as well as the artifact shear stress (green curves) linked to the friction on the particles located out of the expected rough surfaces are also plotted. Note that in both figures the contribution of the particles out of the expected surface can be neglected compared with the total shear stress and hence, the simulations of shear tests can be considered as representative.

Figure 4-7 Results of the direct shear simulations for the rigid model. Surfaces DEM S1-S8 where sheared in two perpendicular directions x (continuous line) and y (broken line). The normal stress applied was 14 MPa.  The shear stress and the dilatancy are very similar for both directions of shearing in all the surfaces S1-S8, especially before the shear peak; maybe with the exception of surfaces DEM S3 and DEM S4, where the dilatancy in the y direction (broken line) is greater than in the x direction (solid line). The similar behavior in the x and y directions can be explained because the roughness of all the surfaces is isotropic [START_REF] Yang | An index for describing the anisotropy of joint surfaces[END_REF], Riss (1997[START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF], [START_REF] Park | A constitutive model for shear behavior of rock joints based on three-dimensional quantification of joint roughness[END_REF].  As the variance of the surface increases, the shear peak and the dilatancy increase; see surfaces DEM S1 compared to surface DEM S2 (surface DEM S2 has lower variance), DEM S3 vs. DEM S4, DEM S5 vs. DEM S6 (surface DEM S6 has lower variance than surface DEM S5) and DEM S7 vs. DEM S8.  As the correlation length increases, the shear peak and the dilatancy decrease; see surfaces DEM S1 vs. DEM S3 (surface DEM S3 has higher correlation length), and surface DEM S5 vs. DEM S7 (surface DEM S7 has higher correlation length than surface DEM S5).  As the roughness exponent decreases, the shear peak and the dilatancy increase , see surface DEM S1 compared to surface DEM S5 (surface S5 has lower roughness exponent than surface DEM S1), [START_REF] Xie | Fractal effects of surface roughness on the mechanical behavior of rock joints[END_REF], [START_REF] Kwafniewiski | Surface roughness evolution and mechanical behavior of rock joints under shear[END_REF]. Note that the influence of the roughness exponent on the shear peak value and the dilatancy is more pronounced than those of the correlation length and the variance.  The residual shear strength measured at 5mm of shear displacement is very similar for all the DEM S1-S8 surfaces indicating that the residual shear strength is independent of the initial surface geometry, [START_REF] Yang | Effect of asperity order on the shear response of three-dimensional joints by focusing on damage area[END_REF].  For a high value of the roughness exponent which means a low roughness of the surface (case of DEM S1, DEM S2, DEM S3 and DEM S4), the shear peak tends to disappear while for a low roughness exponent, i.e., a high roughness of the surface (case of DEM S5, DEM S6) the shear peak is more pronounced. Moreover, if a high value of the roughness exponent is combined with a high value of the self-affine correlation length and a low value of the height variance (as DEM S4: H H , Lc H , σ 2 L ) then the roughness of the surface is very low and in this case the post peak response is very flat (magnification of the effect of a high value of the roughness exponent). In the same way, if a low value of the roughness exponent is combined with a low value of the self-affine correlation length and a high value of the height variance (as DEM S5: H L , Lc L , σ 2 H ), the roughness of the surface is very strong and leads to a very pronounced peak value of the shear stress (magnification of the effect of a low value of the roughness exponent), [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF].  For a low surface roughness (DEM S4), associated to a very flat post-peak response, there is less dilatancy during shear while for a higher (DEM S5) surface roughness, associated to a very pronounced peak of shear stress-shear displacement response, the shear dilatancy is higher.  Finally note that for all the DEM S1-S8 surfaces, the shear stress exhibits a kind of periodic behavior after the shear peak is reached. It is argued that this phenomenon is related to the periodicities of the roughness of the joint's surfaces in contact. As we will see later in this chapter, these undulations in the shear stress tend to disappear when the fracture behavior is activated. The shear and dilatancy responses of the DEM S4 (very low roughness of the surface) and DEM S5 (very high roughness of the surface) appear respectively as a low boundary and high boundary of whole responses of the DEM surfaces.

Another interesting remark is that the initial shear stiffness of the joints, expressed by the slope of the shear stress vs shear displacement at the onset of the curves, is approximately the same for all the DEM surfaces, indicating that this stiffness is quasi-independent of the roughness of the joint. However, the initial shear stiffness changes with the level of compression stress (Figure 4-10), i.e., as the compression stress increases, the shear stiffness also increases, [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF]. Finally, note that there is a direct relation between the compression stress, the peak shear stress and the residual shear stress. A high compression stress leads to a high shear peak and high residual strength. The contrary is true for the dilatancy: a high compression stress leads to a low dilatancy. These behaviors are in agreement with those observed experimentally by [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF].

Figure 4-11 shows experimental results obtained from a direct shear test of a joint. The joint samples were made from a plaster material and the roughness of the joint was produced by suppressing various asperities order from an original joint surface [START_REF] Yang | Effect of asperity order on the shear response of three-dimensional joints by focusing on damage area[END_REF]. Changing the asperity order on a rough surface leads to a change of the height variance and of the correlation length of the surface. The roughness exponent of the surface is affected also but not as much as the height variance and the correlation length. Thus, Figure 4-11 exhibits the combined influence of the height variance and the correlation length on the shear response of the joint. The numerical results obtained in this thesis are in good accordance to the experimental ones mentioned above.

The dilatancy responses corresponding to the previously described shear test (Figure 4-11) are plotted in Figure 4-12a. It can be observed in that some joints contract at the onset of the test instead of dilate. The explanation of this phenomenon is multifactorial involving the degree of matching of the joint's surfaces Zhao (1997), pore closure of the sample as a function of the compression stress, and fracturing of the joint's asperities, among others reasons [START_REF] Archambault | The evolution of void spaces (permeability) in relation with rock joint shear behavior[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF]. In the present study, no appreciable contractancy was observed for all the DEM joints (Figure 4 From Figure 4-13 and Figure 4-14 it is possible to see that, at the onset of shear test (point P1) as well as the shear peak (point P2), the contact points distribution is almost homogenous in the joint plane.

As the shearing test progress, the contact points tend to localize in narrow regions that coincide with the highest zones (highest 2D colatitudes, [START_REF] Gentier | Influence of fracture geometry on shear behavior[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF]) of the joint (red zones in the digital elevation map of the surface). However, for a large shear displacement (point P4), the contact zones do not necessarily correspond to all the higher areas of the joint.

Moreover, it is interesting to note that, during the progress of the shear, the contact points tend to form clusters that are aligned perpendicularly to the shear direction, forming in this way a kind of channels in which there are not contact points. This behavior has also been observed both experimentally and numerically by several authors [START_REF] Riss | Binary images of sheared rock joints: characterization of damaged zones[END_REF], [START_REF] Nemoto | Direct measurement of contact area and stress dependence of anisotropic flow through rock fracture with heterogeneous aperture distribution[END_REF], [START_REF] Park | Numerical method for the determination of contact areas of a rock joint under normal and shear loads[END_REF]. For the rigid joint model, when the shear displacement is around 3 mm and up to the end of the numerical shear test (i.e., at 5mm of shear displacement), only a little percentage of the surfaces is in contact. These contact areas are distributed across the joint in some clusters of contact points (Figure 4-16a), and, as will see later in this chapter, these locations of contact points coincides with the zones of detached particles when using the elastic-fracture model of the joint. This behavior has been observed experimentally in several studies Riss (1997), [START_REF] Gentier | Influence of fracture geometry on shear behavior[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF].

Elastic model

Without fracture

In the previous section, the results of the shear test simulations obtained in the case of a rigid behavior of the model were analyzed. These simulations were used to highlight the influence of the surface morphology and of the friction phenomenon on the shear response of the joints. In the present section, the relative displacements between particles are allowed and as consequence, the model can represent an elastic body made of discrete elements. The elastic properties used for the simulations are the same as those estimated for the REV.

As for the rigid model, the DEM surfaces S1 to S8 are submitted to two different compression stresses (14 and 21 MPa) and sheared along the x-axis and the y-axis of the surfaces. The solid line curves correspond to the (1) shear response (black), (2) the dilatancy (red) and (3) the contribution of the particles located out of the expected rough surface on the shear stress (green) obtained for a shear movement along the x-axis while the dashed curves are associated to a shear along the y-axis (artifact shear stress).

The results of these simulations for a compression stress of 14 MPa are shown in Figure 4-17 for all the DEM surfaces S1-S8, while the results obtained for a compression stress of 21 MPa are plotted in Figure 4-18. Note that the contribution of particles out of the expected rough surface in the shear stress is completely negligible (artifact shear stress). Thus, it seems that elasticity decreases the contribution of these particles compared to the case of the rigid behavior.

a) b)

Figure 4-17 Shear response of the elastic model as function of surface geometry when a normal stress of 14 MPa was applied. The first important result of these simulations is that the elastic model leads to an increase of the shear peak stresses compared to those obtained with the rigid joint model. This can be seen more clearly in Table 4-1 that resumes the average shear peak stresses for all the surfaces S1-S8 for the rigid and the elastic models for a normal stress equal to 14 MPa. Note also that the initial shear stiffnesses of the elastic joints are greater than those obtained for the rigid ones. These results are in agreement with those obtained by [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF] and [START_REF] Ogilvy | Numerical simulation of friction between contacting rough surfaces[END_REF], who suggest that at higher values of the Young Modulus, the contact area between the joints can diminish and as consequence, the shear stiffness of the joints can also diminish (Figure 4 A normal stress of 14 MPa was applied to both joints and then sheared the same distance. Note that the curve of the elastic joint model is shifted upwards with respect to the curve of the rigid model. This means that the work that must be done for the elastic joint model is greater than the one required for the rigid joint model.

On the other hand, the same observation can be made for the residual shear strength: the shear strength obtained with the elastic model is greater than those obtained with the rigid model. One possible explanation of this behavior is that when using the rigid joint model, an external work W R must be done on the system only to overcome the friction between joints due to the roughness of the samples, while in the case of the elastic joint model, an external work W E required to shear the joints must be equal to the energy linked to the friction phenomenon plus the one associated to the elastic deformation of the joint. This implies that the work done for the shear force, in the case of the elastic joint model, must be larger than the one done in the case of the rigid model, for the same level of compression stress; W E >W R . The area under the curve of shear stress vs shear displacement is a measure of the work done by the shear force, as a consequence, when the same compression stress is applied to both models and when both models are sheared at the same distance and in the same direction, the relation W E >W R implies that along the curve of shear stress vs shear displacement of the elastic model is shifted upwards in comparison to the curve of the rigid model. This shift can be seen in Figure 4-19b.

The undulations observed after the shear peak for the rigid model are also present for the elastic model. Nevertheless, the fluctuations of the shear stress obtained with the elastic model are less pronounced as it is possible to see in Figure 4-20 which shows a comparison between both models. Figure 4-21 shows a comparison of the dilatancy of the joints obtained for the rigid and the elastic model. The dilatancy resulting of the elastic model is lower than the one obtained for the rigid model. This behavior can be explained by considering that, in the case of the rigid model, the upper joint must follows the surface roughness of the interface between the joints during shearing, while in the case of the elastic joints, the upper and lower samples tends to deform at the interface and inside the volume as a function of the shear and normal loads applied. The influence of the surface morphology on the shear behavior obtained for the elastic joint model is analogous to the one discussed in the previous section for the rigid model:

 The shear stress and the dilatancy are very similar for both directions of shearing in all the surfaces S1-S8. The similar behavior in the x and y directions can be explained because the roughness of all the surfaces is isotropic [START_REF] Yang | An index for describing the anisotropy of joint surfaces[END_REF], Riss (1997), [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF], [START_REF] Park | A constitutive model for shear behavior of rock joints based on three-dimensional quantification of joint roughness[END_REF].  As the variance of the surface increases, the shear peak and the dilatancy increase; see surfaces DEM S1 compared to surface DEM S2 (surface DEM S2 has lower variance), DEM S3 vs DEM S4, DEM S5 vs DEM S6 (surface DEM S6 has lower variance than surface DEM S5) and DEM S7 vs DEM S8.  As the correlation length increases, the shear peak and the dilatancy decreases; see surfaces DEM S1 vs DEM S3 (surface DEM S3 has higher correlation length), and surface DEM S5 vs DEM S7 (surface DEM S7 has higher correlation length than surface DEM S5).  As the roughness exponent decreases, the shear peak and the dilatancy increase, see surface DEM S1 compared to surface DEM S5 (surface S5 has lower roughness exponent than surface DEM S1), [START_REF] Xie | Fractal effects of surface roughness on the mechanical behavior of rock joints[END_REF], [START_REF] Kwafniewiski | Surface roughness evolution and mechanical behavior of rock joints under shear[END_REF]. Note that the influence of the roughness exponent on the shear peak value and the dilatancy is more pronounced than those of the correlation length and the variance.  For a high value of the roughness exponent which means a low roughness of the surface (case of DEM S1, DEM S2, DEM S3 and DEM S4), the shear peak tends to disappear while for a low roughness exponent, i.e., a high roughness of the surface (case of DEM S5, DEM S6) the shear peak is more pronounced. Moreover, if a high value of the roughness exponent is combined with a high value of the self-affine correlation length and a low value of the height variance (as DEM S4: H H , Lc H , σ 2 L ) then the roughness of the surface is very low and in this case the post peak response is very flat (magnification of the effect of a high value of the roughness exponent). In the same way, if a low value of the roughness exponent is combined with a low value of the self-affine correlation length and a high value of the height variance (as DEM S5: H L , Lc L , σ 2 H ), the roughness of the surface is very strong and leads to a very pronounced peak value of the shear stress (magnification of the effect of a low value of the roughness exponent) [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF].  For a low surface roughness (DEM S4), associated to a post-peak response very flat, there is less dilatancy during shear while for a higher (DEM S5) surface roughness, associated to a very pronounced peak of shear stress-shear displacement response, the shear dilatancy is higher.

As for the rigid model, the residual shear stress obtained for the elastic joint model tends to the same value independently of the original roughness of the joints in contact as shown in Figure 4-22.

Figure 4-22 Effect of the surface roughness and the compression stress on the shear behavior of the elastic joints.

Note that the compression stress applied on the joint has also the same effect that for the rigid joints; at higher compression stress, higher shear peak, higher residual strength, higher shear stiffness and lower dilatancy, as shown in Figure 4 Finally we can see on Figure 4-24 and Figure 4-25 that the distribution of contact points in the elastic joint cases has similar appearance that the one obtained for the rigid joints. For the elastic model there is also formation of channels perpendicular to the shear displacement (in the present case for shear displacements around 1mm) and for very high shear displacements (bigger or equal to 3 mm) the contact points between the upper and lower joints tends to localize only in some areas that coincide whit the highest zones (highest 2D colatitudes, [START_REF] Gentier | Influence of fracture geometry on shear behavior[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF] of the surfaces (red zones in the DEM surface). The main difference observed for the contact points is that there are more contact points or contact area in the case of the elastic joints compared to the rigid one (Figure 4-26). Indeed, as expected intuitively the joints being able to deform (elastic model), there are more probability that they will be in contact in different locations. This behavior is in agreement with the suggestion that for large values of the Young Modulus there is a decrease of the contact area between the joints [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF].

Elastic Rigid Figure 4-26 Comparison of the distribution of contact points for the rigid (rigth) and elastic (left) joint model. Note that there are more contact points when using the elastic model than the rigid one, but the general appearence between both models are very similar (formation of channels and localization of contact points at large shear displacements). The contact points localization observed for large shear displacements coincides with the highest zones of the digital elevation maps of the surfaces S1-S8.

Figure 4-27 Results of the shear simulations for the joints DEM S1-S8 when a compression stress of 14 MPa was applied. The joint's model used is an elastic one and is capable to break.  The undulations of the shear stress after the peak previously observed in the rigid and elastic models, tend to disappear in the elastic-fracture model, reinforcing the idea that the undulations appear when the upper sample follows the morphology of the joint during shearing, that is the case for the rigid model and at lower degree in the elastic model and the elastic-fracture model.  The shear stress and the dilatancy are very similar for both directions of shearing in all the surfaces S1-S8. The similar behavior in the x and y directions can be explained because the roughness of all the surfaces is isotropic [START_REF] Yang | An index for describing the anisotropy of joint surfaces[END_REF], Riss (1997), [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF], [START_REF] Park | A constitutive model for shear behavior of rock joints based on three-dimensional quantification of joint roughness[END_REF].  As the variance of the surface increases, the shear peak and the dilatancy increase; see surfaces DEM S1 compared to surface DEM S2 (surface DEM S2 has lower variance), DEM S3 vs. DEM S4, DEM S5 vs. DEM S6 (surface DEM S6 has lower variance than surface DEM S5) and DEM S7 vs. DEM S8.  As the correlation length increases, the shear peak and the dilatancy decreases; see surfaces DEM S1 vs. DEM S3 (surface DEM S3 has higher correlation length), and surface DEM S5 vs. DEM S7 (surface DEM S7 has higher correlation length than surface DEM S5).  As the roughness exponent decreases, the shear peak and the dilatancy increase , see surface DEM S1 compared to surface DEM S5 (surface S5 has lower roughness exponent than surface DEM S1), [START_REF] Xie | Fractal effects of surface roughness on the mechanical behavior of rock joints[END_REF], [START_REF] Kwafniewiski | Surface roughness evolution and mechanical behavior of rock joints under shear[END_REF]. Note that the influence of the roughness exponent on the shear peak value and the dilatancy is more pronounced than those of the correlation length and the variance.  For a high value of the roughness exponent which means a low roughness of the surface (case of DEM S1, DEM S2, DEM S3 and DEM S4), the shear peak tends to disappear while for a low roughness exponent, i.e., a high roughness of the surface (case of DEM S5, DEM S6) the shear peak is more pronounced. Moreover, if a high value of the roughness exponent is combined with a high value of the self-affine correlation length and a low value of the height variance (as DEM S4: H H , Lc H , σ 2 L ) then the roughness of the surface is very low and in this case the post peak response is very flat (magnification of the effect of a high value of the roughness exponent). In the same way, if a low value of the roughness exponent is combined with a low value of the self-affine correlation length and a high value of the height variance (as DEM S5: H L , Lc L , σ 2 H ), the roughness of the surface is very strong and leads to a very pronounced peak value of the shear stress (magnification of the effect of a low value of the roughness exponent) [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF].  For a low surface roughness (DEM S4), associated to a post-peak response very flat, there is less dilatancy during shear while for a higher (DEM S5) surface roughness, associated to a very pronounced peak of shear stress-shear displacement response, the shear dilatancy is higher.

Figure 4-29 shows the influence of the compression stress on the shear behavior of the elasticfracture joints (shear along the x-axis). The main effect of the compression stress level on the shear behavior is an increase of the magnitude of the shear stress vs. shear displacement response with perhaps a more marked peak. In other words, the shear stress curve shifts upwards when the compression stress increases and conversely. The dilatancy exhibits the contrary effect, i.e., an increase of the compression stress decreases the dilatancy of the joints. These both behaviors have been reported experimentally by several authors [START_REF] Schneider | The friction and deformation behavior of rock joints[END_REF], [START_REF] Bandis | Experimental studies of scale effects on the shear behavior of rock joints[END_REF], [START_REF] Xie | Fractal effects of surface roughness on the mechanical behavior of rock joints[END_REF][START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF] as observed previously.

The compression stress level has also an influence on the fracture behavior of the joints during shear. Indeed, as shown in Figure 4-30, the number of fractures (fractures of bonds) is greater under a stress level of 21 MPa than under 14 MPa, and this is valid whatever the shear displacement.

Figure 4-30 Influence of compression level in fracturing, the surface DEM S7 was sheared in the -x direction and two levels of compression stress were applied: 14 MPa and 21 MPa. The increase of the number of fractures or damaged zones with the compression stress level has also been observed by Gentier (1997), [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF], after Flamand, 2000). Finally, as mentioned previously for the rigid and elastic joint models, the compression stress also influence the initial shear stiffness of the joints, indeed, as the compression stress increases, the initial shear stiffness also increases . Before comparing the results obtained from the DEM shear simulations for the three mechanical models (rigid, elastic and elastic-fracture), let us compare the results obtained from the elasticfracture model with experimental data. Table 4-4 shows data obtained from direct shear experiments extracted from [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF]. Flamand used cylindrical mortar replicas of a natural rock joint which were submitted to normal stresses of 14 and 21 MPa and then sheared (direct shear test). The roughness parameters of the same mortar replicas estimated for [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF] are shown in Table 4-3. [START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF]. *DEM = values obtained from the DEM simulation of the friction of rock joint DEM S1.

Note from Table 4-4 that, despite of the differences of the roughness parameters between the experimental [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF] and DEM joints (Table 4-3), the peak shear stress, residual shear stress and shear peak displacement obtained experimentally and by DEM simulations are very similar. However, the shear stiffness estimated from the DEM simulations is higher than the one obtained experimentally. A possible explanation for this difference could be the mismatching between the upper and lower parts of the joints in the experiments. .

Another interesting aspect related with the shear behavior of rock joints, is the evolution of the contact area as a function of the shear displacement. It is known that during shear, only a fraction of the joint area enters into contact, and is the concentration of such contact areas that affects several aspects of the joints behavior, i.e., channels formation, fracture behavior of asperities or damages zones of the surfaces, electrical and thermal conductivity, propagation of elastic waves, etc. As mentioned before, one approach to study such contact areas in DEM is to study the evolution of contact points between both surfaces of a joint as function of the shear displacement (Figure 4-32).

Figure 4-32 Evolution of contact points during shear for all the DEM surfaces S1-S8, the surfaces were sheared in the -x direction and a constant normal stress of 14MPa was applied.

Figure 4-32 shows the evolution of the contact points during the shear test and this for all the DEM S1 to S8 surfaces (at 14MPa of compression stress). At the onset of the simulation there exist certain quantity of contact points resulting of the compression of the joints again each other. A little shear displacement after, there is a strong increase of the number of contact points [START_REF] Ladanyi | Shear strength and deformability of filled indented joints[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF]. Then, as the shear proceeds, the number of contact points between the surfaces tends to diminish following a power law function of the form = where is the number of contact points as a function of the shear displacement , is the number of contact points at the beginning of the shear simulation and is the exponent of the power law. When the shear displacement is around 3 mm, the number of contact points tends to a constant value for all the joints DEM S1-S8. This shear displacement coincides with the shear displacement for which the residual stress is observable. Note that, as previously observed for the preceding parameters, the evolution of the contact points for all the DEM surfaces are ranged between two extreme behaviors corresponding to the roughest (DEM S5) and the smoothest (DEM S4) surfaces. For all shear displacements with the exception of those corresponding to the residual shear stress regime, the number of contact points of the roughest surface is lesser than the one of the smoothest surface.

The decrease of the contact points observed at the end of the simulation for the rigid and the elastic model is in agreement with the experimental result that at the end of the shear test the contact between surfaces tends to concentrate in some specific locations. As shown previously, the contact The concentration of damages zones during shear has been studied by several authors Riss (1997), [START_REF] Xie | Fractal effects of surface roughness on the mechanical behavior of rock joints[END_REF][START_REF] Gentier | Influence of fracture geometry on shear behavior[END_REF], [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF]. Figure 4-34 shows an example of experimental damaged zones obtained by image analysis treatment by [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF], after Flamand, 2000) when a mortar joint was sheared at constant normal load.

Figure 4-34 Damaged zones of a mortar joint at the end of a direct shear test obtained from [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF], after Flamand, 2000). The black arrow indicates the shear direction. Figure 4-35 shows the evolution of broken bonds or fractures as a function of the shear displacement of the DEM S1 to S8. These evolutions are ranged between two extreme behaviors of fracturing given by the surfaces DEM S4 and DEM S5, corresponding to the flatness surface and roughest surface respectively. This implies that as the joint roughness increases (from DEM S5 to DEM S4) there are more fractures [START_REF] Belem | Fractal analysis of shear joint roughness[END_REF], [START_REF] Xie | Fractal effects of surface roughness on the mechanical behavior of rock joints[END_REF]. Moreover, for each joint, the evolution of broken bonds appears transitional between two asymptotic behavior plotted as dashed lines in Figure 4-35; the onset of the second asymptotic regime seems to correspond with the onset of the plateau of the residual shear stress (Figure 4-29). On the other hand, during the shear tests, the number of fractures is recorded every 100 cycles. The time step dt and the shear velocity v being constant during the shear simulations (as explained in Chapter 3), it can be considered that the number of fractures are recorded as a function of the discrete shear displacement = . In rock mechanics it is common to consider the microfractures as sources of micro seismic events (or acoustic events) because when fractures occurs among the elastic energy stored in the volume of the material surrounding the fracture crack, a part is used to create the new fracture surfaces and the excess of elastic energy is converted in the form of elastic waves (acoustic waves) that propagates through the medium. These elastic waves can be detected experimentally by piezoelectric transducers [START_REF] Moradian | Evaluating damage during shear tests of rock joints using acoustic emissions[END_REF], as shown in Figure 4-36b were the evolution of acoustic emission recorded during a direct shear test is plotted. Figure 4-36a shows the evolution of fractures as a function of the shear displacement for DEM surface S7. [START_REF] Moradian | Evaluating damage during shear tests of rock joints using acoustic emissions[END_REF] Note that the experimental evolution of acoustic emission plotted in Figure 4-36b, [START_REF] Moradian | Evaluating damage during shear tests of rock joints using acoustic emissions[END_REF] is not only related to the micro fractures but also to the friction phenomenon of the asperities. In this sense, if the curves obtained from the numerical simulation and the experimental ones are very similar, they do not reflect necessarily the same underlying physical mechanisms. Another interesting aspect of the fracture phenomenon during shear concerns the distribution of the fracture events (number of micro-fractures). Figure 4-37a shows the frequency of the fracture events obtained for all the DEM surfaces S1-S8 when they are sheared at a compression level of 14MPa. The first thing to note from Figure 4-39 is that most of the external energy supplied to the joints is used to overcome the frictional forces between the discrete elements. As we will see later, this frictional force is strongly related to the joint's roughness. The second major energetic factor is the work done by the compression force; note that this quantity is negative. Since the joints dilate during the shear, i.e., the angle between the compression force and the dilatancy is 180°. The work done by the compression force depends also strongly on the roughness of the joints as detailed in the following. The third most important energetic factor is the elastic energy stored in the bonds and contacts . Note that at the onset of the shear simulation the elastic energy tend to growth and remain approximately constant for large shear displacements; this behavior can be seen more clearly in Figure 4-43b. Finally, the kinetic energy exhibits a very low value compared with the other forms of energy, this can be understand since all the simulations are performed under quasistatic condition and the kinetic energy is mainly dissipated by local damping.

The energy balance of the joint can be written as:

- = ∆ + ∆ 4-4 where
is the energy that enters the system (shear work on Figure 4-39), is the energy that leaves the system (dilatancy work), ∆ is the change in gravitational potential energy of the joint (∆ ≈ 0), and ∆ is the change of the internal energy of the system. As mentioned before, ∆ can be expressed as formed by 4 components: In the present thesis, the internal work estimated at the end of the simulation is labeled as IW end . Figure 4-40 shows the evolution of the different components of the internal energy ∆ as function of the shear displacement for the surface DEM S1 normalized by the internal work IW end . Note from Figure 4-40 that around 80% of the internal energy is expended as friction between particles, while the damping energy, elastic energy and kinetic energy of the particles cover the last 20% of the total internal energy.

∆ = ∆ + ∆ + ∆ + ∆ 4-5 -4.0E-
Among the amount of energy dissipated by the local damping, the energy related to the fracture of bonds represents the major part while the one linked to the kinetic energy dissipation of the other particles (i.e., the particles that are not concerned by the fractures) corresponds to the complement (this later energy is low because the shear test is performed under quasi-static set up). Note that, when a bond breaks, the elastic energy stored in it should be converted theorically in form of kinetic energy of the two particles linked by the bond but, the local damping used in DEM dissipates a large part of this energy. Thus, the energy dissipated by local damping during the simulations can be seen as a picture of the energy dissipated in form of fracture energy, heat and acoustic waves (kinetic energy of the particles) related to the fracture and friction phenomena during an experimental shear test. As mentioned before, the friction energy is strongly dependent on the joint roughness. Indeed, as expected intuitively and shown in Figure 4-42, the energy linked to the friction phenomenon increases as a function of the surface roughness. Moreover as previously observed, the boundary evolutions of the friction energy are associated to smoothest surface (DEM S4, lower limit) and the roughest one (DEM S5, upper limit). Figure 4-43a shows the evolution of the work linked to the normal force applied to the system (normal stress of 14MPa) for all the DEM surfaces S1-S8 (surfaces sheared in the -x direction). As for the friction, the compression work linked to the roughest surface (DEM S5) and the smoothest one (DEM S4), represents the limits of the different evolutions. Indeed, from Figure 4-43a we can see that as the roughness of the joints increases, the work linked to the compression load increases. Moreover, Figure 4-43a reflects also the degree of dilatancy of the joints because the work of the compression load is directly linked to the dilatancy displacement. Finally, the evolution of the elastic energy stored in bonds and contacts is plotted in Figure 4-43b as a function of the shear displacement for all the joints (DEM S1 to S8). For small shear displacements (< 3mm) the elastic energy stored in the bonds and contacts appears dependent on the morphology of the joints. The evolutions related to the smoothest surface (DEM S4) and the roughest one (DEM S5) appear again as the limit evolutions of the energy stored. Nevertheless, for larger shear displacements, this dependency is not very clear. In fact it seems that the elastic energy tends to the same residual value for large shear displacements (> 3mm) corresponding approximately to those for which the shear stress reaches the residual shear stress.

Note in Figure 4-43b that the evolution of the elastic energy during shearing is not a smooth curve but has fluctuation. A source of these fluctuations could be the fracturing of the bonds and the rearranging in the positions of the discrete elements during shearing. Thus, the fluctuations in elastic energy δE elastic must be linked to the fluctuations of the fractures δF (fracture events). Moreover, because the elastic energy stored in the bonds is converted into kinetic energy of the particles when the bond breaks, the fluctuations in kinetic energy δE kinetic must be linked to the fluctuations of elastic energy δE elastic . Figure 4-44a and b show the distribution of δE kinetic and δE elastic for all the DEM surfaces S1-S8, and Figure 4-44c, d ande the corresponding correlation functions for the DEM joint S4 between δE kinetic vs. δF, δE elastic vs. δF and δE kinetic vs. δE elastic respectively. All the fluctuations mentioned above are estimated for two consecutive time steps (or consecutive shear displacement intervals).

Remark from Figure 4-44 that, the distribution of elastic energy and kinetic energy fluctuations, both follow a power law behavior with the joints DEM S5 and DEM S4 being the upper and lower bounds of the distributions behavior. Also note (Figure 4-44c) that δE kinetic and δF are correlated for few consecutive shear displacement steps (one lag = one step of shear displacement between both signals), in fact, for lag = 0 , the correlation coefficient between δE kinetic and δF is 0.136, indicating that a positive change in the number of fractures leads a positive change in the kinetic energy of the set of particles.

Correlation coefficient

Lag distance between signals (µm) Correlation between kinetic fluctuations and elastic fluctuations

On the other hand, the correlation function between δE elastic and δF (Figure 4-44d) takes negative values for small lags, for lag = 0, the correlation function takes the value of -0.181, indicating that a positive change in the number of fractures leads a negative change on the elastic energy of the joints. Finally note that the correlation function between δE kinetic and δE elastic takes the value of -0.139 for a lag = 0 meaning that a decrease in the elastic energy leads an increase of the kinetic energy of the set of particles. Note that the correlation coefficients cf mentioned above are small (cf = ±1 for perfectly correlated signals and cf=0 for uncorrelated signals). However, it is expected that such correlation coefficient increase if the different kinds of energy and number of fractures are register at a lesser number of cycles during the shear simulations (i.e., every 10 cycles or less).

On the other hand, as suggested before, the damping energy is linked with the fracture of bonds during shearing and thus, the fluctuations of fractures δF (fracture events) must be related to the fluctuations of the damping energy δE damping

Correlation coefficient

Lag distance between signals (µm) Correlation between damping fluctuations and fracture fluctuations

Figure 4-45a and b show that the joints with more fractures (i.e., DEM S5) also exhibits more damped energy, also note that joints DEM S4 and DEM S5 correspond to the lower and upper limits of the damped energy evolution of the set of joints. Moreover, similarly to the kinetic and elastic energy fluctuations, the damping fluctuations also present a power law behavior (Figure 4-45c). The correlation function between δE damping and δF (Figure 4-45d) take positive values for small lags, in particular, for lag=0, the correlation function takes the value of 0.368, indicating that a positive fluctuation of the fractures leads a positive fluctuation of the damped energy. Finally, Figure 4-46 shows the correlation function between the kinetic energy fluctuations δE kinetic and the damped energy fluctuations δE damping . As expected, there is a positive correlation between δE kinetic and δE damping (for lag = 0, the correlation function is 0.309) indicating that as the fluctuations of the kinetic energy in the system increases, also increases the proportion of energy dissipated by the local damping. As mentioned before, the correlation coefficient between δE damping and δF and δE kinetic and δE damping could be greater if the sampling steps (number of cycles) for the different kinds of energy and fractures during the simulations diminishes.

Comparison between models (rigid, elastic and elastic-fracture)

Let us now compare the behavior observed for the different models and see the influence of the roughness, elasticity and fracture in a more general way. Figure 4-47 and Figure 4-48 show the shear stress and dilatancy curves of the three models (rigid, elastic and elastic-fracture) sheared along the x-axis under a compression stress of 14 MPa and 21 MPa. From Figure 4-47 and Figure 4-48 it is possible to see that the elastic and rigid models represents the upper and lower limits respectively of the shear stress vs. shear displacement responses. The response of the elastic-fracture models lies between these two limits. The elastic and the elastic-fracture responses are analogous until both reach the shear peak. After the shear peak, both curves start to differentiate; this mean that the deformation of the joints can be considered as elastic until the shear stress peak is reached which is in agreement with the ideas expressed by [START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF]. Note also that, as expected, the initial stiffness of the elastic and the elastic-fracture model are the same, while the initial stiffness of the rigid model is lesser than the one of the latter two models. Moreover, the magnitude of the undulations of the shear stress decreases with the decrease of the Young modulus of the sample (comparison of the rigid and elastic models) and tends to be blurred when the model is able to activate fractures.

-

Correlation coefficient

Lag distance between signals (µm) Correlation between kinetic fluctuations and damping fluctuations

Figure 4-47 Comparisson of the shear behavior observed for the rigid, elastic and elastic-fracture models. The joints were sheared in the -x direction under a compression stress of 14 MPa . In the other hand, the dilatancy of the rigid model and the elastic-fracture model represents the upper and lower limits respectively. The dilatancy of the elastic model lies between these two extreme behaviors. As observed for the shear stress, the difference in dilatancy between the elastic and elastic-fracture model is negligible for short shear displacements <1mm and this difference tends to increase for larger shear displacements. As expected intuitively, the elastic-fracture model leads to a dilatancy lesser than the one of the elastic model, the fractures inducing a decrease of the dilatancy of the joint. In addition, the difference in dilatancy between models is influenced by the roughness of the joints. As the roughness increases, the difference in dilatancy between the models also increases.

From an energetic point of view it is possible to say that, the work done by the shear force decreases as the Young modulus of the joint increases (the area under the curve of shear stress vs. shear displacement is a measure of the work done by the shear force). This means that, for the same roughness and compression stress, rigid rocks require less work to be sheared than softer rocks. Moreover, the work done by the shear force decreases if the model is able to fracture. Thus, a perfectly elastic joint is able to store more elastic energy than an elastic model that can develop fractures, the fracture behavior leading to a limit the quantity of elastic energy that the model/material can store when external loads are applied to the system. So for rocks with the same Young modulus but different compression or tension strength, it is expected that the material corresponding to the lower strength can be sheared easily compared to the material with highest strength.

Finally, the work required by the dilatancy of the joint increases as the Young modulus of the material increases but the fracture tends to decrease the dilatancy work. Thus, a rock mass composed of a soft rock with high strength will require more energy to be sheared than a rock mass composed of a stiff rock with low strength.

Conclusions and future work

In this thesis, the mechanical behavior of rock joints is studied numerically. A methodology to create numerical rock joints having elastic and fracture behaviors equivalent to experimental data (Flamad, 2000) was established using a Discrete Element approach. The numerical joints also have rough selfaffine surfaces, which where created using "numerical molds" produced by the spectral method with different roughness parameters: Hurst exponent H, correlation length L c and height variance σ 2 . Once the elastic and fracture properties of the DEM joints were calibrated, the set of joints were submitted to a direct shear tests at constant normal load. Two levels of compression stress were applied (14 MPa and 21 MPa) and three different mechanical behaviors were implemented (rigid, elastic and elastic fracture). The influences of the roughness, compression stress, shear direction and mechanical behavior on the dilatancy and shear behaviors are analyzed and lead to several conclusions:

 From the numerical method based on the Fourier Transform proposed in Chapter 2, it is possible to generate isotropic and anisotropic self-affine surfaces characterized by the roughness exponent, the correlation length and the height variance of the surfaces. The anisotropy of the surfaces can by produced by changing the values of the roughness exponent, the correlation length and the height variance. These numerical self-affine surfaces can be used as molds to generate self-affine surfaces from DEM.

 DEM can be used to model the elastic and fracture behavior of quasi-brittle materials, such as rocks, mortar, etc. Particularly, it is possible to define a Representative Elementary Volume (REV) formed by approximately 1000 particles to calibrate the elastic and fracture properties of the numerical joints which have a symmetric elastic behavior in compression and tension (same Young modulus and Poisson ratio in compression and in tension). Moreover, the REV also behaves as an isotropic material, meaning that, its Young modulus E, Poisson ratio v and its shear modulus G are related by: = ⁄ (2(1 + ).

 The numerical simulations are performed on the basis of strains equivalent to the experimental ones but using a Young modulus much smaller that the real rock. Since the Young modulus of the REV and the critical time step of integration during the simulations are inversely related, the DEM simulations can be made in relatively short time (3 to 4 days).

 Because the discrete character of the DEM, the models produced have an intrinsic porosity that it is linked to the size distribution of the particles used. In general, as the size of the particles diminishes, it is possible to construct models with lower porosity. However, there are limits to the porosity that can be achieved and the decreasing of the porosity (i.e., the increase of the number of particles) lead to extensive computation time. At some values of porosity, the particles can be arranged without touching; at other values, there will be overlap between particles. In extreme cases of DEM models with very low porosity, the inter-particle forces can be so strong that could influence the fracture behavior of the model in the sense that, when the particles are bonded together and when a bond breaks, the elastic energy stored in the bond and then transformed in kinetic energy of the particles, can be high enough to produce a large avalanche of fractured bonds, an so, fracturing all the model in unrealistic manner.

 To calibrate the fracture properties of the REV, an elliptical fracture model based on strengths is implemented. The elliptical fracture model has only one independent variable (tension strength of the bond) which was calibrated taking into account the internal stresses in the REV resulting from the its genesis process. As a result of the implementation of the elliptical fracture model, the compression and tension responses of the REV exhibit behaviors in agreement with those expected for a quasi-brittle material. Moreover, the ratio between compression strength and the tensile strength obtained from the REV is around 10 as is expected in natural rocks. In tension, the broken bonds are located within a more or less narrow region band oriented perpendicular to the tensile loading axis which is in agreement with the preferential orientation of the micro-cracking in mortar. In compression, the microcracks spread in the whole volume without apparent preferential orientation.

 It is possible to generate self-affine surfaces in DEM using as molds the numerical self-affine surfaces produced by the Fourier Transform method. However, because the intrinsic discrete character of the model used, there is an intrinsic porosity that modifies the roughness of the DEM joints compared to the original surfaces used as molds. In particular, the DEM self affine surfaces produced in such a way have approximately the same height variance and correlation length that the surface used as mold, but, the roughness exponent of the DEM surfaces is lower than the one of the surface used as mold. Indeed, in order to study the influence of the roughness exponent on the mechanical response of the joint, the particles located out of the expected rough surface (and mainly linked to the porosity of the discrete model) are identified so as to estimate their contribution in the whole contact and friction works. However, in this study, the contribution of these "porosity stresses" is negligible.

 In the compression or closure test, as expected intuitively, the elastic-fracture models deform (or close) more than the elastic models which, in turn, deform (or close), more than the rigid model. A possible explanation for that is that the fracture of bonds in elasticfracture models allows a rearrangement of particles (relative displacements between neighbors particles) which lead to an increase of the joint closure compared to the elastic model. Moreover, for each joint, the joint normal stiffness Kn is higher for the rigid model than for the elastic and elastic-fracture models, for which the joint stiffnesses are almost the same. Finally, the contact points between both surfaces of the joint are uniformly distributed in the whole joint surfaces. A possible explanation for this behavior could be that upper and lower joint surfaces are very similar (i.e., they match very well).

 Concerning the shear simulations, the main conclusion is that the roughness of joints (described by three roughness parameters; roughness exponent, correlation length and height variance) plays an important role in the mechanical response of the joints. Moreover, when the joint roughness is isotropic, its mechanical behavior also is isotropic (i.e., independent of the shear direction). The influence of the roughness exponent, correlation length and height variance is resumed in the following points and these conclusions apply for the rigid, elastic and elastic-fracture models:

 The influence of the self-affine roughness exponent on the shear peak value and the dilatancy is more pronounced than those of the correlation length and of the height variance.

The shape of the shear stress vs. shear displacement curve, and especially, the existence of a shear peak, depends mainly on the value of the self-affine roughness exponent, i.e., for a high value of the roughness exponent (a low roughness of the surface), the shear peak tends to disappear while for a low roughness exponent (a high roughness of the surface), the shear peak is more pronounced. If a high value of the self-affine roughness exponent is combined with a high value of the self-affine correlation length and a low value of the height variance, then, the roughness of the surface is very low (flatness surface) and then, the post peak response in shear is very flat: the combination of a high value of the roughness exponent with a high correlation length and a low height variance magnify the relevant effect o a high roughness exponent. On the contrary, if a low value of the roughness exponent is combined with a low value of the self-affine correlation length and high value of the height variance, then, the roughness of the surface is very strong (roughest surface) and then, leads to a very pronounced peak value of the shear stress: the combination of a low value of the roughness exponent with a low correlation length and a high height variance magnify the effect of a low roughness exponent.

 The height variance has a direct relationship with the shear response of the joints, i.e., as the variance of the surface increases, the shear peak and the dilatancy. However, its effect is less pronounced that the effect of the roughness exponent.

 The self-affine correlation length has an inverse relationship with the shear response of the joints, i.e., as the correlation length increases, the shear peak and the dilatancy decreases. As for the height variances, the effect of the self-affine correlation length is less pronounced that the effect of the roughness exponent.

 The shear stress vs. shear displacement curves obtained for the rigid model exhibits periodic fluctuations in the post peak regime. These undulations tens to disappear as the Young modulus of the joints decreases (elastic model) and are no longer noticeable when the joints are able to fracture (elastic-fracture model). We argued that this phenomenon can be related to the average size of the discrete particles.

 The residual shear strength (measured at 5mm of shear displacement) is very similar for all the DEM S1-S8 surfaces indicating that the residual shear strength is independent of the initial surface geometry. This behavior was observed for the three mechanical models: rigid, elastic and elastic-fracture model.

 For the elastic and elastic-fracture models, the initial shear stiffness of the joints, expressed by the slope of the shear stress vs. shear displacement at the onset of the curves, is approximately the same for all the DEM surfaces, indicating that this stiffness is quasiindependent of the roughness of the joint. However, when the rigid model is used, the initial shear stiffnesses are lesser than the corresponding ones obtained for the elastic and elastic fracture model. For the rigid model, the initial stiffness is linked to the joint roughness.

 In all the shear tests simulations, no appreciable contractancy was observed for all the DEM joints, this independently of the compression stress. We argue that this behavior is due to the fact that the joint's surfaces match at almost all the scales of observation.

 At the onset of the shear simulations, the contact points distribution is almost homogenous in the joint plane. For small shear displacements, the number of the contact points between the joint faces increases, then, as the shear displacement increases, the number of contact points decrease following a power law behavior until a minimum average number of contact points. This decrease in contact points leads to the formation of channels perpendicular to the shear direction. At moderate and long shear displacements (more than 2 mm), the contact points tend to localize in narrow regions that coincide with the highest zones of the joint. However, the contact zones do not necessarily correspond to all the higher areas of the joint. The main difference observed for the contact points is that there are more contact points or contact area for joints having the elastic and elastic-fracture model compared to the rigid model joints. Indeed, as expected intuitively the joints being able to deform (elastic model and elastic-fracture model), the probability that they will be in contact in different locations increase. The evolution of the contact points for all the DEM surfaces are ranged between two extreme behaviors corresponding to the roughest (DEM S5) and the smoothest (DEM S4) surfaces. For all shear displacements with the exception of those corresponding to the residual shear stress regime, the number of the contact points of the roughest surface is lesser than the one of the smoothest surface.

 The effect of the normal load on the shear response can be summarized in the following points:  There is a direct relationship between the compression stress, the shear peak , the residual shear stress, and the number of contact points (or contact area) i.e., as the normal load applied to the joint increases, the number of contact points between both faces of the joint increases and the shear peak and the residual shear stress also increases. Moreover, when the joint is able to fracture, an increase in the level of compression stress leads to an increase in the number of fractures in the joint during shear.  The dilatancy and the level of compression stress exhibit an inverse relationship: as the normal stress applied to the joint increases, the dilatancy of the joints decreases.

 The Young modulus of the joints affects its shear behavior. For a fixed roughness and level of normal load, a perfectly elastic joint (without damage) is more resistant to shearing than a perfectly rigid joint (for which the value of the Young modulus tends to infinity). Indeed, as the Young modulus of the joints increases, the shear stress vs. shear displacement curve tends to offset downwards. Moreover, when the joint is able to fracture, its shear behavior lies between the elastic and rigid behavior. In fact, the development of fractures (breaking of bonds) imposes an upper limit of the amount of elastic energy than can be stored in the joints.

 The dilatancy is also affected by the elastic-fracture behavior of the joints. As the Young Modulus of the joints increases (the rigid joint being an extreme case), its dilatancy also increases. If the joint can fracture, its dilatancy is less compared with a perfectly elastic joint with the same roughness.

 When the elastic-fracture model is implemented, the joints develop fractures in the entire simulations. The evolution of broken bonds or fractures as a function of the shear displacement are ranged between two extreme behaviors of fracturing given by the surfaces DEM S4 and DEM S5, corresponding to the flatness surface and roughest surface respectively. This implies that as the joint roughness increases (from DEM S5 to DEM S4) the number of fracture also increases. Moreover, for each DEM joint, the evolution of broken bonds appears transitional between two asymptotic behaviors; the onset of the second asymptotic regime seems to correspond with the onset of the plateau regime of the shear stress (i.e., residual shear stress).

 If the breaking of bonds is seeing as source of acoustic events, the probability density function of the breaking bonds follows an exponential function similar to the Gutenberg-Richter Law. Moreover, the probability density function of breaking bonds obtained from the present simulations are very similar to the experimental acoustic emission obtained in concrete [START_REF] Saliba | Apport de l'émission acoustique dans la compréhension et la modélisation du couplage fluage-endommagement du béton[END_REF] besides the boundary conditions of the experiments and the simulations are different (direct shear test at constant normal load, and, uni-axial compression respectively). Another interesting issue is that the probability density functions of fracture events obtained from the shear simulations of the surfaces DEM S1-S8 tends to collapse in one single curve. However, for small magnitudes of fracture events, the distribution of the events is independent of the surface roughness, for bigger magnitudes of events, the magnitude of the fracture events increases as a function of the surface roughness. A possible explanation of this behavior is that, an increase of the roughness could involve (1) an increase of the asperities impeding the shear displacement and consequently (2) an increase of the fractures as the shear progresses.

 From an energetic point of view, the energy supplied to the joints by the shear and normal loads is transformed in the following forms of energy during shear: work done by the frictional forces between the discrete elements , elastic energy stored in the contacts and bonds , kinetic energy of the system of particles and the energy dissipated by the model through the local damping .

 From the energy balance of the joints during shear, we observe that most of the external energy supplied to the joints is used to overcome the frictional forces between the discrete elements. This frictional force is strongly related to the joint's roughness. In fact, the energy linked to the friction phenomenon increases as a function of the surface roughness. The boundary evolutions of the friction energy are once again associated to smoothest surface (DEM S4, lower limit) and the roughest one (DEM S5, upper limit).

 The second major energetic factor is the work done by the compression force. This quantity is negative since the joints dilate during the shear, i.e., the angle between the compression force and the dilatancy is 180°. The work done by the compression force depends also strongly on the roughness of the joints and is directly related to the dilatancy response.

 The third most important energetic factor is the elastic energy stored in the bonds and contacts. For small shear displacements (< 3mm), the elastic energy stored in the bonds and contacts appears dependent on the morphology of the joints. The evolutions related to the smoothest surface (DEM S4) and the roughest one (DEM S5) appear as the limit evolutions of the elastic energy stored. Nevertheless, for larger shear displacements, this dependency is not very clear. In fact it seems that the elastic energy tends to the same residual value for large shear displacements (> 3mm) corresponding approximately to those for which the shear stress reaches its residual value.

 Finally, the kinetic energy exhibits a very low value compared with the other forms of energy. This can be understood since all the simulations are performed under quasi-static condition and the kinetic energy is mainly dissipated by local damping of the discrete model behavior.

Among the amount of energy dissipated by the local damping, the energy related to the fracture of bonds represents the major part while the one linked to the kinetic energy dissipation of the other particles (i.e., the particles that are not concerned by the fractures) corresponds to the complement (this later energy is low because the shear test is performed under quasi-static set up).

 The damped energy (which is in its major parts linked to the fracture of bonds) increases monotonically during the whole simulations. This phenomenon emphasizes that the energy linked to fracture does not exhibit stationary regime. Moreover, the damped energy appears dependent on the joint roughness. The joints DEM S5 and DEM S4 correspond to the upper and lower limits respectively of the energy damped evolution of the set of joints DEM S1-DEM S8.

 The fracture events, the fluctuations of the damped energy, the kinetic energy fluctuations and elastic energy fluctuations are interrelated. Indeed, the fracture events lead to a positive increase of the damped energy (fracture energy), an increase of the kinetic energy and a decrease of the elastic energy of the set of particles in the joints. In the other hand, the probability density function of the damped energy fluctuations, kinetic energy fluctuations and elastic energy fluctuations follows a power law behavior for which the surfaces DEM S4 and DEM S5 appear as the lower and upper limits.

Perspectives

It was shown in this thesis that DEM can be used to reproduce the quasi-brittle behavior of rocks or mortar). On this basis, the discrete character of DEM is particularly adapted to study the friction behavior of rock joints. The agreement of the results obtained from the simulations of the closure and friction behaviors of joints with the experimental ones, allow considering the study of another aspects of the fracture and friction phenomena of rock joints in the future:

 Firstly, it is expected that during the friction of rock joints, the rock-material located far away the joint surfaces only deforms elastically (without fracturing) while the material located close to the joint surface is deformed elastically (or plastically) and can be fractured. The extend of the transitional zone between these two mechanical behaviors can be determined using DEM. Once determined the characteristics of such transitional zone (or friction process zone), the elastic deformations of the material located far away the joint surfaces and the rock material located close to the joint surfaces can be studied using a coupled Finite Element Method -Discrete Element Method (FEM-DEM) approach. In the FEM-DEM approach, computational time can be diminished if the material close to the joints surfaces is modeled using discrete particles and the material far away the joint surfaces is modeled using FEM. A similar approach can be used to study the fracture behavior of quasi-brittle material if the material located far away the fracture process zone (FPZ) is modeled using FEM and the material inside the FPZ is modeled using DEM. The heat generation during friction and its propagation throw the rock material could be also studied with a DEM-DEM approach.

 The comprehension of the hydro-mechanical behavior of rock joints (propagation fluids throw the channels formed during shearing) is of capital importance for various civil and engineering projects, i.e., geothermal energy exploitation, extraction of oil and gas from the earth crust, CO 2 storage, nuclear waste disposal, etc. In order to study such phenomena, a FEM-DEM approach could be implemented.

 To generalize the results obtained in this thesis, it necessary to study of the influence of the structure of the joints surface (planar distribution of peaks and valleys) on its shear behavior and on the localization and evolution of contact points during friction. As mentioned in Chapter 2, the joint surface structure is linked to the phase distribution used to generate the numerical self-affine surfaces.

 In the other hand, natural rock joints exhibit normally some anisotropy of its roughness. The study from DEM of the influence of the surface roughness anisotropy of the joints on its shear behavior (anisotropy in roughness exponent, correlation length and height variance) can be also interesting.

 Natural joints can have gouge or deposited material between their surfaces that can affect its shear behavior. Thus, in order to understand the influence of such gauge materials during friction (i.e., on the shear resistance of the joints), gouge material with different characteristics (i.e. material type, granulometry, quantity of gouge, etc.) can be modeled with DEM and carry out friction simulations. Moreover, the influence of the gouge material generated by the detachment of particles from the joint surfaces can be studied by using DEM joints with different compression and tension strengths, which in turn, can be calibrated by adjusting the parameters of the fracture criterion (elliptical criterion) used in this thesis.

 Also it is essential to analyze the shear behavior of DEM joints subjected to several values of normal load. Of particular interest is the estimation of the shear stress peak vs. normal load curves for the DEM joints used in this thesis. In this way, the study of the influence of the surface roughness of joints on shear stress peak vs normal load curves can be also particularly interesting. Moreover, different boundary conditions can be implemented for the friction simulations, such as constant global stiffness or dilatancy imposed instead of normal load boundary conditions. The use of constant global stiffness as boundary conditions will allow the study the friction phenomena of rock joints located dipper from the earth crust.

 Moreover, it is known that because weathering (or shear displacement), the faces of natural rock joints do not match perfectly and the degree of matching affects its shear behavior. DEM joints could be easily shifted to different quantities in order to study the influence of matching of a joint on its shear behavior.

 A detailed study of the scales effects on rock friction can be performed by carried out simulations of shear test from different sizes of DEM joints. In order to diminish the simulation times, these scale-effects study can be performed from coupled FEM-DEM simulations. Furthermore, the influence of the shear speed velocity on the shear behavior of rock joints could be also studied from DEM simulations. The development of better analytical models can achieved using the results obtained from the studies mentioned above. These models could include the combined effects of the elasticity, fracture, surface roughness, joint size, etc.

 Finally, the stability of rock masses is nowadays monitored by studying its acoustic emission.

In order to improve such techniques, the generation and propagation of elastic waves through rock mass, especially, those generated by the joints during shear can also be interesting phenomena studied with DEM. As was exposed in Chapter 4 in this thesis, the acoustic emissions can be linked to the fractures. Nevertheless, dedicated studies are needed to better understand the correlations between fractures and acoustic emission, and especially the connections between the fluctuations of kinetic energy, elastic energy and fractures with the propagation and production of these acoustic emission (or elastic waves).

The accurate understanding of the production and propagation of elastic waves from fracture could be of major interest for the study or rock masses stability but also for earthquake phenomena and mining industry problematic. For instance, in gold (ore) mining industry, it is important to extract the maximum of crushed ore with a minimum of explosives. Thus it could be interesting to study the elastic waves propagation throw rock mass when a punctual load (explosive) is applied and how the propagation of these elastic waves influence the fracture behavior of the rock mass.

Figure 0 - 1

 01 Figure 0-1 Layers of the earth (http://www.mysciencebox.org/).

Figure 0 - 2

 02 Figure 0-2 Tectonic Plates (http://science.taskermilward.org.uk/)

Figure 0 - 3

 03 Figure 0-3 Typical stress-strain curve for an intact rock sample

Figure 0 - 4

 04 Figure 0-4 Massive rock fractured by natural process (http://coloradogeologicalsurvey.org/).

Figure 0 - 5

 05 Figure 0-5 Modes of fracture.

Figure 0 - 6

 06 Figure 0-6 Bedding planes (http://en.wikipedia.org/).

Figure 0 - 7

 07 Figure 0-7 Folds formation near Moruya, New South Wales, Australia (http://en.wikipedia.org/).

Figure 0 - 8

 08 Figure 0-8 The Piqiang Fault, China (http://en.wikipedia.org/).

Figure 0 - 9

 09 Figure 0-9 Shear zone (structuralgeo.wordpress.com).

Figure 0 -

 0 Figure 0-10 Example of a dyke (http://coloradogeologicalsurvey.org/).

Figure 0 -

 0 Figure 0-11 Rock Joints (http://whattherockstellus.blogspot.mx/).

Figure 0 -

 0 Figure 0-12 Quantification of discontinuity occurrence along a sample line.

Figure 1 - 1

 11 Figure 1-1 Stiffness of a rock joint when submitted to a) compression load, and b) shear load,[START_REF] Barton | Strength, deformation and conductivity coupling of rock joints[END_REF].

  Figure 1-2 a) Experimental configuration of the joint's closure test, b) experimental results of a joint's closure,[START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF] after[START_REF] Flamand | Validation d'un modèle de comportement mécanique pour les fractures rocheuses en cisaillement[END_REF].

Figure 1 -

 1 Figure 1-3 a) configuration of the direct shear test, and b) triaxial test to measure the shear strength of rock joints[START_REF] Jaeger | Friction and sliding of joints[END_REF].

Figure 1 - 5

 15 Figure1-5 Size effects on rock joints[START_REF] Bandis | Experimental studies of scale effects on the shear behavior of rock joints[END_REF].

Figure 1 - 6

 16 Figure 1-6 Effect of boundary conditions (dilatancy allowed or not) on shear behavior of rock joints.

Figure 1 - 7

 17 Figure 1-7 Shear peak strength vs. normal stress for a dry smooth rock joint, Jaeger and Rosengren (1969).

Figure 1 -

 1 Figure 1-8 a) Typical shear stress-shear displacement curve obtained for a rough rock joint, and b) evolution of shear peak stress and residual shear strength as function of the normal stress.

Figure 1 -

 1 Figure 1-9 a) Representation of smooth joint, b) smooth joint inclined at an angle i respect to the shear direction, c) rough joint with saw-toot asperities inclined at an angle i and d) rough joint with sinusoidal asperities, after Patton (1966).

Figure 1 -

 1 Figure1-10 Anisotropic shear behavior of rock joints induced by anisotropic surface morphology[START_REF] Brown | Shear strength characteristics of Delabole slates[END_REF].

Figure 1 -

 1 Figure 1-11 Patton model (1966).

Figure 1 -

 1 Figure 1-12 Friction components of rock joints[START_REF] Bandis | Experimental studies of scale effects on the shear behavior of rock joints[END_REF].

Figure 2 - 1

 21 Figure 2-1 Set of joints aligned vertically in a rocky mountain (Chipinque park in Monterrey, México).

Figure 2 - 2

 22 Figure 2-2 Scheme representing two points of a rough profile separated by a distance τ.

  Figure 2-3 Schematic representation of the a) structure function ( ) and b) variogram ( ) of a self-affine profile. 2.2.2.3. Anisotropy

Figure 2 - 4

 24 Figure 2-4 Structure function of an anisotropic surface in Hurst exponent. Results taken from Ponson (2006).

Figure 2

 2 Figure 2-5 a) Weathering of tuff highlighting the joint pattern and producing anisotropic rough surfaces, and b) rough surface of slate showing marked anisotropy caused by the fabric structure of the rock (from http://en.wikipedia.org/).

Figure 2 - 6

 26 Figure 2-6 Schematic representations of Abbott's curve for two different profiles.

Figure 2 - 7

 27 Figure 2-7 Basic relations among stochastic process, covariance function, spectral density function and amplitude spectrum, taken from Pardo-Iguzquiza (1993).

Figure 2 -

 2 Figure 2-10 2D representation of the power spectrum density. The low frequencies (u=0, v=0) are in the center of the matrix. Each point P(u,v) corresponds to a value of the Power spectral density of the surface.

Figure 3

 3 Figure 3-1 (a) Ball-ball contact and (b) ball-wall contact (modified from the User Manual of PFC-3D).

Figure 3 - 2

 32 Figure 3-2 Determination of normal direction for ball-wall contact (modified from the User Manual of PFC-3D).

Figure 3 -

 3 Figure 3-3 force displacement behavior of (a) normal component of contact force and (b)shear component of contact force (modified from the User Manual of PFC-3D).

Figure 3 - 4

 34 Figure 3-4 Schematic representation of parallel bond (modified from the User Manual of PFC-3D).

Figure 3

 3 Figure 3-5 (a) single mass-spring system, (b) infinite serie of springs and masses (modified from the User Manual of PFC-3D).

Figure 3 - 6

 36 Figure 3-6 Size distribution of the spherical particles used in the simulation

Figure 3 - 7

 37 Figure 3-7 Schematic example of REV generation

Figure 3 - 8 Figure 3 -Figure 3 -

 3833 Figure 3-8 Explanation of the protocol of measures, the particles are separated in groups where: set 1= red particles, set 2 = blue, set 3 = pink, set 4 particles located in the opposite face of set 3, set 5 = orange and set 6 particles located in the opposite face of particles in set 5.

Figure 3 -

 3 Figure 3-10 Distribution of a) normal stresses and b) shear stresses developed in the parallel bonds of the REV once the boundary walls are removed and the mechanical equilibrium of the system of particles of the REV is reached.

Figure 3 -

 3 Figure 3-11b exhibits the responses obtained during the compression and tensions simulations according to the fracture criterion defined by Eq. 3-46 and plotted in Figure 3-11a. Note that in agreement with the experimental behavior of a quasi-brittle material such as mortar, the compression and tension responses plotted in Figure 3-11b initiate by a linear part associated to the elastic behavior and are followed by a non linear behavior linked to the progressive damage of bonds until the compression and tension strengths are reached respectively.

Figure 3 -

 3 Figure 3-13 Intact bonds, a) before the simulations, b) at the end of tension test and c) at the end of compression simulation.

  Figure 3-14 Variograms in simple a) and log-log b) plot of the 8 synthetic surfaces S1-S8 used to generate the DEM joints.

Figure 3

 3 Figure 3-15 Elevation maps of synthetic surfaces S1-S8.

Figure 3 -

 3 16 shows the dimension of the joint's DEM model.

Figure 3 -

 3 Figure 3-16 Dimensions of the joint's DEM model.

  of the geometrical boundary of the model interact with the discrete elements by means of their normal and tangential stiffness. 3. The DEM model is cycled until the mechanical equilibrium of the discrete particles is reached. 4. Once at equilibrium, the spherical particles in contact are bonded together by the installation of simple (springs) and parallel (elastic bars) bonds. The mechanical and failure properties used for the simple and parallel bonds are the same of those estimated for the REV. 5. The boundary walls are removed and the model is cycled until the particles reach their equilibrium positions.

Figure 3 -

 3 Figure 3-17 a) Mesh generated from triangular walls. This mesh serves as a mold to generate the rough surfaces in DEM, b) Scheme of the generation of discrete elements in the region bounded by the joint's mesh and the boundary walls.

  Figure 3-18 a) mesh of the synthetic surface used to generate the DEM joint, b) synthetic mesh and DEM joint generated, c) resultant DEM joint, d) regular grid to obtain the height of the DEM surface, e) elevation map of the synthetic surfaces used as a mold, f) elevation map of the DEM surface.

  Figure 3-19 Abbot's modified curve to find the mean plane of the DEM surface S1, a) Relative vertical displacement of synthetic surface and the DEM surface, b) synthetic and DEM surfaces in contact and c) plot of the number of contacts between vs relative vertical displacement of the synthetic surface and the DEM surface

Figure 3 -

 3 Figure 3-21c shows in red color the spherical particles that form the DEM self affine joint surface, and in Figure 3-21b is presented the variogram of the DEM self-affine joint and the synthetic joint with different intensities of the added noise in a log-log plot.

Figure 3 -

 3 Figure 3-22 Elevation maps of DEM surfaces generated with the synthetic surfaces S1-S8 as molds.

Figure 3 -

 3 Figure 3-22 continuation Elevation maps of DEM surfaces generated with the synthetic surfaces S1-S8 as molds.

Figure 3 -

 3 Figure 3-23 Scaling properties of DEM surfaces generated with the synthetic surfaces S1-S8 as molds.

  Figure 4-1 a) schematic representation of the ensemble of the upper and lower samples forming the rock joint, b) numerical joint generated with DEM.

Figure 4 - 2

 42 Figure 4-2 The clumps formed by the balls in blue-color are used to apply the load during the compression simulations.

Figure 4 - 3

 43 Figure 4-3 Results of the DEM simulations of the closure test for the joints DEM S1 to DEM S8. Each joint were simulated using the rigid model, elastic model and elastic-fracture model.

  Figure 4-4 a) Contact points during closure test for joint DEM S2, b) evolution of contact area as function of joint closure for rigid, elastic and elastic-fracture models for the joint DEM S2, c) evolution of broken bonds as function of joint's closure for the joint DEM S2 (elastic-fracture model).

Figure

  Figure 4-6 Plan of simulations for the direct shear test.

Figure 4 - 8

 48 Figure 4-8 Results of the direct shear simulations for the rigid model. Surfaces DEM S1-S8 where sheared in two perpendicular directions x (continuous line) and y (broken line). The normal stress applied was 21 MPa.

  observations that we can make from Figure 4-7 and Figure 4-8:

Figure 4 -

 4 Figure4-9 shows the shear response and dilatancy in the x direction when a compression stress of 14 MPa and 21 MPa are applied. In this figure it is easier to see the influence of the surface roughness in the shear stress and dilatancy responses of the joints.

Figure 4 - 9

 49 Figure 4-9 Influence of geometry in the shear response. The surfaces S1-S8 where sheared in the -x direction at two different compression stresses of 14 MPa and 21 MPa.

Figure 4 -

 4 Figure 4-10 Influence of the normal stress on the initial shear stiffness for the rigid joints DEM S1-S8.

Figure 4 -

 4 Figure 4-11 Experimental results of the shear stress response of a rock joint as function of the geometry when subjected to various compression stress levels Yang (2010).

  Figure 4-12 a) Experimental results of the contractancy-dilatancy of a rock joint subjected to a shear test under various compression stress levels Yang (2010), b) Dilatancy curves obtained for the rigid model when a compression stress of 14 MPa was applied.

Figure 4 -

 4 Figure 4-13 Contact points distribution in the joint as function of the shear displacement for the surface DEM S4 subjected to a compression stress of 21 MPa. The white arrows show the direction of movement of the upper joint surface along the -x axis. The colored figure is the digital elevation map of the surface DEM S4.

Figure 4 -

 4 Figure 4-14 Contact points distribution as function of the shear displacement for the surface DEM S4 (σ= 21 MPa ). The white arrows show the direction of movement of the upper joint surface along the -y axis. The colored figure is the digital elevation map of the surface DEM S4.

Figure 4 -

 4 Figure 4-15 Example of numerical results of the formation of channels as function of shear direction Nemoto (2009).

Figure 4 -

 4 16b shows an example of experimental results obtained by image analysis of sheared surfaces[START_REF] Marache | Comportement mécanique d'une fracture rocheuse sous contraintes normale et tangentielle[END_REF].

Figure 4 -

 4 Figure 4-16 a) Contac point distribution as function of the shear displacement for the surface DEM S4 subjected to 21 MPa of compression stress, b) damaged zones on a mortar sample after 5 mm of shear displacement at which a was applied a compression stress of 21 MPa Marache (2002).

Figure 4 -

 4 Figure 4-18 Shear response of the elastic model as function of surface geometry when a normal stress of 21 MPa was applied.

  Figure 4-19 a) contact area between joints as function of the Young Modulus for a compression of 5 MPa, modified from Marache (2002), b) Comparison the shear response between the rigid model and the elastic model for the joint DEM S1.A normal stress of 14 MPa was applied to both joints and then sheared the same distance. Note that the curve of the elastic joint model is shifted upwards with respect to the curve of the rigid model. This means that the work that must be done for the elastic joint model is greater than the one required for the rigid joint model.

Figure 4 -

 4 Figure 4-20 Comparison of the shear behavior between the rigid model and the elastic model.

Figure 4 -

 4 Figure 4-21 Comparison of the dilatancy of the rigid and the elastic joints models.

  -22 and Figure 4-23.

Figure 4 -

 4 Figure 4-23 Influence of the normal stress on the initial shear stiffness for the elastic joints DEM S1-S8.

Figure 4 -

 4 Figure 4-24 Distribution of contact points obtained with the elastic joint model as function of the shear displacement for the surface DEM S4. The upper joint move in the -x direction and a compression stress of 21 MPa was applied. The white arrows correspond to the direction of the shear along the x-axis.

Figure 4 -

 4 Figure 4-25 Distribution of contact points obtained with the elastic joint model as function of the shear displacement for the surface DEM S4. The upper joint move in the -y direction and a compression stress of 21 MPa was applied. The white arrows correspond to the direction of the shearing displacement along the y-axis.

Figure 4 -

 4 Figure4-28 Results of the shear simulations for the joints DEM S1-S8 when a compression stress of 21 MPa was applied. The joint's model used is an elastic one and is capable to break.

  very similar to both previous cases, can be remarked from Figure4-27 and Figure4-28:

Figure 4 -

 4 Figure 4-29 Influence of compression stress level in the shear behavior of joints.

Figure 4 -

 4 Figure 4-31 Influence of the normal stress on the initial shear stiffness for the elastic-fracure joints DEM S1-S8.

Shear

  S2 DEM S3 DEM S4 DEM S5 DEM S6 DEM S7 DEM S8

  in the highest zones of the surfaces as shown in Figure4-33a and b. As expected intuitively, we can observe a correlation between the location of contact points (Figure4-33a and b) and the location of the detached particles, resulting for the fracture of bonds (Figure4-33c and d).

Figure 4 -

 4 Figure 4-33 Contact points when the joint was sheared in the a) -x direction and b) sheared in the -y direction. Detached particles during shear for the surface DEM S1 subjected to a compression stress of 14 MPa sheared in c) -x direction, d) sheared in -y direction. Dark points denote the location of detached particles.

Figure 4 -

 4 Figure 4-35 Evolution of broken bonds as the shear simulation progress for all the joints DEM S1-S8.

  36 Evolution of fractures as function of the shear displacement, a) results obtained from the shear simulation of the surface DEM S7 when a compression stress of 14 MPa was applied, b) evolution of acoustic emission of an experimental direct shear test obtained by

Figure 4

 4 Figure 4-37 a) Distribution of the fracture events/(number of fractures) during shear tests for the surfaces DEM S1-S8, b) Experimental acoustic emission distribution during the fracture of a concrete sample[START_REF] Saliba | Apport de l'émission acoustique dans la compréhension et la modélisation du couplage fluage-endommagement du béton[END_REF].

Figure 4 -

 4 Figure 4-39 continuation Evolution of the different kinds of energy involved in the shear test.

Figure 4 -

 4 Figure 4-40 Distribution of the internal energy as function of the shear displacement.

Figure 4 -

 4 Figure 4-41 shows the evolution of the damped energy as function of the shear displacement for all the DEM joints S1 to S8. It is possible to see from Figure 4-41 that the damped energy of the joints DEM S5 and DEM S4 corresponds to the upper and lower limits respectively of the energy damped evolution of the set of joints DEM S1-DEM S8.

Figure 4 -

 4 Figure 4-42 Evolution of the friction work as function of the joint's roughness.

Friction

  43 Evolution of a) the work done by the compression force and b) the elastic potential energy stored in the bonds and contacts for different joint's roughness.

  44 a) distribution of kinetic energy fluctuations, b) distribution of elastic energy fluctuations, c) correlation function between kinetic energy fluctuations and fracture events (joint DEM S4), d) correlation function between elastic energy fluctuations and fracture events (joint DEM S4), e) correlation function between kinetic energy fluctuations and elastic energy fluctuations (joint DEM S4). The joints DEM S1-S8 where sheared in the -x direction at 14MPa of normal stress.

Figure 4 -Figure 4 -

 44 Figure 4-45 a) evolution of damped energy and b) fractures as function of the shear displacement, c) distribution of damping energy fluctuations, and d) correlation function between fracture events and damping fluctuations for the joint DEM S4.

Figure 4 -

 4 Figure 4-46 Correlation function between the kinetic energy and damped energy fluctuations for the joint DEM S4.

  -affine isotropic with roughness exponent Hx=Hz=0.8, correlation length of 10% the size of the surface Lcx=Lcz=10%, and height variance of σ 2 =120 arbitrary units; b) Reference self-affine isotropic surface with , H=0.8, Lc =10%, and σ 2 = 30 arbitrary units; c) log-log plot of the structure function computed for the surfaces in a) and b), d) variogram computed for the surfaces in a) and b); e) polar representation of the roughness exponent H computed for various directions for the surface in a). -affine anisotropic with roughness exponent Hx=Hz=0.8, correlation length of 10% and 20% the size of the surface Lcz=10%, Lcx=20%, and height variance of σ 2 =30 arbitrary units; b) Reference self-affine isotropic surface with , H=0.8, Lc =10%, and σ 2 = 30 arbitrary units; c) log-log plot of the structure function computed for the surfaces in a) and b), d) variogram computed for the surfaces in a) and b); e) polar representation of the roughness exponent H computed for various directions for the surface in a). -affine anisotropic with roughness exponent Hx=0.6 and Hz = 0.8 , correlation length of 10% Lcx=Lcz=10%, and height variance of σ 2 =30 arbitrary units; b) Reference self-affine isotropic surface with , H=0.8, Lc =10%, and σ 2 = 30 arbitrary units; c) log-log plot of the structure function computed for the surfaces in a) and b), d) variogram computed for the surfaces in a) and b); e) polar representation of the roughness exponent H computed for various directions for the surface in a).EXPERIMENTO 15 units; b) Reference self-affine isotropic surface with , H=0.8, Lc =10%, and σ 2 = 30 arbitrary units; c) log-log plot of the structure function computed for the surfaces in a) and b), d) variogram computed for the surfaces in a) and b); e) polar representation of the roughness exponent H computed for various directions for the surface in a).EXPERIMENTO 15

  is defined by the following five parameters: normal and shear stiffness, and [stress/displacement]; normal and shear strength, and ̅ [stress]; and bond disk radius, . The total force and moment associated with the parallel bond are denoted by and

	A parallel bond
	Bond	(tension)	(shear)	Bond breaks
	breaks
				Slip behavior
				(Only for contacts
				without bonds)
		(	)	1
		1

Table 3 -1 Experimental data Flamand (2000)

 3

	Parameter	Value
	Young Modulus	30 853 MPa
	Poisson's coefficient	0.187
	Compression strength	74.6 MPa
	Tension strength	6.6 MPa

Relative frequency Max. Normal Stress (Pa) Normal Stresses in the REV (Parallel Bonds)

	14%															30%
	12%															25%
	10%															20%
	8%
																15%
	6%
	4%															10%
	2%															5%
	0%															0%
	-8.0E-06	-2.8E-06	2.5E-06	7.7E-06	1.3E-05	1.8E-05	2.3E-05	2.9E-05	3.4E-05	3.9E-05	4.4E-05	5.0E-05	5.5E-05	6.0E-05	6.5E-05	7.0E-05	0.0E+00	8.8E-07	1.8E-06	2.7E-06	3.5E-06	4.4E-06	5.3E-06	6.2E-06	7.1E-06	8.0E-06	8.8E-06	9.7E-06	1.1E-05	1.1E-05	1.2E-05	1.3E-05	1.4E-05

Relative frequency Max. Shear Stress (Pa) Shear Stresses in the REV (Parallel Bonds)

Table 3 -2a Experimental roughness parameters (Marache, 2002

 3

			, after Flamand (2000)) used to
		generate the surface S5
	Roughness exponent H		0.56
	Correlation length L c		0.16 cm
	Variance of height σ²		0.55 mm²
	Table 3-2 sets of parameters used to generate the self-affine surfaces
	Surface name	Roughness exponent	Correlation length	Height variance
		(H)	(L C )	(σ²)
	S1	H H = 0.75	L C-L = 1.6 cm	σ² H =0.547 mm²
	S2	H H = 0.75	L C-L = 1.6 cm	σ² L =0.328 mm²
	S3	H H = 0.75	L C-H = 2.0 cm	σ² H =0.547 mm²
	S4	H H = 0.75	L C-H = 2.0 cm	σ² L =0.328 mm²
	S5	H L = 0.52	L C-L = 1.6 cm	σ² H =0.547 mm²
	S6	H L = 0.52	L C-L = 1.6 cm	σ² L =0.328 mm²
	S7	H L = 0.52	L C-H = 2.0 cm	σ² H =0.547 mm²
	S8	H L = 0.52	L C-H = 2.0 cm	σ² L =0.328 mm²
	*Sub index H and L stand for high and low value of the parameter.

Table 3 -3b sets of parameters of the self-affine joints created with DEM Surface name Roughness exponent (H) Correlation length (L C )

 3

	Height variance
	(σ²)

Table 3 -2c Comparison of roughness parameters obtained experimentally (Marache, 2002) and roughness parameters of the numerical joint DEM S1

 3

	Experimental	DEM S1

Shear displacement (mm) Influence of geometry in shear (Rigid model)

		40								2.50
	Shear stress (MPa)	10 15 20 25 30 35			σ = 14 MPa σ = 14 MPa σ = 14 MPa		DEM S1 DEM S2 DEM S3 DEM S4 DEM S5 DEM S6 DEM S7 DEM S8	0.50 2.00 1.00 1.50 Dilatancy (mm)
		5
		0								0.00
		0.00	1.00	2.00	3.00	4.00	5.00	6.00	7.00	0.00	1.00	2.00	3.00	4.00	5.00	6.00	7.00

Shear displacement (mm) Influence of geometry in shear (Rigid model)

									DEM S1 DEM S2	σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa
									DEM S3
									DEM S4
									DEM S5
									DEM S6
									DEM S7
									DEM S8
		40
		35
	Shear stress (	10 15 20 25 30
		5
		0
		0.00	1.00	2.00	3.00	4.00	5.00	6.00	7.00

Mpa) Shear displacement (mm) Influence of geometry in shear (Rigid model)

			2.50
	σ = 21 MPa σ = 14 MPa σ = 21 MPa	DEM S1
		DEM S2
		DEM S3 DEM S4	2.00
		DEM S5 DEM S6 DEM S7 DEM S8	1.00 1.50 Dilatancy (mm)
			0.50
			0.00
			0.00	1.00	2.00	3.00	4.00	5.00	6.00	7.00

Shear displacement (mm) Influence of geometry in shear (Rigid model)

	DEM S1 DEM S2	σ = 21 MPa σ = 14 MPa σ = 21 MPa
	DEM S3
	DEM S4
	DEM S5
	DEM S6
	DEM S7
	DEM S8

Table 4 -1 Comparison between rigid and elastic model results ( N =14 MPa) DEM Surface Shear Peak Stress (MPa) Residual Shear Stress (MPa) Initial stiffness (MPa/mm)

 4

		Rigid	Elastic	Rigid	Elastic	Rigid	Elastic
	S1	16.15	18.05	10.24	12.48	56.90	82.83
	S2	13.87	15.24	8.33	10.93	58.13	96.50
	S3	14.25	15.71	9.94	13.10	59.80	97.40
	S4	12.32	13.06	9.07	11.30	53.20	105.30
	S5	20.51	22.34	9.00	12.63	64.10	82.90
	S6	16.67	18.37	7.73	11.45	63.60	83.50
	S7	18.23	20.28	7.96	12.23	68.20	108.90
	S8	15.34	16.53	8.74	10.54	62.50	88.60

Shear displacement (mm) Influence of geometry in shear (Rigid model)

									Influence of geometry in shear (Elastic model)
												σ = 14 MPa σ = 14 MPa		DEM S1 DEM S2
															DEM S3
								Shear stress (MPa)							DEM S4 DEM S5 DEM S6 DEM S7 DEM S8
									0.00	1.00	2.00	3.00	4.00	5.00	6.00	7.00
											Shear displacement (mm)
	2.50	DEM S2 DEM S1		σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa σ = 14 MPa				2.50
	2.00	DEM S3 DEM S4							2.00
	1.00 1.50 Dilatancy (mm)	DEM S5 DEM S6 DEM S7 DEM S8						Dilatancy (mm)	1.00 1.50
	0.50								0.50
	0.00								0.00
	0.00	1.00	2.00	3.00	4.00	5.00	6.00	7.00	0.00	1.00	2.00	3.00	4.00	5.00	6.00	7.00

Shear displacement (mm) Influence of geometry in shear (Elastic model)

	DEM S1 DEM S2	σ = 14 MPa σ = 14 MPa
	DEM S3
	DEM S4
	DEM S5
	DEM S6
	DEM S7
	DEM S8

Influence of geometry in shear (Elastic-fracture model)

		40					2.50
				σ = 14 MPa σ = 14 MPa σ = 14 MPa	DEM S1
		35			DEM S2
	Shear stress (MPa)	10 15 20 25 30			DEM S3 DEM S4 DEM S5 DEM S6 DEM S7 DEM S8	Dilatancy (mm)	2.00 1.00 1.50
							0.50
		5
		0					0.00
		0.00	1.00	2.00 Shear displacement (mm) 3.00 4.00 5.00	6.00	7.00	0.00	1.00	2.00 Shear displacement (mm) 3.00 4.00 5.00	6.00	7.00

Influence of geometry in shear (Elastic-fracture model)

													DEM S1		σ = 14 MPa σ = 14 MPa σ = 14 MPa
													DEM S2
													DEM S3
													DEM S4
													DEM S5
													DEM S6
													DEM S7
													DEM S8
		Influence of geometry in shear (Elastic-Fracture model)			Influence of geometry in shear (Elastic-Fracture model)
	Shear stress (MPa)	10 15 20 25 30 35 40			σ = 14 MPa σ = 21 MPa		DEM S5 DEM S6 DEM S7 DEM S8 DEM S3 DEM S4 DEM S1 DEM S2	Dilatancy (mm)	1.00 1.50 2.00 2.50	DEM S5 DEM S6 DEM S7 DEM S8 DEM S4 DEM S3 DEM S1 DEM S2		σ = 14 MPa σ = 21 MPa
												0.50
		5
		0										0.00
		0.00	1.00	2.00	3.00	4.00	5.00	6.00	7.00		0.00	1.00	2.00	3.00	4.00	5.00	6.00	7.00
				Shear displacement (mm)						Shear displacement (mm)
							40							0.40%
						Shear stress (MPa)	5 10 15 20 25 30 35	σ = 21 Mpa σ = 14 Mpa					0.05% 0.10% 0.15% 0.20% 0.25% 0.30% 0.35%	Cumulative number of fractures/
							0							0.00%
							0.00		2.00	4.00			6.00

Total number of bonds of the joint Influence of compression in fracturing (DEM S7) Shear displacement (mm)

Table 4 -3 Roughness parameters obtained from the experimental joints and simulated DEM joint (DEM S1). Experimental Flamand (2000), Marache (2002)

 4

		DEM
	Height variance = 0.6 mm²	Height variance = 0.55 mm²
	Correlation length= 15 mm	Correlation length= 15 mm
	Roughness exponent = 0.55	Roughness exponent = 0.53

Table 4 -4 Comparison of the results obtained experimentally Flamand (2000) and by DEM simulation (DEM S1).

 4

	Compression stress	Peak shear stress	Residual shear stress	Shear stiffness	Shear peak displacement
	(MPa)		(MPa)	(MPa)		(MPa/mm)	(mm)
	Exp	DEM	Exp	DEM	Exp	DEM	Exp	DEM	Exp	DEM
	14	14	18.37	17.92	11.52	11.68	51	79	0.502	0.72
	21	21	24.25	25.66	17.00	15.86	68	110	0.599	0.66
	*Exp = values obtained experimentally by

δE elastic δE elastic (Joules) Histogram of energy of elastic waves in shear for different morphologies σ = 14 MPa

	10000	σ = 14 MPa			DEM S1 DEM S2				10000						DEM S1 DEM S2
	y = 3E-27x -2.587 R= 0.978 1E-11 δE kinetic (Joules) y = 2E-24x -2.387 R= 0.993 1E-10 11.5 17.3 23.1 28.9 34.6 40.4 46.2 52.0 57.7 63.5 Correlation between kinetic fluctuations and fracture 1E-12 1 10 100 1000 DEM S3 1E-09 DEM S4 DEM S5 DEM S6 DEM S7 DEM S8 -0.80 Frequency δE kinetic -0.60 -0.40 -0.20 0.00 0.20 0.40 0.60 0.80 0.0 5.8 69.3 75.1 80.8 86.6 92.4 Correlation coefficient Frequency -1.00 -1.00 1000 0.1 1 10 100 1E-11 y = 2E-39x -3.998 R= 0.983 0.0 5.8 11.5 17.3 23.1 -0.80 -0.60 -0.40 -0.20 0.00 0.20 0.40 0.60 0.80 Correlation coefficient Correlation between elastic fluctuations and fracture DEM S3 DEM S4 y = 9E-36x -3.817 DEM S5 R= 0.982 DEM S6 1E-10 1E-09 1E-08 DEM S7 DEM S8 28.9 34.6 40.4 46.2 52.0 57.7 63.5 69.3 75.1 80.8 86.6 92.4 1.00 fluctuations 1.00 fluctuations
		Lag distance between signals (µm)									Lag distance between signals (µm)
		1.00
		0.80
		0.60
		0.40
		0.20
		0.00
		-0.20
		-0.40
		-0.60
		-0.80
		-1.00
		0.0	5.8	11.5	17.3	23.1	28.9	34.6	40.4	46.2	52.0	57.7	63.5	69.3	75.1	80.8	86.6	92.4

Cumulative fractures Shear displacement Influence of geometry on fractures (broken bonds)

			Evolution damping energy in shsear
	Damping energy (Joules)	4.0E-08 6.0E-08 8.0E-08 1.0E-07 1.2E-07 1.4E-07 1.6E-07 1.8E-07 2.0E-07	DEM S6 DEM S7 DEM S8 DEM S4 DEM S5 DEM S1 DEM S2 DEM S3	σ= 14 MPa		2,000 3,000 4,000 5,000 6,000 7,000 9,000 8,000	DEM S6 DEM S7 DEM S8 DEM S5 DEM S4 DEM S3 DEM S1 DEM S2				σ = 14 MPa
		2.0E-08
	0.0E+00				1,000
		0.00	2.00 Shear displacement (mm) 4.00	6.00	0 0.00		1.00		2.00	3.00		4.00		5.00		6.00	7.00
	10000
				σ= 14 MPa	DEM S1	1.00
	Frequency δE damping	10 100 1000	y = 5E-31x -3.214 R = 0.987	R = 0.995 y = 6E-25x -2.716	DEM S6 DEM S7 DEM S8 DEM S4 DEM S5 DEM S3 DEM S2	-0.40 -0.20 0.00 0.20 0.40 0.60 0.80
		1				-0.60
						-0.80
		0.1				-1.00
		1E-11	1E-10 δE damping (joules) 1E-09	1E-08	0.0	5.8	11.5	17.3	23.1	28.9	34.6	40.4	46.2	52.0	57.7	63.5	69.3	75.1	80.8	86.6	92.4

Remerciements

Elastic-fracture model

In this final section, the elastic-fracture model defined in Chapter 3 is used in the shear simulations of the joints DEM S1 to DEM S8. As for the rigid and elastic models, two different compression stresses are applied to the elastic-fracture model (14 MPa and 21 MPa), and a shear displacement is imposed in two perpendicular directions -x and -y. During the simulations, the shear stress, dilatancy, broken bonds, kinetic energy and elastic potential energy of the model are recorded. Figure 4-27 and Figure 4-28 show the results obtained for the shear simulations of the surfaces DEM S1-S8 when a compression stress of 14 MPa and 21 MPa is applied, respectively.

The solid lines in Figure 4-27 and Figure 4-28 represent the values of shear stress, the dilatancy and the friction related to particles located out of the expected rough surface obtained for the -x direction of shear, the dotted lines represent these parameters when the model is sheared in the -y direction. Note that, as for the rigid and elastic models, the contribution in the friction phenomenon of the particles located out of the expected rough surfaces (named "artifact shear stress" and represented by the solid and dashed green lines in Figure 4-27 and Figure 4-28) is completely negligible compared to the total shear stress in these simulations.

Table 4-1 resumes the average shear peak stresses, residual shear stress and initial shear stiffness for all the surfaces S1-S8 for the rigid, elastic, and elastic-fracture models for a normal stress equal to 14 MPa. Note from Table 4-2 that the use of the elastic model for the friction of the DEM joints, leads to an increase of the shear peak stresses and the residual shear stresses compared to those obtained with the use of the rigid joint model. The values obtained when the elastic-fracture model is implemented, lie between the elastic and the rigid ones. The same argument applies for the initial shear stiffness.

If we assume that each fracture of bonds leads to the same amount of elastic energy converted in the form of elastic (acoustic) waves then the distribution of fractures plotted in Figure 4-37a can represent an approximate picture of the distribution of the amplitude of the acoustic emissions associated to the micro-fractures. For comparison, Figure 4-37b shows the distribution of the acoustic emission intensity obtained from the fracturing of concrete samples obtained from [START_REF] Saliba | Apport de l'émission acoustique dans la compréhension et la modélisation du couplage fluage-endommagement du béton[END_REF].

It is possible to see from Figure 4-37 that the probability of the acoustic events obtained from the shear simulations are very similar to the experimental acoustic emission obtained by [START_REF] Saliba | Apport de l'émission acoustique dans la compréhension et la modélisation du couplage fluage-endommagement du béton[END_REF] despite the boundary conditions of the experiments and the simulations are different (direct shear test at constant normal load, and, uni-axial compression respectively).

Another interesting issue is that the curves distribution of fracture events obtained from the shear simulations of the surfaces DEM S1-S8 tends to collapse in one single curve . A detailed study of these curves indicates that the behavior of the distribution of fracture events as function of the joint roughness can be divided roughly in two parts separated by the vertical dotted line in Figure 4-38. Indeed, they are superimposed when the magnitude of the fracture events is low (1-3 fractures), while as the magnitude of the events increases (more than 4 fractures), the roughest surface DEM S5 tends to produce fracture events of major magnitude than the smoothness surfaces DEM S4. Again, all the curves are ranged between the boundary responses corresponding to the roughest and the smoothest surfaces (DEM S5 and DEM S4 respectively). Thus, for small magnitudes of fracture events, the distribution of the events is independent of the surface roughness while, for bigger magnitudes of events, the magnitude of the fracture events increases as a function of the surface roughness. A possible explanation of this behavior is that, an increase of the roughness could involve (1) an increase of the asperities impeding the shear displacement and consequently (2) an increase of the fractures as the shear progresses.

Note that such behavior is also observed for the earthquakes and is expressed by the Gutenberg-Richter law:

where N is the number of earthquakes with magnitude equal or greater than M, and a and b are constants determined by the seismic activity of the region of interest. In equation 4-3, M is a measure of the amount of energy released during the earthquakes.

Friction of rock joints from an energetic point of view

The shear phenomena can also be studied from an energetic point of view.