REMERCIEMENTS

Le travail de thèse présenté dans ce mémoire a été effectué au sein de l'équipe géotechnique (CERMES) du laboratoire Navier, et financé par une allocation de recherche de l' École Normale Supérieure de Cachan.

Je tiens à remercier chaleureusement mes deux directeurs de thèse, Behrouz Gatmiri et Ahmad Pouya. Behrouz a été à l'origine de cette thèse et je le remercie de m'avoir fait confiance et de m'avoir laissé une grande liberté dans mes choix scientifiques. Je remercie Ahmad pour ses conseils avisés et sa présence lors de mes derniers mois de rédaction.

Je tiens également à remercier les membres de mon jury de thèse, Djimédo Kondo et Farid Benboudjema pour avoir accepté de rapporter ma thèse, et Antonio Gens pour être venu de Barcelone participer à ma soutenance.

Je remercie bien sûr vivement Chloé Arson pour son soutien et son aide scientifique tout au long de ma thèse et plus particulièrement lors de mon séjour à GeorgiaTech. J'espère que ce ne sont que les débuts d'une longue collaboration. J'ai également beaucoup apprécié l'environnement de travail au CERMES. Merci à tous les chercheurs, doctorants et post-doctorants pour l'ambiance chaleureuse et les interactions scientifiques constructives. Je tiens à remercier plus particulièrement Jean-Michel Pereira pour ses conseils scientifiques et son aide lors de ma recherche d'emploi.

La thèse, ce n'est pas tous les jours facile. Je tiens donc à remercier mes collègues doctorants Valérie, Jean-François, Simona et Vinh, avec qui j'ai partagé des bureaux, des inquiétudes et des moments de détente. C'était un plaisir de faire ce bout de chemin ensemble.

Merci enfin à mes amis et à ma famille pour les bons moments passés en dehors du travail, et bien sûr à Antoine, de m'avoir soutenu et d'avoir toujours cru en moi. i This thesis work objective is the development of a hydro-mechanical constitutive model which accounts for damage-plasticity couplings in unsaturated geomaterials. The coupling of this model with hydraulic transfers is made possible by its implementation into the Finite Element code Θ-Stock. In order to achieve this implementation, a specific stresspoint algorithm has been developed. Fully coupled hydro-mechanical problems have been simulated, such as the creation of the Excavation Damaged Zone around a tunnel and the initiation of damage due to desiccation and humidification.

A double effective stress incorporating both the effect of suction and damage (assumed isotropic) is defined based on thermodynamical considerations. The advantage of this approach is that it results in a unique stress variable being thermodynamically conjugated to elastic strains. A pressure-dependent hyperelastic formulation is used to describe the behaviour inside the elastic domain. The evolution of elastic rigidities with damage is then studied. Two hypotheses are compared, the principle of strain equivalence and the principle of equivalent energy.

, itself based on the Cam-Clay model. Two distinct criteria are defined for damage and plasticity, which can be activated either independently or simultaneously. Their formulation in terms of effective stress and suction allows them to evolve in the total stress space with suction and damage changes. This leads to a direct coupling between damage and plasticity and allows the model to capture the ductile/brittle behaviour transition occurring when clays are drying.

A specific explicit algorithm has been developed to handle the association of the two dissipative phenomena. The implementation of the constitutive model into the Finite Element code Θ-Stock allowed for the simulation of fully coupled hydro-mechanical problems. The hydraulic transfer laws also consider the saturation state.

This fully coupled model is first applied to simulate the development of micro-cracks during dessication of a soil sample. Damage initiation is explained mainly by the important iii iv

The role of a mathematical model for the description of soil behaviour is twofold: firstly, as a simple qualitative framework against which soil behaviour may be assessed and, secondly, as a complete quantitative model for analysis and design. In its first role the model is of use both in judging the quality and consistency of data, and in predicting the character and trends of behaviour of a soil under a variety of circumstances. For this purpose the model must represent the essential features of soil behaviour in as simple a manner as possible.

A simple model is also desirable for the second role; this is because, in a detailed numerical calculation, an excessively complex model may be difficult to implement and expensive to use. Very often the quality of the available data about the actual soil from the field will not warrant the additional precision of calculation achieved using a complex model. If, however, the model is to be used for mathematical calculations, it is important that it should be properly founded in the theory of continuum mechanics and should be internally consistent.

A simple model will require few parameters for the description of a soil, and each of these should have a definite physical significance so that an assessment can easily be made of its importance and the likely results of any changes in its value. The parameters should be measurable

INTRODUCTION

Geotechnical engineering's primary purpose has been to provide methods to engineers so as to predict how the ground would react to man-made constructions, with a focus on settlements (small-strain behaviour) and rupture.

With the rise of ecological and energy concerns, research in geotechnical engineering now goes beyond soil-structure problems, and has been extended to a wide range of soil-environment interactions issues [START_REF] Gens | Soil-environment interactions in geotechnical engineering[END_REF]. Among the various environmental components which have an impact on soil, water plays a major role. Water invades or leaves the ground during rainfall and evaporation, sometimes provoking damage such as landslides or desiccation cracks. Underground flows also have their importance when digging deep galleries. Nearing the contact with atmosphere, soil often is in an unsaturated state, which affects the material mechanical properties as well as the fluid transfer laws. Great progress has been made in the last decades in a view to provide modelling frameworks capable of reproducing soil behaviour in all their states, saturated or not. This research field is still very active, with many more issues to cope with [START_REF] Sheng | Unanswered questions in unsaturated soil mechanics[END_REF].

In parallel, degradation of geomaterials has also received much interest. Soils and rocks properties, such as rigidity, strength, or permeability, are known to be subject to changes after being submitted to hydric or mechanical solicitations. Several approaches have been used to model this degradation. Some works assume the elastic properties to evolve with plastic strains [START_REF] Sulem | Elasto-plastic modelling of Red Wildmoor sandstone[END_REF][START_REF] Kavvadas | A constitutive model for structured soils[END_REF][START_REF] Rouainia | A kinematic hardening constitutive model for natural clays with loss of structure[END_REF][START_REF] Baudet | A constitutive model for structured clays[END_REF]. This type of model is sometimes called structured model. Another approach, the one that is used in this thesis, is to use the framework of Continuum Damage Mechanics, first developed for metals and later extended to concrete and rocks. This approach assumes that the degradation of material properties is due to the initiation and propagation of microcracks. Contrary to the previous approach, it allows damage and plastic strains to be developed independently. The main modelling difference is the postulate of an effective stress, applied on the intact solid matrix, which controls the global material behaviour. Concerning granular materials, a new family of models consider the crushing, or breakage, of grains submitted to compressive loads [START_REF] Einav | Breakage mechanics-Part I: Theory[END_REF].

Up to now, few attempts have been made to model damage in unsaturated geomaterials. Some models have been developed which consider damage in unsaturated geomaterials [START_REF] Arson | A mixed damage model for unsaturated porous media[END_REF], damage-plasticity couplings in saturated geomaterials [START_REF] Chiarelli | Modeling of elastoplastic damage behavior of a claystone[END_REF][START_REF] Conil | Poroplastic damage model for claystones[END_REF] and even damage-plasticity couplings in unsaturated geomaterials [START_REF] Hoxha | Saturated and unsaturated behaviour modelling of Meuse-Haute/Marne argillite[END_REF][START_REF] Jia | Elastoplastic damage modelling of argillite in partially saturated condition and application[END_REF]. However, these models, initially formulated for rock, ignore some specific important features of clay soil behaviour, such as the dependence of elastic moduli to pressure. Moreover, damage-plasticity models proposed for rock so far fail at predicting the transition between ductile and brittle behaviour associated with suction increase.

The proposed work aims to develop a framework to couple damage and plasticity in unsaturated porous media, for which clay minerals in the solid matrix are expected to play an important role.

Two approaches have to be considered when constructing a modelling framework. One can either make assumptions about the nature of phenomena taking place at the microscale, develop micromechanical behaviour laws and upscale them to obtain the macroscopic behaviour. This is the micromechanical approach. The other strategy is to define macroscopic laws, which aim at reproducing some of the observed phenomena. These laws depend on some state variables, chosen for their sound physical meaning, and are developed in order to respect thermodynamic principles. This is the phenomenological approach, which is mainly based on energy considerations, and which is used in the present work.

An accurate prediction of all aspects of a material behaviour leads to very complex models requiring a great number of parameters, which makes calibration difficult. Since this work mainly aims to develop a flexible constitutive framework and to highlight its main characteristics concerning damage-plasticity couplings and non saturation, we chose not to account for certain features, such as intrinsic and damage-induced anisotropy of the material.

Chapters 1 to 3 present the mechanical constitutive model assumptions and equations. They develop how damage and suction enter in the definition of the double effective stress, which is the stress considered to be applied on the solid skeleton and to drive the overall behaviour (chapter 1). It is then shown how this double effective stress, coupled with a hyperelastic pressure dependent formulation, results in a dependence of apparent rigidities on suction and damage (chapter 2). Finally, damage and plasticity criteria are presented.

Their expression in terms of the double effective stress results in an automatic dependence on suction and damage when compared to total stresses (chapter 3).

Chapter 4 develops the air and water transfer laws needed to simulate fully coupled hydromechanical problems. These transfer laws depend on the saturation state, which makes the problem non-linear, not only for the mechanical part, but also for the hydraulic part.

In chapter 5 the strong and weak forms of the partial differential equation system are developed. Spatial and time discretisations are presented, as well as the resulting matricial form of the problem. The global algorithm, used to solve the boundary value problem, consists in dividing the load in small load increments, and to perform iterations, using the Newton-Raphson algorithm, in order to find the displacement and pore pressure fields respecting equilibrium. A specific local algorithm has been developed for the present CHAPTER 1

DEFINITION OF THE DOUBLE EFFECTIVE STRESS

This chapter aims at introducing a double effective stress incorporating both suction and damage effects. The formulation of this stress is based on thermodynamical considerations. Elastic strains will be directly related to this double effective stress through elastic potentials.

A strong emphasis is put on the soil water retention properties which have a great impact on this double effective stress.

Ce chapitre a pour but d'introduire une contrainte doublement effective incorporant les effets aussi bien de la succion que de l'endommagement. Sa formulation est basée sur des considérations thermodynamiques. Cette contrainte doublement effective sera directement liée aux déformations élastiques à travers des potentiels élastiques. Une attention particulière est apportée aux propriétés de rétention du sol, celles-ci ayant une influence importante sur la contrainte doublement effective.

Unsaturated soils

This work deals with unsaturated geomaterials. Unsaturated soils have first been treated as specific materials different from the saturated ones. However, it is clear nowadays that, although suction has a great impact on unsaturated porous media behaviour, non saturation should only be considered as a specific state of the material.

Unsaturated soils are multiphasic. They are constituted of a solid phase and two fluid phases, the water liquid phase and the air gaseous phase (see figure 1.1). More complex models can consider dissolved air into the liquid phase and water vapour into the gaseous phase.

Water meniscus Air Fig. 1.1: Schematic view of an unsaturated medium

The total volume, V tot , is the sum of the volume occupied by the solid phase, V s , and the volume of the pores, V v , itself divided between a volume filled by water, V w , and a volume filled by gas, V g (see figure 1.2). The relative proportions of these phases can be described by the porosity, φ, and the saturation degree, S l :

Solid

Liquid Gas

φ = V v V tot (1.1) S l = V w V v (1.2)
When water evaporates from the saturated soil, a difference between the air and water pressures arises. This pressure difference, s = p gp l , is called suction. In the case in which air remains equal to the atmospheric pressure, water pressure is negative and suction takes a positive value. The air-water interface (called meniscus) starts to curve when suction increases. The radius of the meniscus decreases when suction increases, and once it becomes as small as the pore throats, air can invade the porous structure, which becomes unsaturated. Fine soils can be submitted to a high suction level before starting to desaturate. The suction for which air breaks through the soil mass is called the air entry value.

The combination of the water surface tension, F t , and the negative pore water pressure, p l , results in a force that tends to pull the soil grains towards one another. The resulting force, F c , on the solid skeleton is similar to a compressive stress (see figure 1.3).

Fig. 1.3: Air-water-solid interaction for two spherical particles An increase in suction will therefore lead to a decrease of the total volume (shrinkage), and wetting soils (i.e. decreasing suction) will usually make them swell. Suction also contributes to stiffen the soil against external loading thanks to the bounding between grains created by water menisci in tension. The additional component of normal force at the contact will also prevent slippage between grains and thus increase the external force needed to cause plastic strains. However, when wetting a soil under constant mechanical loading, the resaturation destroys the bounds due to water menisci and may induce an irrecoverable volumetric compression (called collapse). These main characteristics of unsaturated soils mechanical behaviour are represented on figure 1.4. Changes in suction may also induce irreversible processes (plasticity or damage) during a drying process.

Suction is closely related to the degree of saturation through the water retention curve (see section 1.1.3).

Unsaturated soil models are usually extensions of saturated soil ones. The most widely used of them is the Cam-clay model, first developed by [START_REF] Roscoe | On the yielding of soils[END_REF] and later modified by [START_REF] Roscoe | On the generalized stress-strain behaviour of 'wet' clay[END_REF]. However, the debate still goes on about which stress framework to use to extend saturated soil models to account for unsaturated states.

Modelling approaches in the literature

A comprehensive review of the existing stress frameworks can be found in the paper of [START_REF] Nuth | Effective stress concept in unsaturated soils: Clarification and validation of a unified framework[END_REF].

In the following we use the convention of soil mechanics which means that compressive stresses are positive.

Bishop's effective stress [START_REF] Bishop | The effective stress principle[END_REF] attempted to extend Terzaghi's effective stress principle for saturated soils to unsaturated soils. He proposed the following form for the effective stress:

σ * = σ -[(1 -χ)p g + χp l ]I (1.3)
in which χ is a parameter weighting the contributions of liquid and gas pressures on the effective stress. The relative effect of each phase on the overall behaviour depends on the saturation state of the material. Indeed, at very low saturation states, even if suction is high, the contact surface between water and grains is so small that the resulting effect will be negligible. χ is therefore often considered as being a function of the degree of saturation, ranging from 0 for dry states to 1 for saturated states.

The stress defined in equation (1.3), sometimes called Bishop's stress, generalised effective stress, or constitutive stress, can also be expressed in terms of net stress and suction:

σ * = σ -[p g -χ(p g -p l )]I = σ net + χsI (1.4)
Early experimental and theoretical works have attempted to find a unique relationship between the χ parameter and the degree of saturation. Some experimental data are shown in figure 1.5. Table 1.1: Formulas used in the literature for χ = χ(S l )

The difficulty of finding a unique relationship led to some doubts in the scientific community about the validity of Bishop's effective stress approach. Moreover, using Bishop's stress alone within an elastic framework doesn't allow the collapse phenomenon described previously to be captured (figure 1.4). Indeed, decreasing suction induces a decrease in the effective stress, which can not be associated to opposite volume changes depending on the external applied load (swelling for low stress and collapse for high stress). This appeared to be the main limitation to the use of an effective stress for partially saturated soils at that time (until it was incorporated into an elasto-plastic framework).

Independent stress variables

As a response to the limitations of the single effective stress approach, [START_REF] Coleman | Stress/strain relations for partly saturated soil[END_REF] was the first to introduce the idea of using two independent variables to describe unsaturated soils behaviour. [START_REF] Fredlund | Stress state variables for unsaturated soils[END_REF] proposed to use one of the following combinations of the total stress σ, the water pressure p l and the air pressure p g :

(σ -p g I) and (p g -p l ) (1.5) (σ -p l I) and (p g -p l ) (1.6) (σ -p g I) and (σ -p l I) (1.7)
The first of these combinations is the most frequently used. However, the use of net stress, σ net = σp g I, and suction to model unsaturated soils behaviour fails to provide a smooth transition between unsaturated and saturated states when suction becomes zero. Indeed, if the constitutive model is expressed in terms of net stress and suction, as soon as suction becomes equal to zero, the net stress remains the only variable controlling the soil behaviour, which means that elastic strains and shear strength will depend on net stress only, and that the effect of a positive water pressure is not accounted for. It would be of course possible to switch to a saturated soil model when suction reaches zero, but it would not be efficient in terms of numerical implementation, when one wants to be able to consider positive water pressures as well as suction.

Bishop's effective stress combined with an extra variable

Due to the limitations expressed above, recent works have been focusing on combining an effective stress following the definition of Bishop with an extra variable, usually linked to suction. Collapse can then be accounted for by using an elasto-plastic framework.

As it is shown in the next section, the effective stress and suction related variables can be derived from thermodynamic principles.

Constitutive stress based on energy considerations

References about the thermodynamics of open porous media can be found in [START_REF] Coussy | [END_REF], [START_REF] Collins | Application of Thermomechanical Principles to the Modelling of Geotechnical Materials[END_REF], [START_REF] Hansen | A thermodynamically consistent framework for theories of elastoplasticity coupled with damage[END_REF] and [START_REF] Coussy | Approche énergétique du comportement des sols non saturés[END_REF].

The energy balance can be written as (chemical or electrical energies are neglected):

K + U = P ext + Q (1.8)
with:

❼ K, the kinetic energy ❼ U , the specific internal energy ❼ P ext , the power of external force ❼ Q, the heat flow into the system

The kinetic energy theorem gives [START_REF] Coussy | Mécanique des milieux poreux[END_REF]:

K + P def + P f = P ext (1.9)
with ❼ P def = Ẇ , the strain work rate ❼ P f , the power of inertia forces (relative movement of fluids to solid skeleton)

So we have the variation of internal energy

U = Ẇ + Q + P f (1.10)
which, expressed locally becomes:

u = ẇ -∇ • q + p f (1.11)
The local strain work rate is given by (deformation of the solid skeleton):

ẇ = σ : ε + p l φl + p g φg (1.12)
in which φ l is the relative liquid porosity:

φ l = φS l (1.13) φl = φS l + φ Ṡl (1.14)
and φ g is the relative gas porosity.

φ g = φ -φS l (1.15) φg = φ -φS l + φ Ṡl (1.16)
Equation 1.12 can be rewritten as

ẇ = σ : ε + p l S l φ + p l φ Ṡl + p g φ -p g S l φ -p g φ Ṡl = σ : ε + (S l p l + (1 -S l )p g ) φ + (p l -p g )φ Ṡl (1.17)
If we assume the incompressibility of the solid skeleton, the porosity change is directly related to the volumetric strain, φ = -I : ε, which gives:

ẇ = [σ -(S l p l + (1 -S l )p g )I] : ε -(p g -p l )φ Ṡl = σ * : ε -s * Ṡl (1.18) with σ * = σ -(S l p l + (1 -S l )p g )I = σ -p g I + (p g -p l )S l I = σ net + sS l I s * = (p g -p l )φ = sφ (1.19)
σ * can be called the average skeleton stress tensor [START_REF] Jommi | Remarks on the constitutive modelling of unsaturated soils[END_REF], the constitutive stress [START_REF] Sheng | Finite element formulation and algorithms for unsaturated soils. Part I: Theory[END_REF] or the generalised effective stress [START_REF] Laloui | On the use of the generalised effective stress in the constitutive modelling of unsaturated soils[END_REF].

s * is called the modified suction.

Constitutive equations can be derived from Clausius-Duhem Inequality (see part 3, equations 3.13 and 3.15 ):

σ * = ∂ψ ∂ε e
(1.20)

s * = - ∂ψ ∂S l (1.21)
in which ψ is the Helmholtz free energy.

A possible extension of the model could be to consider a compressible solid skeleton and to take the biot coefficient b = 1.

We would then have:

σ * = [σ -b(S l p l + (1 -S l )p g )I] (1.22)
Equation 1.19 is identical to one of the combinations given by [START_REF] Houlsby | The work input to an unsaturated granular material[END_REF]. We will therefore use from now the following set of state variables:

σ * = σ net + sS l I (1.23) s * = φs (1.24)
When the material becomes saturated, i.e. when S l = 1, equation 1.23 becomes

σ * = σ -p l I = σ ′ (1.25)
which is the saturated effective stress of Terzaghi.

With this stress choice, it is easy to model the behaviour of a soil submitted to negative as well as positive water pressures, without any discontinuity between these two cases. This is a real advantage for numerical implementation.

1.1.3 Water retention properties and their impact on constitutive stress

The wetting-drying process in a single pore is illustrated in figure 1.6. During drying, suction increases and the meniscus radius decreases until it becomes as small as the pore throat and air is able to invade the pore. This process is not reversible since a small reduction of suction will not induce an immediate filling of the pore by water. During resaturation, water will entirely refill the pore when the meniscus radius will coincide with the largest pore dimension. [START_REF] Miller | Physical Theory for Capillary Flow Phenomena[END_REF])

A typical drying and wetting curve resulting from this hysteretic process is shown in figure 1.7.

The air-entry value of the soil is the suction for which air starts to enter the largest pores in the soil. The residual saturation degree is the saturation degree for which a large suction change is required to remove additional water from the soil. The adsorption curve differs from the desorption curve as a result of the hysteretic process explained before.

Water retention curve basic formulas

The effective degree of saturation is defined as

S e = S l -S r 1 -S r (1.26)
with S r the residual degree of saturation.

For the sake of simplicity, hysteresis effects will be neglected in the following study.

In this case the degree of saturation can be expressed as a bijective function of suction.

S e = f (s) (1.27)
The two most commonly used models are the one of [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF],

S e = s e s λ (1.28)
in which s e is the air entry suction value, and the one of [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF].

S e = 1 1 + (α vg s) nvg mvg (1.29)
for which the number of parameters is sometimes reduced by stating

n vg = 1 1 -m vg (1.30)
The formulation of Brooks and Corey as the advantage of needing less parameters and being mathematically simpler. However, the Van Genuchten formulation usually fits better to experimental data and we will therefore use it in the next developments.

When the porosity changes or when damage occurs, this relationship changes due to the change of pores size and connectivity. Moreover, to respect the thermodynamic framework exposed in the previous section, the degree of saturation should be a function of the modified suction. The following section will therefore deal with the incorporation of porosity into equation 1.29.

Effect of porosity on the water retention curve

Water retention curves (WRC) are usually obtained under constant stress. Strains due to drying or wetting are therefore already taken into account. However, larger changes in porosity due changes in the stress state can notably modify the water retention properties since it depends on the pore size distribution. Several recent works have studied the dependency of water retention curves on porosity and proposed ways to incorporate porosity into WRC formulas (Gallipoli et al., 2003a;[START_REF] Sun | Elastoplastic modelling of hydraulic and stress-strain behaviour of unsaturated soils[END_REF][START_REF] Tarantino | A water retention model for deformable soils[END_REF]Masin, 2010;[START_REF] Sheng | Coupling hydraulic with mechanical models for unsaturated soils[END_REF][START_REF] Buscarnera | The yielding of brittle unsaturated granular soils[END_REF][START_REF] Le | Stochastic analysis of unsaturated seepage through randomly heterogeneous earth embankments[END_REF].

Our primary goal is to respect thermodynamical consistency. We will therefore simply replace s by s * = φs in Van Genuchten equation and see how it affects water retention properties.

This is equivalent as considering that the α vg coefficient of equation 1.29 depends on porosity:

α(φ) = α vg ref φ φ ref (1.31) S e =     1 1 + (α vg ref φ φ ref s) nvg     mvg (1.32)
Suction as a function of the effective degree of saturation is expressed as:

s = φ ref φα vg ref S -1 mvg e -1 1 nvg (1.33)
This is the simplest way to incorporate modified suction into the WRC formula. Moreover this approach has the advantage of not requiring new parameters. We will now compare WRCs obtained by this formula with porosity dependent WRCs found in the literature.

Boom clay: Water retention curves used by [START_REF] Delahaye | Soil heterogeneity and preferential paths for gas migration[END_REF] We can see that our proposition underestimates the influence of porosity on water retention properties, especially for porosities larger than the reference porosity. However, since the calculated WRCs are situated between the reference porosity WRC and the ones calculated by Delahaye and Alonso, we can consider that it is pertinent to incorporate suction effects into WRCs by replacing s by s * , especially since this solution doesn't require the determination of additional parameters.

Compacted Speswhite kaolin: Gallipoli et al. (2003b) proposed the following relationship:

S e = 1 1 + (α vg (v -1) β s) nvg mvg (1.34) in which v -1 = e = φ 1 -φ
with the following values for compacted speswhite kaolin :

m vg = 0.03586, n vg = 3.746, α vg = 0.02691 kPa -1 , β = 8.433
We will take as a reference the values for a porosity φ ref = 0.5, which gives

α vg ref = α vg φ ref φ ref -1 β = 0.093.
Figure 1.10 shows the WRCs obtained by equation 1.34 (Gallipoli et al., 2003b), which can be compared with the ones given by our proposition (equation 1.32) shown in figure 1.11.

With our proposal, we can see that the changes in water retention curves due to small variations of porosity is seriously underestimated. However, since the induced curve shift follows the logical direction and that it helps ensuring thermodynamical consistency, we will keep this formulation in the following work. The constitutive stress, σ * , varies with the product S l s. It can therefore be interesting to study how this component evolves with suction.

Van Genuchten parameters for different materials are given in table 1.2 [START_REF] Leong | Review of Soil-Water Characteristic Curve Equations[END_REF].

Material

S r α vg (kPa

(φ = φ ref )
The suction component of the constitutive stress as well as the WRC are represented in figures figures 1.12 to 1.15 for the different materials.

For coarse materials, such as mine tailings for instance, this component reaches a maximum for an intermediate saturation state during desaturation, and then decreases when the dry state is approached. It will therefore never reach high values and suction will thus have a low influence on the material mechanical behaviour.

Concerning fine-grained soils, such as clays, however, this component tends to infinity for low saturation states. [START_REF] Alonso | A microstructurally based effective stress for unsaturated soils[END_REF] saw this as a major drawback from this simple expression for the constitutive stress since, as suction increases, it predicts an unrealistic compression. However, this issue only concerns very high suctions and we will consider that the previously chosen approach remains valid as long as the degree of saturation does not take low values. Introduced by [START_REF] Kachanov | Time of the rupture process under creep conditions[END_REF], the effective stress in the sense of damage mechanics is based on the fact that the resisting section decreases when micro-cracks develop.

Since the applied force, F , at the sample scale is not affected by damage, the effective stress, σ, applied on the remaining section, S ef f , increases (see figure 1.16)

F = σS = σS ef f (1.

Introduction of a double effective stress

The two previous sections allowed us to introduce two quantities describing the stress applied on the solid matrix. On one hand, the constitutive stress, in unsaturated soils, takes into account the effect of water menisci in tension, acting like a compressive stress on the solid matrix. On the other hand, the effective stress, in the sense of damage mechanics, allows us to account for the decreasing material surface sustaining mechanical loads, resulting from the creation of micro-cracks. There is a need to define a new quantity, representing the stress applied on the solid matrix when the material is affected by both suction and damage simultaneously. This quantity will be called the double effective stress and is assumed to control the porous material mechanical behaviour.

Two simple combinations of the previous effective stresses can be imagined to incorporate both damage and suction into this double effective stress, σ * :

σ * 1 = σ -p g I + S l sI 1 -d = σ * 1 -d (1.39) σ * 2 = σ 1 -d -p g I + S l sI = σ -p g I + S l sI (1.40)
To choose between these two expressions, we consider that a damaged sample submitted to a change in suction should behave differently compared to the intact sample. This hypothesis has been considered by other authors, such as [START_REF] Carmeliet | Poromechanical modelling of shrinkage and damage processes in unsaturated porous media[END_REF], who consider that damaged materials experience more swelling when wetted that intact ones.

Considering that the total applied stress, the gas pressure, as well as damage are kept constant, the change in the double effective stress due to a suction increment would be:

σ * 1 = ( Ṡl s + S l ṡ)I 1 -d (1.41) σ * 2 = ( Ṡl s + S l ṡ)I (1.42)
Since elastic strains are directly related to the double effective stress, for the second expression, the strain change due to suction change would be the same for an intact and a damaged sample.

We will thus choose the first expression for the double effective stress:

σ * = σ -p g I + S l sI 1 -d (1.43)
We define the following quantities:

❼ Mean stress: p = 1 3 tr(σ)
❼ Deviatoric stress tensor:

σ d = σ -pI ❼ Deviatoric stress: q = 3 2 σ d : σ d
Then the double effective triaxial variables are:

p * = p -p g + S l s 1 -d (1.44) q * = q = q 1 -d (1.45)
It can be noted that with this definition of the double effective stress, suction effects are isotropic, and thus don't have any impact on the deviatoric stress.

The existence of a double effective stress, in which suction and damage effects on mechanical behaviour are included, is a key assumption in the following modelling developments. It will be shown in the next chapters how this double effective stress allows for damage and suction effects on elastic and dissipative behaviours to be reproduced. Of course its expression could be refined based on micromechanical studies. It will be shown, however, that this expression already allows us to capture well enough the main behaviour features observed during coupled hydromechanical loadings.

Usual elasto-plastic constitutive models for soils assume that soil behaves elastically for small strains. Although this assumption is idealistic, it is still reasonable, especially for overconsolidated clays, to consider an elastic domain in which the mechanical behaviour is reversible. Even for soils undergoing large plastic strains, the influence of the elastic part of the constitutive model should not be underestimated. Indeed, when modelling boundary value problems, the elastic component of the elastoplastic constitutive model determines the stress state of the soil at the onset of plasticity, which has an non negligible influence on the final plastic strains.

The small-strain stiffness of soils is known to be highly non-linear and has been experimentally studied extensively [START_REF] Viggiani | Stiffness of fine-grained soil at very small strains[END_REF][START_REF] Jovicic | Stiffness of coarse-grained soils at small strains[END_REF][START_REF] Rampello | Small-strain stiffness of reconstituted clay compressed along constant triaxial effective stress ratio paths[END_REF][START_REF] Atkinson | Non-linear soil stiffness in routine design[END_REF][START_REF] Callisto | Shear strength and small-strain stiffness of a natural clay under general stress conditions[END_REF]. Dynamics methods are usually used to deduce the elastic shear stiffness from shear waves velocities. The resonant column test consists in measuring the response of a cylindrical sample subjected to forced harmonic torsional vibrations. Another technique consists in transmitting and receiving shear waves using bender elements (electro-mechanical transducers). The small-strain shear modulus is usually assumed dependent on the stress state, the porosity (through the specific volume or the void ratio) and the loading history (through the overconsolidation ratio, the preconsolidation pressure or plastic strains).

In this chapter, we consider the initial elastic stiffness to depend on the double effective stress state. It will be shown how relating the double effective stress to elastic strains through constitutive equations results in a dependence of the apparent stiffness on damage and suction.

Preliminaries

Triaxial variables: In the following, potentials are written in terms of triaxial variables, so that σ :

dε = pdε v + qdε s = pdε v + σ d : dε d (2.1) p = 1 3 tr(σ) (2.2) q = 3 2 σ d : σ d (2.3) ε v = tr(ε) (2.4) ε s = 2 3 ε d : ε d (2.5) σ d = σ -pI (2.6) ε d = ε - 1 3 ε v I (2.7)
These relationships are also true for elastic strains (ε e ) and effective stresses (σ * , σ and σ * ).

Elasticity definitions: Different forms can be used to define elasticity which can be summarised as follows [START_REF] Houlsby | Principles of hyperplasticity: an approach to plasticity theory based on thermodynamic principles[END_REF]:

A material is said to be elastic if the stress can be expressed as a single-valued function of the strain:

σ = f (ε e ) (2.8)
The incremental stress-strain relationship is then written : σ = ∂f (ε e ) ∂ε e εe (2.9)

If the stress-strain relationship is originally expressed in its incremental form, the material is described as hypoelastic:

σ = f (ε e ) εe (2.10)
If the stresses can be derived from a strain energy potential, the material is said to be hyperelastic: 2.11) which gives the following stress-strain relationship :

σ = ∂f (ε e ) ∂ε e (
σ = ∂ 2 f (ε e ) ∂ε e 2 εe (2.
12)

It should be noted that hyperelastic materials automatically respect also the elasticity and hypoelasticity definitions. Elastic materials are also hypoelastic but not necessarily hyperelastic. When starting from the definition of the incremental law for hypoelastic materials, it is not always possible to find a potential from which the hypoelastic formulation could be derived.

Elastic potentials: Different energy functions can be defined depending on which variables are used as independent state parameters.

The elastic strain energy (Helmholtz free energy) used in equation 2.11 is written as a function of strains and will be noted from now ψ e 0 (ε e ). Deriving the incremental constitutive equations from this potential results in a fourth-order rigidity tensor, D e (ε e ), being a function of strains. (2.17)

Linear elasticity: For linear elasticity, the elastic potential takes the form:

ψ e 0 (ε e ) = K 2 (ε e v ) 2 + 3G 2 (ε e s ) 2 = pε e v 2 + qε e s 2 = p 2 2K + q 2 6G = G e 0 (σ) (2.18)
which after derivation gives,

p q = K 0 0 3G ε e v ε e s (2.19)
and its incremental form:

dp dq = K 0 0 3G dε e v dε e s (2.20)

Pressure dependence of moduli in Cam-Clay

Experimental evidence shows that, for geomaterials, bulk and shear moduli vary with confining pressure. This non-linearity is sometimes attributed to the closure of micro-cracks or to imperfect bonding between grains [START_REF] Zimmerman | Compressibility of sandstones[END_REF].

It is usually taken into account in Cam-Clay based models [START_REF] Roscoe | On the generalized stress-strain behaviour of 'wet' clay[END_REF]) by choosing the following expression for the bulk modulus:

K = 1 + e κ CC p (2.21)
To complete the stress-strain relationship, either the shear modulus of Poisson's ratio are assumed to take a constant value.

Assuming a constant Poisson's ratio gives a pressure dependant shear modulus which is in accordance with experimental data [START_REF] Baldi | Developments in modelling of thermohydro-geomechanical behaviour of Boom clay and clay-based buffer materials[END_REF].

G = 3(1 -2ν) 2(1 + ν) K = 3(1 -2ν) 2(1 + ν) 1 + e κ CC p (2.22)
However, this hypoelastic formulation does not derive from a thermodynamic potential. It has been demonstrated by [START_REF] Zytynski | On modelling the unloadingreloading behaviour of soils[END_REF] that this model is non-conservative since it can produce permanent shear strains during a closed stress cycle (see figure 2.1). One could argue that this would have an impact only under cyclic loading. However, [START_REF] Einav | Pressure-Dependent Elasticity and Energy Conservation in Elastoplastic Models for Soils[END_REF] compared tunnelling induced settlement predicted with either an hypoelastic model or a hyperelastic model and highlighted significant differences between both results. ε a (%) ε r (%) Fig. 2.1: Permanent strains resulting from a closed stress cycle in the elastic domain assuming a constant Poisson's ratio within a hypoelastic formulation On the contrary, choosing a constant shear modulus gives a conservative model. However, this does not agree with experimental observations and Poisson's ratio may take negative values for low confining stresses which makes this choice not pertinent.

Pressure dependent hyperelastic models

In order to ensure thermodynamical consistency and energy conservation, several researchers have been working on developing hyperelastic potentials with shear moduli dependent on pressure [START_REF] Houlsby | The use of a variable shear modulus in elastic-plastic models for clays[END_REF][START_REF] Borja | Coupling plasticity and energy-conserving elasticity models for clays[END_REF][START_REF] Sulem | Elasto-plastic modelling of Red Wildmoor sandstone[END_REF][START_REF] Einav | Pressure-Dependent Elasticity and Energy Conservation in Elastoplastic Models for Soils[END_REF][START_REF] Houlsby | Elastic moduli of soils dependent on pressure: a hyperelastic formulation[END_REF].

Deriving constitutive equations from potentials, considering pressure dependent bulk and shear moduli, necessarily leads to a rigidity matrix comporting an additional coupling modulus, J:

dp dq = K J J 3G dε e v dε e s (2.23)
The presence of the coupling modulus, J, results in stress induced anisotropy. This means that, even if the material properties are initially isotropic, a change in mean stress can create deviatoric strains, and that a change in deviatoric stress only can lead to volumetric strains. [START_REF] Houlsby | The use of a variable shear modulus in elastic-plastic models for clays[END_REF], [START_REF] Borja | Coupling plasticity and energy-conserving elasticity models for clays[END_REF] and [START_REF] Houlsby | Elastic moduli of soils dependent on pressure: a hyperelastic formulation[END_REF] drew constant volume (undrained) stress paths in the pq plane. It gives curved lines instead of the straight lines usually obtained when taking J * = 1. This form of the undrained stress path is in accordance with experimental results [START_REF] Shaw | Behaviour of dry granular materials under repeated load biaxial and triaxial stress conditions[END_REF][START_REF] Borja | Coupling plasticity and energy-conserving elasticity models for clays[END_REF], which advocates for the use of this form for the rigidity matrix.

Models are usually derived in a way to recover the following expression for the bulk modulus along the isotropic axis:

K p r = 1 κ p p r n = k p p r n (2.24)
This expression 2.24 is slightly different from equation 2.21, since, for n = 1, it gives straight swelling lines in the plane ln v -ln p * , instead of in the plane v -ln p * as commonly assumed. [START_REF] Butterfield | A natural compression law for soils (an advance on e -log p')[END_REF] showed that this gives a better adequation with experimental results, and [START_REF] Hashiguchi | On the linear relations of v-ln p and ln v-ln p for isotropic consolidation of soils[END_REF] advanced several advantages of this choice. This expression also makes the use of the natural (logarithmic) strains which allows for an easier extension to finite strains.

The two main hyperelastic formulations found in the literature are presented in table 2.1.

Houlsby (1985)

Houlsby et al. ( 2005)

ψ e 0 p r κ exp ε e v κ + 3 2 αp r exp ε e v κ ε e d : ε e d p r [k(1 -n)] 2-n 1-n k(2 -n) (ε * v ) 2 + 3gε e s 2 k(1 -n) 2-n 2-2n ε * v = ε e v + 1 k(1 -n) G e 0 p (2-n) e p (1-n) r k(1 -n)(2 -n) - p k(1 -n) p 2 e = p 2 + k(1 -n)q 2 3g K p κ 1 - n(1 -n)kq 2 3gp 2 e kp (1-n) r p n e G q √ ε d : ε d   1 1 -n(1-n)kq 2 3g(p 2 e -np 2 )   gp (1-n) r p n e J q κ nkpqp (1-n) r p (n-2) e
Table 2.1: Hyperelastic pressure dependent models from the literature

The first formulation of [START_REF] Houlsby | The use of a variable shear modulus in elastic-plastic models for clays[END_REF] (also used by [START_REF] Borja | Coupling plasticity and energy-conserving elasticity models for clays[END_REF]) exhibits a maximum attainable stress ratio (q/p) max = 3ακ/2 (the rigidity matrix becomes singular).

Moreover, no simple form of the corresponding Gibbs energy exists, and the rigidity matrix can therefore not be expressed in terms of stresses only. [START_REF] Houlsby | Elastic moduli of soils dependent on pressure: a hyperelastic formulation[END_REF] derived a new expression for which both Gibbs and Helmholtz energies are given. This formulation also has the advantage of letting the freedom to choose any value between 0 and 1 for the mean pressure exponent.

The resulting deviatoric and volumetric strains occurring during a closed stress cycle are represented in figure 2.2.

Since the second potential [START_REF] Houlsby | Elastic moduli of soils dependent on pressure: a hyperelastic formulation[END_REF] is expressed in terms of stresses and allows the pressure exponent, n, to be chosen in the whole range from 0 to 1, we will use this one in the following. For general stress states, the rigidity tensor is expressed as:

D e ijkl = p r p e p r n nk σ ij σ kl p 2 e + k(1 -n)δ ij δ kl + 2g δ ik δ jl - 1 3 δ kl δ ij (2.25)

Effect of damage on elastic behaviour

Damage modelling is based on the existence of an effective stress, as defined by [START_REF] Kachanov | Time of the rupture process under creep conditions[END_REF], such as

σ = σ 1 -d (2.26)
To derive the elastic behaviour of the damaged materials, two main assumptions have been used in previous works:

❼ The principle of strain equivalence, which states that the strain associated with a damaged state under the applied stress is equivalent to the strain associated with its undamaged state under the effective stress [START_REF] Lemaitre | Aspect phénoménologique de la rupture par endommagement[END_REF]. This principle is represented on figure 2.3.

❼ The principle of equivalent elastic energy which stipulates that the elastic energy of the damaged material is the same in form as that of an undamaged material except that the stress is replaced by the effective stress [START_REF] Cordebois | Endommagement anisotrope en élasticité et plasticité[END_REF].

Damaged material

Intact matrix

Classical elasticity and plasticity laws Hansen and Schreyer (1994) compared these two assumptions. In figure 2.4, damage is initiated at point 2 of the stress path 1-2-3. Point 7 denotes the effective state using the principle of strain equivalence. The stored energy associated with the final state, given by the area 1-3-6, is equated to the area in the effective space using the undamaged modulus, which is given by the area 1-5-4. This results in an effective state defined by an effective stress and an effective strain (point 5). The effective stress obtained by the principle of strain equivalence is much higher than the one for the principle of equivalent elastic energy. The choice between these two assumptions will therefore strongly influence the model behaviour.

In the numerical algorithm (see chapter 5.6), the stress-strain behaviour is solved within the double effective stress space and conversion back to total stresses is the last step of the algorithmic process. Apparent rigidities are therefore not explicitly needed.

However, to show the effect of damage and pressure dependency on apparent moduli, we will study them along the isotropic compression line (i.e. for q = 0) for both principles.

Along the isotropic axis, apparent moduli reduce to:

K app (σ, d = 0) = kp r p p r n (2.27) G app (σ, d = 0) = gp r p p r n
(2.28)

J app (σ, d = 0) = 0 (2.29)
A graphical definition of these moduli is given in figures 2.5 and 2.6. The ratios, R K and R G , between the apparent damaged and intact moduli are defined as:

R K = K app (p, d) K app (p, d = 0) (2.30) R G = G app (p, d) G app (p, d = 0) (2.31)

Principle of strain equivalence

The principle of strain equivalence can be expressed as

ε e (σ, d) = ε e ( σ) = ∂G e 0 ( σ) ∂ σ (2.32)
and the incremental form,

∂ε e ∂σ = ∂ε e ∂ σ ∂ σ ∂σ = 1 1 -d ∂ 2 G e 0 ( σ) ∂ σ2 (2.33)
which gives the relationship between the damaged and the intact compliance matrices,

C e (σ, d) = 1 1 -d C e 0 ( σ) (2.34)
and the rigidity matrix,

D e (σ, d) = [C e (σ, d)] -1 = (1 -d) [C e 0 ( σ)] -1 = (1 -d)D e 0 ( σ) (2.35)
hence the apparent moduli:

K app = (1 -d)K σ 1 -d = (1 -d)kp r p (1 -d)p r n = (1 -d) 1-n kp r p p r n (2.36) G app = (1 -d)G σ 1 -d = (1 -d)gp r p (1 -d)p r n = (1 -d) 1-n gp r p p r n (2.37)
Bulk and shear moduli are affected similarly by damage and by the pressure exponent, n.

A single ratio can therefore be defined to represent rigidity degradation.

R = R K = R G = (1 -d) (1-n) (2.38)
The evolution of R with damage for different values of n is given in figure 2.7.

Principle of equivalent elastic energy

The principle of equivalent elastic energy can be expressed as

G e (σ, d) = G e 0 ( σ) (2.39)
which gives after derivation,

ε e = ∂G e 0 ( σ) ∂σ = ∂G e 0 ( σ) ∂ σ ∂ σ ∂σ = 1 1 -d ∂G e 0 ( σ) ∂ σ (2.40)
and the incremental form,

∂ε e ∂σ = ∂ε e ∂ σ ∂ σ ∂σ = 1 (1 -d) 2 ∂ 2 G e 0 ( σ) ∂ σ2 (2.41)
which gives the relationship between the damaged and the intact compliance matrix,

C e (σ, d) = 1 (1 -d) 2 C e 0 ( σ) (2.42)
and the rigidity matrix:

D e (σ, d) = [C e (σ, d)] -1 = (1 -d) 2 [C e 0 ( σ)] -1 = (1 -d) 2 D e 0 ( σ) (2.43)
hence the apparent moduli

K app = (1 -d) 2 K σ 1 -d = (1 -d) 2 kp r p (1 -d)p r n = (1 -d) 2-n kp r p p r n (2.44) G app = (1 -d) 2 G σ 1 -d = (1 -d) 2 gp r p (1 -d)p r n = (1 -d) 2-n gp r p p r n (2.45)
Bulk and shear moduli are again affected similarly by damage and by the pressure exponent, n. The ratio between damaged and intact rigidities is therefore:

R = R K = R G = (1 -d) (2-n) (2.46)
The evolution of R with damage for different values of n is given in figure 2.8.

Comparison of both assumptions

Figures 2.7 and 2.8 both show that the influence of damage on apparent moduli decreases when n is higher. It is due to the fact that when effective stresses increase, the rigidity of the intact matrix increases too since it is pressure dependent. This phenomenon partly compensates the rigidity degradation due to the diminution of effective surface on which loads are applied.

We can also notice that this effect is much more important for the hypothesis of strain equivalence.

Figure 2.9 shows the evolution of radial stiffness during triaxial compression tests as a function of radial strain obtained by [START_REF] Chiarelli | Modeling of elastoplastic damage behavior of a claystone[END_REF]. The radial stiffness was obtained by calculating the slope of unloading-reloading loops.

Since axial strains include both elastic and plastic strains, it is difficult at this point to determine which of the two hypotheses would better fit experimental results. Indeed, Fig. 2.9: Evolution of radial stiffness during triaxial compression tests as a function of radial strain [START_REF] Chiarelli | Modeling of elastoplastic damage behavior of a claystone[END_REF] since stiffness degradation can be chosen as a way to calibrate the damage evolution law, both assumptions should be able to represent the same stiffness degradation with different damage evolution parameters.

However, when damage reaches 50 %, the stress-strain behaviour becomes unstable under increasing stress loading and reaches the softening part of the stress-strain curve. One should therefore check that the loss of stiffness obtained experimentally can be attained by the model for damage values under 0.5.

In the following developments we will use the hypothesis of strain equivalence since it provides a straightforward way to couple damage and plasticity phenomena.

Unsaturated elastic behaviour

This section aims to study how suction affects the elastic response to mechanical loads, and also how hydric loading, represented by a suction change, induces elastic strains.

The incremental stress-strain relationship expressed in terms of the constitutive stress is:

σ * = D e (σ * ) εe (2.47)
with

D e ijkl (σ * ) = p r p * e p r n nk σ * ij σ * kl p * e 2 + k(1 -n)δ ij δ kl + 2g δ ik δ jl - 1 3 δ kl δ ij (2.48)
To simplify the following expressions, we assume a constant null gas pressure (p g = 0).

Effect of suction on the mechanical elastic rigidity

The increment of total stress due to a mechanical load at constant suction is:

σ = σ * = D e (σ * ) εe (2.49)
Let us now compare the unsaturated rigidity matrix with the saturated one. The apparent rigidity matrix will be studied for isotropic stress states (i.e. when q = 0). Replacing stresses by constitutive stresses, the formula of table 2.1 become:

p * e = p * = p + Sls (2.50)

K app unsat = kp (1-n) r p * n (2.51) G app unsat = gp (1-n) r p * n (2.52)
Then the ratio between unsaturated and saturated moduli is:

R unsat = p * p n = p + S l s p n (2.53) Figures 2.10 to 2.
12 show the evolution of this ratio with mean pressure and with S l s. It can first be noticed that soil stiffen against mechanical loading when suction increases, except in the case where linear elasticity is assumed. The rigidity change due to suction is higher for high values of the pressure exponent, n. Moreover, the influence of suction on rigidity decreases when the mean pressure increases. Indeed, mean pressure effects then become predominant compared to suction effects. Fig. 2.12: p = 1e7 Pa

Elastic strains due to suction changes

The increment of elastic strain due to a suction increment under constant total stress is:

εe = D e (σ * ) -1 σ * = D e (σ * ) -1 ( Ṡl s + S l ṡ)I (2.54)
A significant advantage of the effective stress approach is the ability to capture suction induced strains without the need of extra parameters in addition to mechanical parameters.

This approach is validated against an oedometer swelling test performed by [START_REF] Volckaert | Thermalhydraulic-mechanical and geochemical behaviour of the clay barrier in radioactive waste repositories (model development and validation)[END_REF] Experimental values of the resulting swelling strains are represented in figure 2.14. Swelling strains are also computed thanks to equation 2.47, giving the relationship between the constitutive stress increment and the elastic strain increment. The numerical simulation is in good agreement with experimental data. Knowledge of the water retention properties in addition to mechanical rigidity parameters therefore allows us to reproduce adequately the elastic swelling behaviour observed during wetting.

General elasticity law

The two previous sections have highlighted how suction and damage, independently, affect the elastic stiffness.

Combining the definition of the double effective stress with the principle of strain equivalence gives the following elasticity law, accounting for both damage and suction effects:

σ * = D e ( σ * ) εe (2.55)
with

D e ijkl ( σ * ) = p r p * e p r n nk σ * ij σ * kl p * 2 e + k(1 -n)δ ij δ kl + 2g δ ik δ jl - 1 3 δ kl δ ij (2.56)
It has been shown that this formulation, combining the mechanical properties of the intact, dry, material, with the definition of a double effective stress representing the stress applied on the solid skeleton, allows us to adequately represent the principal effects of damage and suction on the apparent mechanical behaviour in the reversible domain.

The following chapter will deal with dissipative phenomena, i.e. the evolution of damage and plastic strains under mechanical loading. The principle of strain equivalence, despite its limitations in reproducing rigidity degradation, will be used in the following developments due to its ability to be easily coupled with plasticity.

CHAPTER 3

DISSIPATIVE PHENOMENA (DAMAGE AND PLASTICITY)

The purpose of this chapter is to describe the dissipative part of the constitutive model. Damage criterion and evolution law are introduced as well as the plasticity framework.

Damage and plasticity equations are expressed in terms of the double effective stress. It is shown how damage and plasticity criteria evolve with suction and damage in the double effective stress space as well as in the total stress space.

Le but de ce chapitre est de décrire la partie dissipative du modèle constitutif. The dissipative behaviour of geomaterials can be due to two different phenomena. The first one is plasticity, which causes irreversible strains, the second one is damage, which can be seen as the creation of microcracks and expresses itself by the degradation of elastic properties. These behaviours can appear independently or simultaneously depending on the material. (See figures 3.1 to 3.3)

ε σ E 0 E 0 ε p ε e ε σ E 0 E 0 ε p ε e Fig. 3.1: Plastic ε σ E 0 E 0 E(Ω) ε e ε σ E 0 E 0 E(Ω) ε e Fig. 3.2: Brittle ε σ E 0 E 0 E(Ω) ε p ε e ε σ E 0 E 0 E(Ω) ε p ε e

Fig. 3.3: Coupling of damage and plasticity

Early damage models were developed for purely brittle materials [START_REF] Krajcinovic | The continuous damage theory of brittle materials -Part 1: General Theory[END_REF][START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture du béton de structure[END_REF][START_REF] Ortiz | A constitutive theory for the inelastic behavior of concrete[END_REF]. Some of them considered permanent strains due to residual opening of micro-cracks [START_REF] Halm | An anisotropic model of damage and frictional sliding for brittle materials[END_REF][START_REF] Swoboda | An energy-based damage model of geomaterials-I. Formulation and numerical results[END_REF][START_REF] Arson | A mixed damage model for unsaturated porous media[END_REF]. However, they didn't account for the large plastic strains observed in geomaterials. This led to the development of models combining damage and elastoplasticity, first for concrete which were later extended to rock behaviour. These models are based either on phenomenological [START_REF] Simo | Strain-and stress-based continuum damage models-I. Formulation[END_REF][START_REF] Ju | On energy-based coupled elastoplastic damage theories: constitutive modeling and computational aspects[END_REF][START_REF] Yazdani | Combined plasticity and damage mechanics model for plain concrete[END_REF][START_REF] Hansen | A thermodynamically consistent framework for theories of elastoplasticity coupled with damage[END_REF][START_REF] Al-Rub | On the coupling of anisotropic damage and plasticity models for ductile materials[END_REF][START_REF] Grassl | Damage-plastic model for concrete failure[END_REF][START_REF] Einav | Coupled damage and plasticity models derived from energy and dissipation potentials[END_REF] or on micromechanical approaches [START_REF] Zhu | Micromechanical analysis of coupling between anisotropic damage and friction in quasi brittle materials: Role of the homogenization scheme[END_REF][START_REF] Guéry | A micromechanical model of elastoplastic and damage behavior of a cohesive geomaterial[END_REF][START_REF] Zhu | A micromechanics-based thermodynamic formulation of isotropic damage with unilateral and friction effects[END_REF].

The following developments are based on thermodynamics, whose basic principles are recalled in the next section.

Thermodynamics

Second law

The second law of thermodynamics states that the change of entropy, S, is such that

Ṡ Q T (3.1)
which gives, with local quantities:

η + ∇ • q T 0 (3.2) ∇ • q T = ∇ • q T - q • ∇T T 2 (3.3)
Incorporation of the expression of internal energy (equation 1.11) and strain work rate (equation 1.18) leads to:

T η - q • ∇T T -∇ • q = u -σ * : ε + s * Ṡl -p f (3.4) T η -u + σ * : ε -s * Ṡl - q • ∇T T + p f 0 (3.5)

Helmholtz free energy

It is more convenient to use temperature as a state variable instead of entropy. We will thus use the specific Helmholtz free energy defined through a Legendre transform as:

ψ = u -T η (3.6) ψ = u -T η -η Ṫ (3.7)
We assume that the material state is described by the values of the following state variables:

The elastic strain, ε e , the degree of saturation, S l , temperature, T , damage, d, and a hardening variable, χ. Helmholtz free energy can therefore be written as:

ψ = ψ(ε e , S l , T, d, χ) (3.8) ψ = ∂ψ ∂ε e εe + ∂ψ ∂S l Ṡl + ∂ψ ∂T Ṫ + ∂ψ ∂d ḋ + ∂ψ ∂χ χ (3.9)

Clausius-Duhem Inequality -dissipation

We assume the following strain partition:

ε = ε e + ε p (3.10)
Combining equations 3.5 and 3.7 leads to the Clausius-Duhem inequality:

σ * : ε -s * Ṡl -( ψ + η Ṫ ) - q • ∇T T + p f 0 (3.11)
which gives, introducing equation 3.9:

σ * - ∂ψ ∂ε e : εe -η + ∂ψ ∂T Ṫ -s * + ∂ψ ∂S l Ṡl + σ * : εp - ∂ψ ∂d ḋ - ∂ψ ∂χ χ - q T • ∇T + p f 0 (3.12)
We can obtain the following constitutive equations by constructing specific loading paths, with εp = 0, ḋ = 0, χ = 0, ∇T = 0 and p f = 0: By analogy with the previous relationships, the following thermodynamic conjugate forces are defined:

Y d = - ∂ψ ∂d ξ = - ∂ψ ∂χ (3.16)
The pairs of conjugated variables are summarised in table 3.1.

Strain-like variables

Work conjugate variables ε e Elastic strain tensor σ * Average skeleton stress tensor S l Degree of saturation s * Modified suction

T Temperature η Entropy d Damage Y d χ Hardening variables ξ Table 3.1: List of variables
The dissipation is then given by

D = σ * : εp + Y d ḋ + ξ : χ - q T • ∇T + p f 0 (3.17)
We can assume that there exist several dissipation potentials, which are decoupled so that each of them is positive:

D = D M + D T + D f (3.18)
with:

❼ The mechanical (intrinsic) dissipation:

D M = σ * : εp + Y d ḋ + ξ : χ 0 ❼ The thermal dissipation: D T = - q T • ∇T 0 ❼ The dissipation due to fluid flow: D f = p f 0
The following work will be focused on the mechanical dissipation D M .

Hypotheses on the form of Helmholtz free energy

In this work, isothermal processes are considered.

The following form is assumed for the Helmholtz free energy:

ψ = ψ(ε e , S l , χ, d) = ψ e (ε e , d) + ψ l (S l ) + ψ p (ε p , χ) (3.19)
This choice leads to certain couplings between the different phenomena involved.

Concerning the damage-elastic part of Helmholtz free energy, we chose to use the form proposed by [START_REF] Ju | On energy-based coupled elastoplastic damage theories: constitutive modeling and computational aspects[END_REF]:

ψ e (ε e , d) = ψ e 0 (ε e )(1 -d) (3.20)
This form is in accordance with the principle of strain equivalence combined with the existence of an effective stress:

σ * = ∂ψ ∂ε e = ∂ψ e 0 ∂ε e (1 -d) (3.21) σ * = σ * 1 -d = ∂ψ e 0 ∂ε e (3.22)

Damage criterion and evolution

The thermodynamic force conjugated to damage is:

Y d = - ∂ψ ∂d = - ∂ψ e ∂d = ψ e 0 (ε e ) (3.23)
This expression of Helmholtz free energy implies that the force conjugated to damage only depends on elastic strains. This is equivalent as saying that it depends on the double effective stress state since they are related through the constitutive equation 3.13. This is a strong assumption since it implies that the initiation and evolution of damage is driven by elastic strains. Some researches proposed that damage should be created by plastic strains. However, our approach decoupling both dissipative phenomena allows both of them to appear independently, and thus to simulate brittle material as well as plastic materials. Moreover, we will see later that plastic processes can still indirectly influence damage by modifying the stress state (see part 7.2.1).

General framework

To describe damage initiation and evolution, a damage criterion, f d = 0, is defined. The equation f d = 0 defines a surface, which separates the double effective stress space into a domain in which no dissipative phenomena occur (inside) and a domain of non-attainable stress states (outside). This surface is assumed to be convex.

Ideally, this function should depend on the thermodynamic force conjugated to damage:

f d (Y d (ε e ), d) = 0 (3.24)
Checking the sign of this function allows to determine in which state the material is in relation to damage.

If the stress state lies inside the damage surface, or in the case in which the stress state lies on the damage yield surface but the point goes towards the inside of the yield surface (unloading), then there is no evolution of damage.

If the stress state lies on the damage criterion and the loading is increasing, damage increases and hardening causes the damage surface to extend so that the stress state always remains on the surface.

These conditions can be summarised as follow:

if f d < 0 or ḟd < 0 then ḋ = 0 if f d = 0 and ḟd = 0 then ḋ 0 (3.25)
The stress state has to remain on the damage criterion during loading, which means that, at any time, f d = 0 and ḟd = 0. Deriving equation 3.24 gives the consistency condition:

0 = ḟd = ∂f d ∂Y d Ẏd + ∂f d ∂d ḋ (3.26)
from which the damage evolution law can be deduced:

ḋ = - Ẏd ∂f d ∂Y d ∂f d ∂d (3.27)

Specific functions

Damage criterion based on the thermodynamic force conjugated to damage

The expression of the thermodynamic force conjugated to damage given by the hyperelastic potential of [START_REF] Houlsby | Elastic moduli of soils dependent on pressure: a hyperelastic formulation[END_REF] can be derived from equation 3.23 and table 2.1:

Y d = ψ e 0 (ε e ) = p r [k(1 -n)] 2-n 1-n k(2 -n) (ε * v ) 2 + 3gε e s 2 k(1 -n) 2-n 2-2n (3.28)
This expression can also be expressed in terms of stresses. Deriving Gibbs free energy gives:

ε e v = ∂G e 0 ( σ * ) ∂ p * = 1 k(1 -n) p * p (1-n) r p * n 0 -1 (3.29) ε * v = ε e v + 1 k(1 -n) = 1 k(1 -n) p * p (1-n) r p * n 0 (3.30) ε e s = ∂G e 0 ( σ * ) ∂q = q 3gp (1-n) r p * n 0 (3.31)
Introducing equation 3.30 and 3.31 into equation 3.28 leads to:

Y d ( σ * ) = p (n-1) r k(2 -n) p * ( 2-n) 0 = p (n-1) r k(2 -n) p * 2 + k(1 -n)q * 2 3g (2-n) 2
(3.32)

A simple choice for the damage criterion would be:

f d = Y d -C 0 = 0 (3.33)
The shape of this criterion is given in figure 3.4, in the plane of the normalised variables:

p = p * [C 0 k(2 -n)p (1-n) r ] 1 2-n (3.34) q = q * k(1-n) 3g [C 0 k(2 -n)p (1-n) r ] 1 2-n (3.35) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 p ¯q ¯¯Fig. 3.4: Shape of damage criterion Y d -C 0 = 0
The choice of this criterion has two disadvantages. The first one is that it makes it possible to damage the material by isotropic compression, which is arguable. The second and main issue is that damage will initiate for a lower deviatoric stress when the confining pressure is higher, which is in disagreement with experimental evidence. Indeed, geomaterials are known to be more brittle at low confining pressure and more plastic at high confining pressures.

Damage criterion based on Drucker-Prager

For this reason, we choose to use a more empirical damage criterion which will enable us to represent the fact that geomaterials are less damageable at high confining pressure.

We propose to base this criterion on Drucker-Prager, which gives the following simple form for the damage loading function:

f d (Y * d , d) = Y * d -C 0 -C 1 d (3.36) with Y * d = 3J 2 -C 2 I 1 = q -C 2 p * (3.37)
which gives, for the expression of the yield surface,

q -C 2 p * -C 0 -C 1 d = 0 (3.38)
in which C 1 is a hardening parameter. The lower C 1 , the faster d will increase with deviatoric stress.

The shape of the damage criterion in the effective stress space is given in figure 3.5 for different values of damage. It can be seen that when damage increases, the material is hardening with respect to effective stresses.

q 1 -d -C 2 p + sS l 1 -d -C 0 -C 1 d = 0 (3.39) q -C 2 (p + sS l ) -(1 -d)(C 0 + C 1 d) = 0 (3.40)
The shape of the corresponding damage criterion in the total stress space is given in figures 3.6 and 3.7.

Figure 3.6 shows the evolution of the damage criterion with damage. It can been seen that, although the intact fraction of the material is hardening, an apparent softening behaviour appears after a certain value of damage is reached, when considering total stresses.

Figure 3.7 shows the evolution of the damage criterion with suction. Although suction doesn't have an effect on the damage criteria in the effective stress space, when considering total stresses, suction increases the size of the undamaged domain.

∂f d ∂d = -C 1 (3.41) ∂f d ∂Y * d = 1 (3.42) Ẏ * d = ∂Y * d ∂ σ * σ * = ∂Y * d ∂ p * I 3 + ∂Y * d ∂ q 3 σd 2q : σ * = - C 2 3 I + σd q : σ * (3.43)
from which the damage evolution law can be deduced:

ḋ = - Ẏ * d ∂f d ∂Y * d ∂f d ∂d = 1 C 1 - C 2 3 I + σd q : σ * (3.44) If we define A d = 1 C 1 - C 2 3 I + σd q
, the damage evolution law becomes:

ḋ = A d ( σ) : σ * (3.45)
A typical evolution of damage with deviatoric stress during a triaxial loading is given in figure 3.8.

Plasticity

According to [START_REF] Jommi | Remarks on the constitutive modelling of unsaturated soils[END_REF], a possible way to construct a model for unsaturated soils, starting from a saturated one, may be summarised in two steps :

❼ the substitution of the average skeleton stress for effective stress ❼ introduction in the basic saturated elastoplastic model of the modifications necessary to take into account the effects of the interfaces on the overall mechanical behaviour According to [START_REF] Ju | On energy-based coupled elastoplastic damage theories: constitutive modeling and computational aspects[END_REF], it appears reasonable to state that the plastic flows occur only in the undamaged material by means of effective quantities. Therefore, the characterisation of the plastic response should be formulated in the damaged effective stress space and the stress tensor is replaced by the damaged stress tensor, σ, into the equations of plasticity. This follows the principle of strain equivalence.

We combine these two approaches and choose as a basis a model for saturated soils, in which we will replace the saturated effective stress by our double effective stress. A dependence of the yield criterion on suction will be added.

We will therefore need to define: ❼ A yield criterion, which determines the boundary of the elastic domain:

f p = f p ( σ * , χ) (3.46)
such as if f p < 0 or ḟp < 0 then εp = 0 if f p = 0 and ḟp = 0 then εp 0 (3.47)

❼ A plastic potential:

g p = g p ( σ * , χ) (3.48)
❼ A plastic flow rule (if the plastic potential is taken equal to the yield criterion, the flow rule is said to be associated):

εp = Λp ∂g p ∂ σ * (3.49)
Fig. 3.9: Plastic potential and direction of plastic flow ❼ A hardening variable, χ = χ(χ 0 , s), depending on suction and on χ 0 , the hardening parameter for the saturated state

❼ A hardening law: χ0 = f ( εp ) (3.50)
The saturated soil model that we chose to develop to highlight how our modelling assumptions affect the plastic behaviour is the Cam-Clay model (CC). It was at first proposed by [START_REF] Roscoe | Yielding of clays in states wetter than critical[END_REF] for clays and later modified by [START_REF] Roscoe | On the generalized stress-strain behaviour of 'wet' clay[END_REF]. The different components of this model are described below. In this way, it is similar to the most widely used model for unsaturated soils, the Barcelona Basic Model (BBM) [START_REF] Alonso | A constitutive model for partially saturated soils[END_REF]. However some slight modifications are introduced, most of them being already used by other researchers. The main difference is the fact that the BBM model is expressed in terms of net stresses, whereas ours is expressed in term of the unsaturated effective stress.

Any other plasticity model could be implemented following the same approach in order to fit better the behaviour of different materials.

Definition of the preconsolidation pressure and hardening law

The typical volumetric response for an isotropic compression on a saturated sample is given in figure 3.10. The specific volume, v = 1 + e, is often used to represent soil volumetric behaviour. Along path AB, the soil behaves elastically and we recover the stress-strain relationship already seen in chapter 2.3, for the case n = 1.

εv = εe v = - v v = κ ṗ p (3.51)
When the mean pressure reaches the preconsolidation pressure p 0 , the materials yields and undergoes plastic strains. The preconsolidation pressure is chosen to be the hardening parameter. Along path BC, the volumetric stress-strain relationship is:

εv = εe v + εp v = - v v = λ ṗ p (3.52)
Integration of these equations leads to straight swelling (equation 3.51) and consolidation (equation 3.52) lines in the plane ln v -ln p. This represents a slight alteration of the Cam-Clay model for which these lines were taken as straight in the v -ln p plane. This modification is supported by experimental evidence [START_REF] Butterfield | A natural compression law for soils (an advance on e -log p')[END_REF], and has been used by many authors [START_REF] Wroth | A critical state model for predicting the behaviour of clays[END_REF][START_REF] Borja | Coupling plasticity and energy-conserving elasticity models for clays[END_REF][START_REF] Sheng | A constitutive model for unsaturated soils : thermomechanical and computational aspects[END_REF]. [START_REF] Hashiguchi | On the linear relations of v-ln p and ln v-ln p for isotropic consolidation of soils[END_REF] presented physical arguments in favour of a bilogarithmic law.

When yielding along path BC, the preconsolidation pressure increases and remains equal to the mean pressure. Equations 3.51 and 3.52 therefore become:

εe v = κ ṗ0 p 0 (3.53) εv = εe v + εp v = λ ṗ0 p 0 (3.54)
from which the hardening law can be deduced:

εp v = εv -εe v = (λ -κ) ṗ0 p 0 (3.55) ṗ0 = p 0 λ -κ εp v (3.56)
Extension to n = 1: To consistently derive the hardening law in the case in which the pressure exponent is not equal to 1, the previous equation needs to be modified. It is assumed that the dependence of the normal consolidation line on mean pressure follows the same law as for the swelling line.

Swelling line:

εv = εe v = - v v = κ p (1-n) r ṗ p n (3.57)
Normal consolidation line:

εv = εe v + εp v = - v v = λ p (1-n) r ṗ p n (3.58)
It should be noted that these lines are not straight lines anymore (see figure 3.11). The hardening law becomes:

εp v = εv -εe v = (λ -κ) p (1-n) r ṗ0 p n 0 (3.59) ṗ0 = p n 0 p (1-n) r λ -κ εp v (3.60)

Dependence of the preconsolidation pressure on suction

To be able to reproduce the extension of the elastic domain with suction, the preconsolidation pressure has to depend on suction and on the saturated preconsolidation pressure, p 0 :

p * c = p * c (p 0 , s) (3.61)
When drawn in the plane p -s, this curve is called the Loading-Collapse (LC) curve.

Many expressions have been proposed in the literature. Some of them are summarised in table 3.2.

Reference [START_REF] Alonso | A constitutive model for partially saturated soils[END_REF]

) (BBM) p net c p r = p 0 p r λ(0)-κ λ(s)-κ (3.62) λ(s) = λ(0)[(1 -r) exp(-βs) + r] (3.63) Jommi (2000) p * c p r = p 0 p r λ(0)-κ λ(s)-κ + S l s (3.64) λ(s) = λ 0 (p * -S l s) p * (3.65) Wheeler et al. (2003) p * c = p 0 (3.66) Sheng et al. (2004) p * c = p r p 0 p r λ(0)-κ λ(s)-κ + S l s (3.67) λ(s) = λ(0)[(1 -r) exp(-βs) + r] (3.68) Sun et al. (2008) p * c p r = p 0 p r λ(0)-κ λ(s)-κ (3.69) λ(s) = λ(0) + λ s s p g + s (3.70)
Table 3.2: Different expressions used for the LC curve in the literature

In the following we will use the proposition of [START_REF] Sheng | A constitutive model for unsaturated soils : thermomechanical and computational aspects[END_REF], because this is the closest one to that used in BBM. Since most of experimental works are based on the BBM, the values of the different parameters obtained from experimental studies are usually given for this expression of the LC curve.

The shape of the LC curve in both total and effective stress space is given in figure 3.12.

Yield surface

Cam-clay models have been developed in the framework of Critical State Soil Mechanics [START_REF] Roscoe | On the yielding of soils[END_REF][START_REF] Schofield | Critical state soil mechanics[END_REF]). The critical state concept states that soils and other granular materials, if continuously distorted until they flow as a frictional fluid, will come into a well-defined critical state. At the onset of the critical state, shear distortions occur without any further changes in mean effective stress, deviatoric stress or void ratio. The critical state is described in the (p * , q) plane by the line of equation:

q = M p * (3.71)
Original Cam-Clay model: The original Cam-Clay model [START_REF] Roscoe | Yielding of clays in states wetter than critical[END_REF] is based on the assumption that the deviatoric strain is entirely plastic, which is why no shear modulus is defined in the original Cam-Clay model.

It also assumes that, for any stress state lying on the yield surface, the dissipated energy is independent of the volumetric strain and is equal to:

Ẇ p = M p εp s (3.72)
The yield surface is derived by applying the normality condition (in the original versions of CC, the flow rule is assumed to be associated), which gives:

f p = q + M p * ln p * p 0 (3.73)
in which p 0 is the preconsolidation pressure. The corresponding yield surface is represented in figure 3.14. The main disadvantage of this expression of the yield surface is the fact that its intersection with the isotropic axis (q = 0) is angular. It leads to numerical issues since the derivative of f p is not unique at this point. Moreover, plastic deviatoric strains are created during isotropic loading, which is counter-intuitive.

Modified Cam-Clay model: [START_REF] Burland | Correspondence on The yielding and dilation of clay[END_REF] modified the aforementioned model by stating that the assumption of isotropy requires that under isotropic stress (q = 0), there is no distortion ( εp s = 0), which leads to the following expression for the dissipated work:

Ẇ p = p ( εp v ) 2 + (M εp s ) 2 (3.74)
Applying the normality condition to equation 3.74 gives the equation of the yield surface:

f p = q2 -M 2 p * (p * c (p 0 , s) -p * ) (3.75)
in which p * c is the preconsolidation pressure. The shape of the yield criterion in the effective stress space is given in figure 3.15. The yield surface in the effective stress space does depend on suction, but not on damage.

q 2 -M 2 (p + S l s) [(1 -d)p * c -(p + S l s)] = 0 (3.76)
and the equation of the critical state line becomes (equation 3.71)

q = M (p + S l s) (3.77)
The shape of the yield surface and of the critical state line in the total stress space is given in figures 3.16 and 3.17.

Figure 3.16 shows the evolution of the yield surface with suction. We can see that, with respect to total stresses, the elastic domain increases with suction. Suction also induces an apparent cohesion. Fig. 3.17: Effect of damage on yield surface -Total stress space Figure 3.17 shows the evolution of the yield surface with damage. We can see that damage has a softening effect on plastic behaviour. Although plastic and damage dissipative potentials were assumed to be decoupled, we can see that the introduction of the double effective stress into plasticity equations leads to a direct coupling between damage and plasticity.

Plastic flow rule

Cam-clay models have first been designed to consider an associated flow rule. However, it used to overestimate K 0 values and many models (for example BBM) use a non-associated flow rule.

The plastic potential is defined as:

g p = ζ q2 -M 2 p * (p * c -p * ) (3.78)
in which ζ is calculated to respect the condition of zero lateral strain for K 0 stress states in the usual Cam-Clay model.

ζ =   1 1 - κ λ   M (M -9)(M -3) 9(6 -M ) (3.79)
The shape of the plastic potential is shown in figure 3.18. We can see that the non-associated flow rule will give less deviatoric plastic strains that the associated flow rule.

Derivation of the incremental law

In the following, we derive the incremental stress-strain relationship for the case in which both damage and plastic strains occur at the same time.

Following the method of [START_REF] Sheng | Finite element formulation and algorithms for unsaturated soils. Part I: Theory[END_REF], pore fluid pressures are treated like strain components so that we have:

σ = f ( ε, ṗg , ṗl ) (3.82)
We first express the relationship between the increment of double effective stress and the increment of total strains and suction. Damage is not involved while staying in the effective stress space.

The hypothesis of strain partition, combined with the plastic flow rule (equation 3.80), gives:

ε = εe + εp = εe + Λp ∂g p ∂ σ * (3.83)
The elasticity law, defined in chapter 2. We define the following quantities:

A f = ∂f p ∂ σ * (3.89) A g = ∂g p ∂ σ * (3.90) A = ∂f p ∂p * c ∂p * c ∂p 0 ∂p 0 ∂ε p v ∂g p ∂ p * (3.91) B = ∂p 0 ∂ε p v ∂g p ∂ p * (3.92) C = ∂f p ∂p * c ∂p * c ∂s (3.93)
the detailed expressions of the different derivatives are given in appendix A.

We obtain the expressions for the plastic multiplier,

Λp = A f : D e : ε + C ṡ A f : D e : A g -A (3.94) the plastic flow, εp = A g ⊗ (D e : A f ) A f : D e : A g -A ε + CA g A f : D e : A g -A ṡ (3.95) the hardening law, ṗ0 = B A f : D e : ε + C ṡ A f : D e : A g -A (3.96)
and the incremental double effective stress -strain relationship.

σ * = D ep ( σ * ) ε + W ep ( σ * ) ṡ (3.97)
with

D ep = D e 1 - A g ⊗ (D e : A f ) A f : D e : A g -A (3.98) W ep = -D e CA g A f : D e : A g -A (3.99)
To be able to use our model with the Newton-Raphson method, we need to express the increment of total stress as a function of the increments of total strains and suction.

Deriving the relationships between effective stresses and total stresses gives

σ * = (1 -d) σ * (3.102) σ * = (1 -d) σ * -σ * ḋ = (1 -d) σ * -σ * (A d : σ * ) (3.103)
which lead to the expression of total stress rate as a function of total strain rate and pore fluid pressures rate:

σ = [(1 -d)I -σ * ⊗ A d ](D ep ε + W ep ṡ) + ṗg I -S l (s) + ∂S l ∂s s I ṡ (3.104) σ = D dep ( σ * , d) ε + ṗg I -F s ( σ * , d) ṡ (3.105) with D dep = [(1 -d)I -σ * ⊗ A d ]D ep (3.106) F s = S l (s) + ∂S l ∂s s I -[(1 -d)I -σ * ⊗ A d ]W ep (3.107)
The incremental equation is then:

σ = D dep ( σ * , d) ε + F s ( σ * , d) ṗl + (I -F s ( σ * , d)) ṗg (3.108)
with the different coefficient taking the following values for the four possible loading cases:

❼ Elastic loading:

D dep = (1 -d)D e (3.109) F s = S l (s) + ∂S l ∂s s I (3.110)
❼ Damage loading:

D dep = [(1 -d)I -σ * ⊗ A d ]D e (3.111) F s = S l (s) + ∂S l ∂s s I (3.112)
❼ Plastic loading:

D dep = (1 -d)D ep (3.113) F s = S l (s) + ∂S l ∂s s I -(1 -d)W ep (3.114)
❼ Damage-Plastic loading:

D dep = [(1 -d)I -σ * ⊗ A d ]D ep (3.115) F s = S l (s) + ∂S l ∂s s I -[(1 -d)I -σ * ⊗ A d ]W ep (3.116)
In this chapter, damage and plasticity evolution laws have been presented. The assumption of a double effective stress, associated with the principle of strain equivalence, allows for a direct damage-plasticity coupling. Indeed, although damage and plasticity criteria are expressed in terms of the double effective stress, and consequently do not depend explicitly on damage and suction, they evolve with damage and suction in the total stress space.

This approach, consisting in defining two dissipative criteria in the double effective stress space, appears to be flexible. Each component of the model can be easily modified without any need for heavy changes in the numerical implementation. All that is needed is a definition for the double effective stress, and the derivatives of damage and plasticity functions (developed in appendix A). The numerical implementation will therefore make use of explicit algorithms in order to keep this flexibility.

The following chapter deals with fluid transfer equations, which are required to model fully coupled hydro-mechanical problems.

CHAPTER 4

HYDRAULIC TRANSFERS

In order to allow for the representation of fully coupled hydro-mechanical problems, transfer laws have to be defined. The fluid transfer properties of porous media are highly dependent on the microstructure and the saturation state, and the resulting problem is therefore non-linear. This chapter presents a review of some expressions used in the literature, and the different models chosen to be implemented in θ-Stock.

Les lois de transfert doivent être définies afin de pouvoir étudier des problèmes hydromécaniques complètement couplés. Les propriétés de transfert de fluides dans les milieux non saturés sont fortement dépendants de la microstructure et de l'état de saturation du milieu poreux. Le problème est donc non linéaire. Ce chapitre présente certaines des expressions utilisées dans la littérature, ainsi que les différents modèles choisis pour être implémentés dans θ-Stock. Clays and rocks are known to have very low permeabilities (10 -12 to 10 -21 m/s). Fluid transfers will therefore have a great impact on the materials mechanical behaviour since overpressures are likely to develop under mechanical loading, and inhomogeneous stress states will appear during hydric loading.

Fluid mass conservation equations are developed for water and air constituents, which will enable the derivation of equations for FEM implementation. Two phases, liquid and gaseous, are considered. The liquid phase is assumed to contain liquid water and dissolved air. The gaseous phase is assumed to be composed of air as well as water vapour. Transfer laws are assumed to be function of the porosity as well as the saturation state. Significant works about coupled hydro-mechanical flow in unsaturated porous media include the work of [START_REF] Philip | Moisture movement in porous materials under temperature gradients[END_REF], [START_REF] Sophocleous | Analysis of water and heat flow in unsaturated-saturated porous media[END_REF], [START_REF] Geraminegad | A coupled thermoelastic model for saturatedunsaturated porous media[END_REF], [START_REF] Olivella | Nonisothermal multiphase flow of brine and gas through saline media[END_REF], [START_REF] Thomas | Analysis of coupled heat, moisture and air transfer in a deformable unsaturated soil[END_REF], [START_REF] Schrefler | A coupled model for water flow, airflow and heat flow in deformable porous media[END_REF], ...

The constitutive laws used in θ-Stock have been presented in [START_REF] Gatmiri | θ-STOCK, a powerful tool of thermohydromechanical behaviour and damage modelling of unsaturated porous media[END_REF].

Water phase

The water mass balance equation is written as

ρm + ∇ • Q m = 0 (4.1)
in which ρ m is the average water density of the gas-liquid mixture, and Q m is the total moisture flux. These quantities include liquid water ( w ) and vapour ( vap ) components.

ρ m = ρ w φS l + ρ vap φ(1 -S l ) (4.2)
in which φ is the porosity, S l the degree of saturation, ρ w the water density and ρ vap the vapour density.

Q m = Q w + Q vap = ρ w V w + ρ vap V vap = ρ w (V w + V * vap ) (4.3)
in which V α is the velocity of the α component (α = w, vap), and

V * vap = ρ vap ρ w V vap .

Liquid water flow

Water is assumed to flow under suction gradients and gravity, following the generalised Darcy's law:

V w = -K w ∇ - s γ w + z (4.4)
in which s is the suction, γ w the unit weight of water (γ w = 9.81 kN.m -3 ), and z the vertical coordinate.

Assuming a constant value for γ w , equation 4.4 can be rewritten as follows:

V w = D P l ∇s -K w ∇z (4.5)
with

D P l = 1 γ w K w
The water permeability, K w , is assumed to be isotropic, and to depend on porosity and on the degree of saturation:

K w = k r (S l )K int (φ) (4.6)
in which K int is the intrinsic permeability, depending on the material microstructure, and k r is the water relative permeability, depending on the saturation state.

Intrinsic permeability: Several expressions of the intrinsic permeability as a function of the porosity φ, or the void ratio, e, have been proposed in the literature. K w0 is the water permeability for a reference porosity. [START_REF] Lloret | Consolidation of unsaturated soils including swelling and collapse behaviour[END_REF] proposed to use the following form,

K int (e) = K w0 10 αe (4.7)
which is the expression used in most of the models incorporated in θ-Stock [START_REF] Gatmiri | Nouvelle formulation de la surface d'état en indice des vides pour un modèle non linéaire élastique des sols non saturés[END_REF][START_REF] Gatmiri | θ-STOCK, a powerful tool of thermohydromechanical behaviour and damage modelling of unsaturated porous media[END_REF].

The Kozeny-Carman's formula,

K int (φ) = K w0 φ 3 (1 -φ) 2 (1 -φ 0 ) 2 φ 3 0 (4.8)
in which φ 0 is the reference porosity for which K int = K w0 , has been used by [START_REF] Gens | Analysis of a full scalein situ test simulating repository conditions[END_REF] and [START_REF] Delahaye | Soil heterogeneity and preferential paths for gas migration[END_REF]. [START_REF] Delage | On the thermal consolidation of Boom clay[END_REF] found a linear relationship between the intrinsic permeability and the porosity, for Boom clay (for a porosity between 30 and 40 %).

Although we chose to consider microstructure effects by making K w dependent on porosity only, it is recognised that this is a simplifying assumption since damage may have a great impact on permeability by creating a new porosity network. The dependence of permeability on damage is a complex problem, which deserves to be handled carefully.

Recent research on this subject includes the work of [START_REF] Shao | Coupling between anisotropic damage and permeability variation in brittle rocks[END_REF], [START_REF] Arson | A mixed damage model for unsaturated porous media[END_REF], [START_REF] Maleki | Numerical simulation of damage-Permeability relationship in brittle geomaterials[END_REF], [START_REF] Arson | Influence of damage on pore size distribution and permeability of rocks[END_REF] and [START_REF] Pereira | Retention and permeability properties of damaged porous rocks[END_REF]. These works consider more complex damage models accounting for anisotropy.

Relative permeability: It is well known that water permeability, as well as air permeability, varies strongly with the degree of saturation.

The evolution of relative permeabilities to air and water obtained by [START_REF] Vachaud | Air and water flow during ponded infiltration in a vertical bounded column of soil[END_REF] on a fine sand are shown in figure 4.1.

A very common expression for water relative permeability, following an empirical approach, is:

k r (S l ) = S l -S r 1 -S r m (4.9)
in which m is determined to fit experiment data, and takes different values in several works Fig. 4.1: Relative permeabilities to air and water [START_REF] Vachaud | Air and water flow during ponded infiltration in a vertical bounded column of soil[END_REF] (m = 3 [START_REF] Irmay | On the hydraulic conductivity of unsaturated soils[END_REF], m = 3.5 (Averjanov, 1950), m = 4 [START_REF] Corey | Measurement of Water and Air Permeability in Unsaturated Soil[END_REF]). [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF], on the basis of data from 50 soils, found values of m ranging from 2.5 to 24.5.

Another approach, called statistical by [START_REF] Fredlund | Predicting the permeability function for unsaturated soils using the soil-water characteristic curve[END_REF], consists in deriving the hydraulic conductivity from the water retention curve. This approach is based on the fact that both the permeability function and the water retention curve are primarily determined by the pore-size distribution of the soil. [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF] proposed to express the hydraulic conductivity as a function of the nor-

malised water content Θ = θ -θ r θ sat -θ r : k r (Θ) = √ Θ     Θ 0 1 s(Θ) dΘ 1 0 1 s(Θ) dΘ     2 (4.10)
Assuming that this relationship is also valid when replacing the normalised water content by the effective degree of saturation (equation 1.26), and applying Mualem's approach to Van Genuchten's water retention curve expression (equation 1.33), the corresponding permeability becomes [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF]:

k r (S l ) = S l 1 -(1 -S 1/mvg l ) mvg 2 (4.11)
which can also be expressed in terms of suction:

k r (s) = (1 -(αs) nvg-1 [1 + (α vg s) nvg ] -mvg ) 2
(1 + (α vg s) nvg ) mvg/3 (4.12)

Van Genuchten (1980) obtained very good predictions of the hydraulic conductivity of several soils (sandstone, silts, clay and loam) after fitting the parameters on the water retention curves.

Water vapour flow

The water vapour flow is assumed to follow a Fick's diffusion law in a tortuous medium.

The vapour diffusion is linked to the gradient of vapour density and the vapour flow follows the formulation proposed by of [START_REF] Philip | Moisture movement in porous materials under temperature gradients[END_REF]:

V vap = 1 ρ vap q vap (4.13)
with q vap the vector of vapour flux density (kg m -2 s -1 ):

q vap = -D 0 ν vap α vap φ(1 -S l )∇ρ vap (4.14) in which, ❼ D 0 (m 2 s -1
) is the molecular diffusivity of water vapour in air. The formula used in θ-Stock is the empirical equation given by [START_REF] Geraminegad | A coupled thermoelastic model for saturatedunsaturated porous media[END_REF]. in which T is the temperature (➦C), and p g the gas pressure (Pa).

❼ ν vap is the "mass-flow factor" introduced to allow for the mass flow of vapour arising from the difference in boundary conditions governing the air and vapour components of the diffusion system.

ν vap = p g p g -p vap (4.16)
ν vap is close to 1 at normal soil temperature.

❼ α vap is a tortuosity factor allowing for the extra path length (α vap = 1.024).

❼ ρ vap (kg m -3 ) is the density of water vapour. According to [START_REF] Philip | Moisture movement in porous materials under temperature gradients[END_REF],

ρ vap = ρ 0 vap h (4.17)
ρ 0 vap being the density of saturated water vapour, and h the relative humidity, related to suction through Kelvin's law,

h(s) = exp - sg γ w R vap T (4.18)
in which g is the gravitational constant, and R vap is the specific vapour constant (R vap = 461.5 m 2 s -2 K -1 ).

In [START_REF] Geraminegad | A coupled thermoelastic model for saturatedunsaturated porous media[END_REF], the saturated vapour density depends only on temperature (kg/m 3 ):

ρ 0 vap = 10 -3 exp 19.819 - 4975.9 T + 273.15 (4.19)
Since temperature is considered homogeneous, the spatial derivation of equation 4.17 becomes:

∇ρ vap = ρ 0 vap ∇h = ρ 0 vap ∂h ∂s ∇s = -hρ 0 vap g γ w R vap T ∇s (4.20)
Combining equations 4.13, 4.14 and 4.20 gives the water velocity,

V vap = 1 ρ vap D 0 ν vap α vap φ(1 -S l )hρ 0 vap g γ w R vap T ∇s (4.21)
Since experimental data (from [START_REF] Philip | Moisture movement in porous materials under temperature gradients[END_REF]) do not appear to show the choking of the vapour at high moisture content, [START_REF] Ewen | Heating unsaturated medium sand[END_REF] modified the previous expression, assuming that the vapour flow area is equal to the porosity, and thus replacing φ(1 -S l ) by φ in equation 4.21:

V vap = 1 ρ vap D 0 ν vap α vap φhρ 0 vap g γ w R vap T ∇s (4.22)
The normalised vapour velocity therefore becomes: Let now write φ, Ṡl , ρw and ρvap as functions of strain, suction and fluid pressure increments.

V * vap = ρ vap ρ w V vap = D Pvap ∇s (4.23) with D Pvap = 1.024 ρ w D 0 ν vap φhρ 0 vap g R vap T 1 γ w (4.24)

Water mass conservation equation

Due to the assumption of solid grains incompressibility and with the soil mechanics convention, the porosity increment can be expressed as a function of the strain increment: φ =εv = -I : ε (4.28)

Deriving Van Genuchten's equation (1.29) gives the increment of degree of saturation as a function of the suction increment:

Ṡl = ∂S l ∂s ṡ = -(1 -S r ) m vg n vg α nvg vg s nvg-1 (1 + (α vg s) nvg ) mvg+1 ṡ = g 2 ṡ (4.29)
The water density is assumed to be pressure dependent (temperature dependency is not considered since we only consider isothermal conditions). ρw = ∂ρ w ∂p l ṗl = β p ṗl (4.30) in which β p is the water compressibility with respect to pressure (β p = 5.10 -10 m -2 s 2 ).

The vapour density rate can be derived from equation 4.17: -

ρvap = ρ 0 vap ∂h ∂s ṡ = -hρ 0 vap g γ w R vap T ṡ = A ṡ (4.
[ρ w S l + ρ vap (1 -S l )]I : ε + [-φ(ρ w -ρ vap )g 2 + φS l β p -φ(1 -S l )A] ṗl + [φ(ρ w -ρ vap )g 2 + φ(1 -S l )A] ṗg + ∇ • ρ w (D P ∇s -K w ∇z) = 0 (4.33)

Air phase

The air mass balance equation is written as

ρa + ∇ • Q a = 0 (4.34)
in which ρ a is the average air density of the gas-liquid mixture, and Q a is the total air flux.

ρ a = ρ a φ(1 -S l + HS l ) (4.35)
in which ρ a is the intrinsic air density, and H = 0.02 is Henry's constant and is used to account for dissolved air into water.

In the porous space occupied both by air and water vapour, the vapour flux decreases the air flux. The total air flux therefore becomes:

Q a = ρ a V a + ρ a HV w -ρ w V * vap (4.36)
in which V a is the gaseous air velocity. V w and V * vap are respectively the liquid water velocity and the normalised vapour velocity, already defined by equations 4.5 and 4.23.

Gaseous air transfers are supposed to follow Darcy's law

V a = -K a ∇ p g γ a + z (4.37)
Similarly to water permeability, the air permeability depends on the porosity and on the porosity (or void ratio). [START_REF] Lloret | Consolidation of unsaturated soils including swelling and collapse behaviour[END_REF] proposed the following formula,

K a = c a γ a µ a [e(1 -S l )] αa (4.38)
in which γ a is the specific weight of air, µ a its viscosity, and c a and α a are constants. This expression has been used by [START_REF] Thomas | Analysis of coupled heat, moisture and air transfer in a deformable unsaturated soil[END_REF] and [START_REF] Gatmiri | Udam: A powerful finite element software for the analysis of unsaturated porous media[END_REF] and is implemented in θ-Stock.

The air flux can be expressed as a function of fluid pressure gradients:

Q a = -[ρ a (K a +HK w )∇z]-[ρ a HD P l -ρ w D Pvap ]∇p l + - ρ a K a γ a + ρ a HD P l -ρ w D Pvap ∇p g (4.39)
Deriving equation 4.35, the variation of the average water density can be expressed as

ρa = ρa φ(1 -S l + HS l ) + ρ a φ(1 -S l + HS l ) -ρ a φ(1 -H) Ṡl (4.40)
The increments of porosity and degree of saturation as a function of strain, suction and fluid pressure increments have already been defined by equations 4.28 and 4.29.

Assuming that air is an ideal gas, ρa = 1 R a T ṗg = α p ṗg (4.41)

with R a = 287 J.K -1 .kg -1 , the specific gas constant of dry air.

Incorporating equations 4.28, 4.29 and 4.41 into equation 4.35, the increment of air density becomes:

ρa = -(1-S l +HS l )ρ a I : ε-(H -1)φρ a g 2 ṗl +[(H -1)φρ a g 2 +φ(1-S l +HS l )α p ] ṗg (4.42)
Incorporation of equations 4.39 and 4.42 into equation 4.34 gives the mass conservation equation:

(1 -S l + HS l )ρ a I : ε

+ (H -1)φρ a g 2 ṗl -[(H -1)φρ a g 2 + φ(1 -S l + HS l )α p ] ṗg + ∇ • [ρ a (K a + HK w )∇z] + ∇ • [(ρ a HD P l -ρ w D Pvap )∇p l ] -∇ • [(- ρ a K a γ a + ρ a HD P l -ρ w D Pvap )∇p g ] = 0 (4.43)
The mass conservation equations developed in this chapter together with the mechanical incremental equation derived in chapter 3.4, form a system of partial differential to be solved to determine the solution of boundary value problems. These equations are highly non-linear since permeabilities depend on the saturation state and the porosity of the material. The following chapter will be devoted to the presentation of the numerical methods used to solve this system of equations.

The resolution of bi-dimensional coupled hydro-mechanical problems requires several discretisations. The first spatial discretisation is made possible by means of the use of the Finite Element method. A second discretisation in time is necessary due to the transient fluid flow. Moreover, due to the high non linearity of mechanical and transfer phenomena, loads have to be split into a certain number of increments, and iterative algorithms have to be used to reach a convergent solution.

The solution of the initial boundary value problem will give the fields of displacements, u, water pressure, p l , and air pressure, p g .

Strong form of the coupled problem

Synthesis of the equations

The problem is represented by a system of three equations, corresponding to the solid skeleton equilibrium and to the mass conservation of water and air.

Solid skeleton:

The equilibrium equation can be written as:

∇ • σ + ḃ = 0 (5.1)
in which b is the body force vector.

Incorporating the incremental mechanical law (equation 3.108) into the previous equation gives:

∇ • [D dep ( σ * , d) : ε] + ∇ • [F s ( σ * , d) ṗl ] + ∇ • [(I -F s ( σ * , d)) ṗg ] + ḃ = 0 (5.2)
Water:

The mass conservation equation of pore water is,

ρm + ∇ • Q m = 0 (5.3)
which becomes, expressing the average water density rate in terms of strain and pore pressures rates, and the total moisture flux in terms of vertical coordinates and pore pressure gradients,

-[ρ w S l + ρ vap (1 -S l )]Id : ε + [-φ(ρ w -ρ vap )g 2 + φS l β p -φ(1 -S l )A] ṗl + [φ(ρ w -ρ vap )g 2 + φ(1 -S l )A] ṗg -∇ • (ρ w K w ∇z) -∇ • (ρ w D P ∇p l ) + ∇ • (ρ w D P ∇p g ) = 0
(5.4)

Air:

The mass conservation equation of air is,

ρa + ∇ • Q a = 0 (5.5)
which becomes, expressing the average air density rate in terms of strain and pore pressures rates, and the total air flux in terms of vertical coordinates and pore pressure gradients,

-(1 -S l + HS l )ρ a Id : ε -(H -1)φρ a g 2 ṗl + [(H -1)φρ a g 2 + φ(1 -S l + HS l )α p ] ṗg -∇ • [ρ a (K a + HK w )∇z] -∇ • [(ρ a HD P l -ρ w D Pvap )∇p l ] + ∇ • [(- ρ a K a γ a + ρ a HD P l -ρ w D Pvap )∇p g ] = 0
(5.6)

Boundary conditions

We note Ω the domain and Γ its boundary. Boundary conditions can be either of Dirichlet's type if values of variables are imposed, or of Neumann's type if it is their flux which is given.

Solid skeleton:

On the solid skeleton it is possible to impose either the value of displacements, u,

∀t, ∀x ∈ Γ u , u(x, t) = ũ(x, t) (5.7)
or the value of the normal stress, t.

∀t, ∀x ∈ Γ σ , σ(x, t) • n = t(x, t) (5.8) with Γ u ∪ Γ σ = Γ and Γ u ∩ Γ σ = ∅
In θ-Stock, imposed displacement can only take a null value.

Water and air:

Equations relative to air and water are similar and will be treated together. The α subscript can represent water, w , or air, a .

Given boundary conditions can be either the fluid pore pressure, ∀t > 0, ∀x ∈ Γ pα , p α (x, t) = pα (x, t) (5.9)

or the fluid flux normal to the surface,

∀t > 0, ∀x ∈ Γ Qα , Q α (x, t) • n = qα (x, t) (5.10) with Γ pα ∪ Γ Qα = Γ and Γ pα ∩ Γ Qα = ∅

Initial conditions

The initial state is described by the full field of displacements and pore pressures at t = 0.

∀x ∈ Ω, u(x, t = 0) = u 0 (x) (5.11) ∀x ∈ Ω, p α (x, t = 0) = p α 0 (x) (5.12)

Weak form

The following developments will require some mathematical tools. We consider that second order tensors are represented by capitalised bold letters, A, vectors by italic bold letters, a, and scalars by italic letters, a.

Divergence operator properties:

∇ • (A • f) = A : ∇f + f • (∇ • A) (5.13) ∇ • (af) = f • ∇a + a(∇ • f) (5.14)
Divergence theorem:

Ω (∇ • f)dΩ = Γ (f • n) dΓ (5.15)
in which n is the outward unit normal to the boundary surface Γ.

Displacements

The variational formulation is obtained by multiplying the equilibrium equation by u * , a kinematically admissible displacement field with zero boundary condition (u * |Γu = 0), and integrating over the domain Ω.

∀u * , Ω (∇ • σ) • u * dΩ + Ω ḃ • u * dΩ = 0 (5.16)
Using the symmetry of the stress tensor, the properties of the divergence operator, and the divergence theorem, the variational formulation becomes:

∀u * , Γ ( σ • n) • u * dΓ - Ω σ : ∇u * dΩ + Ω ḃ • u * dΩ = 0 (5.17)
Under the assumption of small strains, strains are defined by

ε(u * ) = 1 2 (∇u * + (∇u * ) T ) (5.18)
The symmetry of the stress tensor gives,

σ : ε(u * ) = σ : 1 2 (∇u * + (∇u * ) T ) = σ : ∇u * (5.19)
Boundary conditions are:

u * = 0 on Γ u (5.20) σ • n = ṫ on Γ σ (5.21)
Introducing boundary conditions and strains into the variational formulation gives:

∀u * , Ω σ : ε(u * )dΩ = Γσ ṫ • u * dΓ + Ω ḃ • u * dΩ (5.22)
We can now introduce the incremental law (equation 1.22) into equation 5.22. The weak form of the problem can be stated as follows:

Find u, p l and p g , such as,

∀u * , Ω (ε(u * )) T : D dep ( σ * , d) : ε(u) dΩ + Ω (ε(u * )) T : F s ( σ * , d) ṗl dΩ + Ω (ε(u * )) T : [(I -F s ( σ * , d)) ṗg dΩ = Γσ (u * ) T • ṫ dΓ + Ω (u * ) T • ḃ dΩ (5.23)

Fluid pressures

The mass conservation of fluid pressure is:

ρα + ∇ • Q α = 0 (5.24)
in which the α subscript is either related to water or air.

The variational formulation is obtained by multiplying the mass conservation equation by p * α , a kinematically admissible pressure field with zero boundary condition (p α * |Γp α = 0), and integrating over the domain Ω.

∀p * α , Ω ρα p * α dΩ + Ω (∇ • Q α )p * α dΩ = 0 (5.25)
Using the properties of the divergence operator and the divergence theorem, the variational formulation becomes:

∀p * α , Ω ρα p * α dΩ + Γ p * α (Q α • n)dΓ - Ω Q α • ∇p * α dΩ = 0 (5.26)
Boundary conditions are: p * α = 0 on Γ pα (5.27)

Q α • n = qα on Γ Qα (5.28)
Introducing boundary conditions into the variational formulation gives:

∀p * α , Ω ρα p * α dΩ + Γ Qα qα p * α dΓ - Ω Q α • ∇p * α dΩ = 0 (5.29)
Introducing the specific expression of the average densities and fluxes of air and water, the weak form of the water mass conservation is, Find u, p l and p g , such as,

∀p * l ∈ V p l (Ω), Ω -[ρ w S l + ρ vap (1 -S l )]Id : ε(u)p * l dΩ + Ω [-φ(ρ w -ρ vap )g 2 + φS l β p -φ(1 -S l )A] ṗl p * l dΩ + Ω [φ(ρ w -ρ vap )g 2 + φ(1 -S l )A] ṗg p * l dΩ = - Ω ρ w K w (∇p * l ) T • ∇zdΩ - Ω ρ w D P (∇p * l ) T • ∇p l dΩ + Ω ρ w D P (∇p * l ) T • ∇p g dΩ - Γ Qw p * l qw dΓ (5.30)
and the weak form of the air mass conservation is, Find u, p l and p g , such as,

∀p * g ∈ V pg (Ω), - Ω p * g (1 -S l + HS l )ρ a Id : ε(u)dΩ - Ω p * g (H -1)φρ a g 2 ṗl dΩ + Ω p * g [(H -1)φρ a g 2 + φ(1 -S l + HS l )α p ] ṗg dΩ = - Ω ρ a (K a + HK w )(∇p * g ) T • ∇zdΩ - Ω (ρ a HD P l -ρ w D Pvap )(∇p * g ) T • ∇p l dΩ + Ω (- K a g + ρ a HD P l -ρ w D Pvap )(∇p * g ) T ∇p g dΩ - Γ Qa p * g qa dΓ
(5.31)

Application of the Finite Element Method

Spatial discretisation

The initial boundary value problem will be treated in two dimensions, which allow us to study to types of problems, using a plane strain or an axisymmetric analysis. Four degrees of freedom (dofs) are considered.

For plane strain conditions, the dofs are the horizontal displacement, u x (x, t), the vertical displacement, u y (x, t), the water pore pressure, p l (x, t), and the air pore pressure, p g (x, t).

x = x y is the coordinates vector.

For axisymmetric conditions, the dofs are the radial displacement, u r (x, t), the axial displacement, u z (x, t), the water pore pressure, p l (x, t), and the air pore pressure, p g (x, t).

x = r z is the coordinates vector.

The following developments are made for the plane strain configuration, and are valid for axisymmetric analysis (by replacing x by r , y by z and z by θ ).

Space is discretised by means of four-nodes quadrilateral isoparametric 2D elements, which means that the dofs are interpolated by the same functions as the geometry. Coordinates are interpolated by bilinear shape functions:

x = N • x = N 1 x 1 + N 2 x 2 + N 3 x 3 + N 3 x 3 (5.32) y = N • y = N 1 y 1 + N 2 y 2 + N 3 y 3 + N 3 y 3 (5.33)
with the shape functions:

N =          N 1 N 2 N 3 N 4          =          (1 -ξ)(1 -η)/4 (1 + ξ)(1 -η)/4 (1 + ξ)(1 + η)/4 (1 -ξ)(1 + η)/4          (5.34)
Displacements and pore pressures are expressed in terms of their values at a finite number of points in space and their values between these points are approximated by using the same bilinear shape functions N:

               u x (x, t) = N • u x (t) u y (x, t) = N • u y (t) p l (x, t) = N • p l (t) p g (x, t) = N • p a (t) (5.35)
For each dof, a vector of its values at nodes is defined:

dof(t) =          dof 1 (t) dof 2 (t) dof 3 (t) dof 4 (t)          (5.36)
Due to the specific configurations studied, only 4 components of the stress and strain tensors are needed. Stress and strain tensors can then be represented by vectors. We define the stress and strain vectors, as well as the identity vector i:

σ =          σ xx σ yy σ xy σ zz          (5.37) ε =          ε xx ε yy 2ε xy ε zz          (5.38) i =          1 1 0 1          (5.39)
In plane strain conditions, ε zz = 0

The fourth order elasticity tensor is represented by a second order tensor:

D e =      D e
              εxx εyy 2 εxy εzz          (5.41)
Since spatial derivatives depend on the coordinates system, the strain operator takes different values for plane strain and axisymmetric analyses:

Plane strain : 5.46) Axisymmetry : 5.47) in which r m is the medium radius of the element.

B =             ∂N 1 ∂x 0 ∂N 2 ∂x 0 ∂N 3 ∂x 0 ∂N 4 ∂x 0 0 ∂N 1 ∂y 0 ∂N 2 ∂y 0 ∂N 3 ∂y 0 ∂N 4 ∂y ∂N 1 ∂y ∂N 1 ∂x ∂N 2 ∂y ∂N 2 ∂x ∂N 3 ∂y ∂N 3 ∂x ∂N 4 ∂y ∂N 4 ∂x 0 0 0 0 0 0 0 0             (
B =             ∂N 1 ∂r 0 ∂N 2 ∂r 0 ∂N 3 ∂r 0 ∂N 4 ∂r 0 0 ∂N 1 ∂z 0 ∂N 2 ∂z 0 ∂N 3 ∂z 0 ∂N 4 ∂z ∂N 1 ∂z ∂N 1 ∂r ∂N 2 ∂z ∂N 2 ∂r ∂N 3 ∂z ∂N 3 ∂r ∂N 4 ∂z ∂N 4 ∂r N 1 r m 0 N 2 r m 0 N 3 r m 0 N r m 0             (
Spatial derivatives of pore pressure will be given by:

   ∇(p l (x, t)) = ∇N p l (t)
∇(p g (x, t)) = ∇N p g (t) (5.48) Plane strain :

∇N =     ∂N 1 ∂x ∂N 2 ∂x ∂N 3 ∂x ∂N 4 ∂x ∂N 1 ∂y ∂N 2 ∂y ∂N 3 ∂y ∂N 4 ∂y     (5.49) Axisymmetric analysis : ∇N =     ∂N 1 ∂r ∂N 2 ∂r ∂N 3 ∂r ∂N 4 ∂r ∂N 1 ∂z ∂N 2 ∂z ∂N 3 ∂z ∂N 4 ∂z    
(5.50)

Numerical integration

Numerical integration is performed by a Gauss-Legendre quadrature. The integrands are evaluated at specific points of the element, then weighted, and summed.

1 -1 1 -1 f (ξ, η)dξdη ≈ n i=1 n j=1 w i w j f (ξ i , η j ) (5.51)
In Θ-Stock, four integration points are used:

n = 2 (5.52) ξ i , η i = ± 1 √ 3 (5.53) w 1 = w 2 = 1 (5.54)
The determinant of the Jacobian matrix J is used to transform the volume element from the cartesian coordinates to the natural coordinates: (5.56)

The numerical integration becomes:

In Plane strains:

Ω f dΩ = f (x, y)dxdy = 1 -1 1 -1 f (ξ, η)det(J)dξdη = f -, 1 √ 3 det(J)
(5.57)

In axisymmetric analysis:

Ω f dΩ = f (r, z)rdrdz = 1 -1 1 -1 f (ξ, η)det(J)rdξdη = f - 1 √ 3 , - 1 √ 3 + f 1 √ 3 , - 1 √ 3 + f - 1 √ 3 , 1 √ 3 + f 1 √ 3 , 1 √ 3 rdet(J)
(5.58)

System to solve

After spatial discretisation and numerical integration, the equations to solve can be written in the following matricial form: (5.59) in which the terms corresponding to the mechanical equilibrium are, (5.63) the terms corresponding to water flow are,

   R uu R uw R ua C wu C ww C wa C au C aw C aa         u(t) ṗl (t) ṗg (t)      +    0 0 0 0 K ww K wa 0 K aw K aa         u(t) p l (t) p g (t)      =      Ḟσ (t) F w (t) F a (t)     
R uu = Ω B T D dep ( σ * , d) B dΩ (5.60) R uw = B T F s ( σ * , d) ⊗ N dΩ (5.61) R ua = B T I d -F s ( σ * , d) ⊗ N dΩ (5.62) Ḟσ (t) = Γσ N u T ṫ dΓ + Ω N u T ḃ dΩ
C wu = - Ω (ρ w S l + ρ vap (1 -S l )) N T ⊗ i B dΩ (5.64) C ww = Ω [-φ(ρ w -ρ vap )g 2 + φS l β p -φ(1 -S l )A] N T ⊗ N dΩ (5.65) C wa = Ω [φ(ρ w -ρ vap )g 2 + φ(1 -S l )A] N T ⊗ N dΩ (5.66) K ww = Ω ρ w D P ∇N T ∇N dΩ (5.67) K wa = - Ω ρ w D P ∇N T ∇N dΩ (5.68) F w (t) = - Ω ρ w K w ∇N T ∇zdΩ - Γ Qw N T
qw dΓ (5.69) and the terms corresponding to the air flow are:

C au = - Ω (1 -S l + HS l )ρ a N T ⊗ i B dΩ (5.70) C aw = - Ω (H -1)φρ a g 2 N T ⊗ N dΩ (5.71) C aa = Ω [(H -1)φρ a g 2 + φ(1 -S l + HS l )α p ] N T ⊗ N dΩ (5.72) K aw = Ω (ρ a HD P l -ρ w D Pvap ) ∇N T ∇N dΩ (5.73) K aa = - Ω (- K a g + ρ a HD P l -ρ w D Pvap ) ∇N T ∇N dΩ (5.74) F a (t) = - Ω ρ a (K a + HK w ) ∇N T ∇zdΩ - Γ Qa N T qa dΓ (5.75)

Discretisation in time

The solution of the coupled system (equation 5.59) requires integration to be performed over a discrete number of time increments. Traditionally, most time-stepping schemes are based on the θ-method. In this procedure, all time-dependent quantities are represented by a weighted sum of the values at the start and end of the time step, respectively, with the rate quantities being approximated by average rates over the time step(see figure 5.3).

The value of the quantity, f , at a certain intermediate time,

t θ = t n + θ(t n+1 -t n ), is
determined by considering a linear interpolation between the two extremities of the time step:

f (t θ ) = (1 -θ)f n + θf n+1 = f n + θ(f n+1 -f n ) = f n + θ∆f (5.76)
It is then assumed that f (t) is constant on [t n ,t n+1 ] and is equal to f (t θ ). Hence the integrals:

t n+1 tn f (t)dt = t n+1 tn f (t θ )dt = (f n + θ∆f )∆t (5.77) t n+1
tn ḟ (t)dt = ∆f (5.78) After integration between t and t n+1 equation 5.59 becomes : (5.79) which can be rewritten in the following form:

   R uu R uw R ua C wu C ww C wa C au C aw C aa         ∆u ∆p w ∆p g      +    0 0 0 0 K ww K wa 0 K aw K aa         (u n + θ∆u)∆t (p wn + θ∆p w )∆t (p g n + θ∆p g )∆t      =      ∆F σ (F wn + θ∆F w )∆t (F an + θ∆F a )∆t     
   R uu R uw R ua C wu C ww + θ∆tK ww C wa + θ∆tK wa C au C aw + θ∆tK aw C aa + θ∆tK aa         ∆u ∆p w ∆p g      =      ∆F σ θ∆t∆F w θ∆t∆F a      + ∆t      0 F wn -K ww p wn -K wa p g n F an -K aw p wn -K aa p g n      (5.80)
The stability and accuracy of the θ-method depends on the weighting parameter

θ = t θ -t n t n+1 -t n (5.81)
Particular integration schemes are recovered when θ takes specific values. θ = 0 corresponds to an explicit (Euler) integration, θ = 1 to an implicit (Backward Euler) integration and θ = 0.5 to a trapezoidal (Crank-Nicolson) integration. In θ-Stock, θ = 2/3. If θ 0.5 the process is unconditionally stable [START_REF] Booker | An investigation of the stability of numerical solutions of Biot's equations of consolidation[END_REF].

For mixed pressure-displacements problems, oscillations can appear in the initial pore pressure solution due to the large difference between the order of magnitude of permeability and mechanical rigidities coefficients, resulting in an ill-conditioned matrix.

The best solution to avoid numerical oscillations would be to choose the polynomial describing the pore pressure rates to be one order lower than the polynomial describing the velocities [START_REF] Zienkiewicz | The Finite Element Method: Its Basis and Fundamentals[END_REF].

Several stabilisation schemes have been proposed to permit equal-order interpolation for displacements and pore pressures fields. They usually consist in the addition of additional terms to the mass balance equation, modifying the incompressibility constraint. Recent works on this topics include the one of [START_REF] Mira | A new stabilized enhanced strain element with equal order of interpolation for soil consolidation problems[END_REF], [START_REF] White | Stabilized low-order finite elements for coupled soliddeformation/fluid-diffusion and their application to fault zone transients[END_REF], [START_REF] Aguilar | Numerical stabilization of Biot's consolidation model by a perturbation on the flow equation[END_REF], [START_REF] Preisig | Stabilization procedures in coupled poromechanics problems: A critical assessment[END_REF], ... We chose to use another option, which does not require strong modifications of the source code. [START_REF] Vermeer | An accuracy condition for consolidation by finite elements[END_REF], determined a lower limit for the time step, given a specific mesh size, ∆t 1 6

(∆h) 2 θC v (5.82) in which C v is the coefficient of consolidation (C v = Kw (1+ν) 3(1-ν) K+ φβp ρw γw
).

Equation ( 5.82) can be approximate by: ∆t 1 6

γ w (∆h) 2 θKK w (5.83)

Iterative resolution of equilibrium equations system

In order to solve equation 5.80, the total applied load is split into several loads increments which are considered to be applied over a period of time ∆t.

For each increment, an iterative algorithm based on the Newton-Raphson method is used to obtain a convergent solution. This section develops the iterative algorithm used to determine the response for one load increment. The smaller the increment size, the smaller the number of iterations needed to reach convergence.

The equation to solve is: K ∆dof = ∆F (5.84) in which

K =    R uu R uw R ua C wu C ww + θ∆tK ww C wa + θ∆tK wa C au C aw + θ∆tK aw C aa + θ∆tK aa    (5.85) ∆dof =      ∆u ∆p w ∆p g      (5.86) ∆F =      ∆F σ θ∆t∆F w θ∆t∆F a      + ∆t      0 F wn -K ww p wn -K wa p g n F an -K aw p wn -K aa p g n      (5.87)
The different steps of the iterative algorithm are described below and represented in figure 2. At iteration n, equation 5.84 is solved, with a rigidity matrix, K n-1 ( σ * n-1 , φ n-1 , s n-1 ), depending on the stress state at the beginning of the increment, and calculated considering an elastic loading (equations 3.109 and 3.110). The method used is therefore not the original Newton-Raphson method (which would consider elastoplastic loading), nor exactly the modified Newton-Raphson method (for which we would have

K n = K 0 ). K n-1 ∆dof n = Ψ 0 (5.89)
Since the problem is non-linear, the dofs values found by solving this system don't respect the mechanical equilibrium and iterations have to be performed to reach mechanical equilibrium.

3. The local stress-point algorithm (see section 5.6) is executed to find the real stress increment ∆σ n corresponding to the increments in strains ∆ε n , (deduced from ∆u n ), and pore pressures, ∆p ln and ∆p gn .

4. Knowing the updated stress state, the internal force increment is calculated:

∆F int = Ω B T ∆σ dΩ (5.90) 5.
The residual (which represents the unbalanced forces) is updated:

Ψ n = Ψ n-1 -∆F n (5.91) ∆F n = ∆F int K W A ∆dof n (5.

Stress point algorithm

The stress-point algorithm gives the increment in stress and plastic strains associated with the given increment of strains and pore pressures. The algorithmic procedure follows the one proposed by [START_REF] Ju | On energy-based coupled elastoplastic damage theories: constitutive modeling and computational aspects[END_REF] for damage-plastic couplings.

Input

Mechanical parameters

The following parameters are required as an input for the local algorithm:

❼ Elasticity parameters: κ, g, n, p r HAR (pressure dependent elasticity)

❼ Plasticity parameters: M , λ, r, β, p r BBM ❼ Retention properties: S r , α vg , m vg , n vg (Van Genuchten) ❼ Damage parameters: C 0 , C 1 , C 2

Initial state

The state at the start of the loading step is defined by the initial total stress, σ 0 , the total and elastic strains, ε 0 and ε e 0 (and thus implicitly the plastic strain), the fluid pore pressures, p g 0 and p l 0 (and therefore the suction s 0 = p g 0p l 0 ), and the values of the preconsolidation pressure, p 0 0 , and the damage parameter, d 0 , which act as hardening parameters.

Loading

As in the work of [START_REF] Sheng | Finite element formulation and algorithms for unsaturated soils. Part I: Theory[END_REF], pore pressures are treated as additional strain components, even though suction is similar to a stress variable in constitutive modelling.

The load input in then composed of the strain increment, ∆ε (deduced from the displacements increment given by the global algorithm), and the pore pressures increments, ∆p g and ∆p l .

The explicit scheme follows the following steps:

1. Calculation of the double effective stress 2. Elastic trial stress state (see section 5.6.2) 3. Plasticity correction (see section 5.6.3) 4. Damage correction (see section 5.6.4) 5. Calculation of the total stress

Trial elastic state

Since plastic and damage criteria are expressed in terms of the double effective stress, the first step is to calculate the initial double effective stress:

σ * 0 = σ + (-p g 0 + S l 0 (s 0 )s 0 )i 1 -d 0 (5.96)
The initial elasticity matrix is calculated with equation (2.56).

D e ( σ * 0 ) (5.97)

Since strain and pore pressures increment are fixed, their final values are:

ε 1 = ε 0 + ∆ε (5.98) p g 1 = p g 0 + ∆p g (5.99)
p l 1 = p l 0 + ∆p l (5.100)

s 1 = s 0 + ∆s (5.101)
The increment of the double effective stress, obtained by assuming that the material behaves elastically, is ∆ σ * = D e ( σ * 0 )∆ε (5.102) from which the trial elastic state is deduced.

σ * tr = σ * 0 + ∆ σ * (5.103) ε e tr = ε e 0 + ∆ε (5.104) It will next be checked if this stress state lies inside the plastic and damage surfaces (in this case the trial elastic state is the final one), or if dissipative phenomena are activated.

Plasticity

The plastic criterion is checked first, since plastic laws depend on suction and double effective stress increments, which are not affected by damage.

Check yield criterion

The yield criterion is checked for a stress state corresponding to the trial double effective stress, the final suction, and the initial preconsolidation pressure: If the stress state was previously lying inside the yield surface ( f p ( σ * 0 , p 0 0 , s 0 ) < 0), it is necessary to determine for which fraction of the total strain the soil behaves elastically before reaching the yield surface. This is equivalent to finding the scalar quantity, α, for which f p ( σ * 0 + αD e ∆ε, p * 0 0 , s 0 + α∆s) = 0 (5.109)

f p ( σ * tr , p 0 0 , s 1 ) = q2 tr -M 2 p * tr (p * c (p 0 0 , s 1 ) -p * tr ) (5.105) ❼ if f p F T OL
The starting stress state for stress integration will then be:

σ * e = σ * 0 + αD e ∆ε (5.110)

s e = s 0 + α∆s (5.111)
Various algorithms can be used to solve equation 5.109. In our case, we used the Pegasus algorithm from [START_REF] Dowell | The "Pegasus" method for computing the root of an equation[END_REF].

Stress integration

An explicit stress integration scheme is used to determine the increment in double effective stress, plastic strains and preconsolidation pressure.

Equations 3.95 to 3.97 give, (replacing fourth order tensors and second order tensors by the equivalent second order tensors and vectors)

∆ε p = a g a T f D e a T f D e a g -A ∆ε n + Ca g a T f D e a g -A ∆s n (5.112) ∆p 0 = Ba T f D e a T f D e a g -A ∆ε n + BC a T f D e a g -A ∆s n (5.113) ∆ σ * = D e i - a g a T f D e a T f D e a g -A ∆ε n + D e Ca g a T f D e a g -A ∆s n (5.114)
Since the problem is highly non-linear, the strain-suction increment is divided in smaller substeps: ∆ε n = ∆T ∆ε (5.115) ∆s n = ∆T ∆s (5.116) and the quantities, D e , a f , a g , A, B and C are updated at the beginning of each subincrement.

The size of the subincrements are adapted to keep the error below a certain tolerance, using the automatic substepping scheme developed by [START_REF] Sloan | Refined explicit integration of elastoplastic models with automatic error control[END_REF] for the saturated case, and extended by [START_REF] Sheng | Finite element formulation and algorithms for unsaturated soils. Part I: Theory[END_REF] to the unsaturated case.

Yield surface drift correction

At the end of stress integration, the stress state may diverge from the yield condition so that |f p | > F T OL. In that case it is necessary to correct the stress state to restore the stresses to the yield surface.

The method proposed by [START_REF] Potts | A critical assessment of methods of correcting for drift from the yield surface in elasto-plastic finite element analysis[END_REF], and later used by So lowski and [START_REF] So Lowski | Explicit stress integration with error control for the Barcelona Basic Model. Part I: Algorithms formulations[END_REF] for unsaturated soils, is described below.

The algorithm calculates a correction of the double effective stress state δ σ * while the increment of strains and suction are kept unchanged. This implies an associated change in the elastic and plastic strains partition:

δ σ * = D e δε e = -D e δε p = -δλD e ∂g p ∂ σ * (5.117) in which the scalar multiplier δλ is obtained from the consistency condition:

δλ = f p a T f D e a g -A (5.118)
The corrected values of the double effective stress, elastic and plastic strains and the hardening parameter are computed. If the stress still lies outside of the yield locus after the first correction, the procedure is repeated until the stress state is mapped back within the set tolerance.

At the end of the plasticity algorithm, the new values of the double effective stress, σ * 1 , the elastic and plastic strain, ε e 1 and ε p 1 , and the hardening parameter, p 0 1 , are known.

Damage

The updated value of the double effective stress, σ * 1 , serves as an input for the damage algorithm, which is checked in a second time.

Check damage criterion

The damage criterion is checked with the final value of the double effective stress and initial damage: 5.121) and the total stress state can be updated (section 5.6.5).

f d 1 = Y * d 1 -C 0 -C 1 d 0 (5.119) with Y * d 1 = Y * d ( σ * 1 ) = q1 -C 2 p * 1 (5.120) ❼ if f d 1 F T OL, there is no damage evolution, d 1 = d 0 (
❼ if f d 1 > F T OL, the damage algorithm is activated (go to step 2).

Intersection with damage criterion

As for the plastic algorithm, the Pegasus algorithm is used to determine which fraction α of the strain is applied before starting to create damage.

f d ( σ * 0 + D e ( σ * 0 )α∆ε e , d 0 ) = 0 (5.122)
Damage starts to increase for the stress state:

σ * c = σ * 0 + D e ( σ * 0 )α∆ε e (5.123)

Damage increment

The damage increment is calculated so that the final stress state lies on the damage surface:

∆d = ∆Y d C 1 = Y d ( σ * 1 ) -Y d ( σ * c ) C 1 (5.124)
The damage algorithms allows to update the value of the damage parameter:

d 1 = d 0 + ∆d (5.125)

Final state

Once all variables are updated, it possible to go back to the total stress state:

σ 1 = (1 -d 1 ) σ * 1 + p g 1 i -S l 1 (s 1 )s 1 i (5.126)
The total stress difference ∆σ = σ 1σ 0 is used in the global Newton-Raphson algorithm (equation 5.90).

The whole set of equations and numerical procedures needed to study fully coupled hydro-mechanical problems have been presented in this chapter. The following chapters consist in applications of the constitutive and numerical frameworks described previously to simulation of various boundary value problems. Chapter 6 makes use of the local stress-point algorithm to simulate purely mechanical problems at different water contents.

Chapter 7 shows the ability of this work to simulate complex problems in which fluid transfers, coupled with mechanics, play a major role.

CHAPTER 6

SIMULATION OF LABORATORY EXPERIMENTS

This chapter presents simulations of mechanical laboratory tests on two materials, Callovo-Oxfordian argillite and Boom Clay. Good agreement is achieved between experimental results and simulations. When no experimental data are available, it is shown that the mechanical behaviour follows expected trends. The argillite exhibits a coupled plasticdamage behaviour with a slight dependence on water content. The high sensibility of Boom clay on suction is well reproduced, with a ductile behaviour at low suctions, and a more brittle behaviour at high suctions.

Ce chapitre présente des résultats de simulation d'essais de laboratoire sur deux matériaux, l'argilite de l'Est et l'argile de Boom. Les résultats de simulations montrent une bonne adéquation avec les résultats expérimentaux. En l'absence de données expérimentales, le comportement mécanique suit les tendances attendues. Le comportement de l'argilite montre un couplage entre endommagement et plasticité, avec une faible dépendance à la teneur en eau. La forte influence de la succion sur le comportement de l'argile de Boom est bien reproduit, le comportement étant très plastique pour les faibles succions et plus à tendance fragile pour les fortes succions.

Callovo-Oxfordian Argillite

Due to their high mechanical strength and their low permeability, argillites are the more serious option to serve as a geological barrier in the context of radioactive waste disposal in France. Numerous research works have been funded in the recent years to study their thermo-hydro-mechanical properties, and many experimental results have been published in the literature. That is why we chose this material to validate our model. Figure 6.1 shows the typical stress-strain curve obtained from triaxial compression tests. It is clear that argillites exhibit plastic strains. Damage is revealed by elastic properties degradation as well as microstructural microcracks. This material is therefore adapted to be studied by a coupled damage-plastic model. The most extensive available experimental study on argillites comes from the PhD thesis of [START_REF] Chiarelli | Étude expérimentale et modélisation du comportement mécanique de l'argilite de l'est[END_REF], from which two papers have been derived [START_REF] Chiarelli | Modeling of elastoplastic damage behavior of a claystone[END_REF][START_REF] Shao | A coupled elastoplastic damage model for semi-brittle materials and extension to unsaturated conditions[END_REF]. She performed mechanical tests at different confining pressures and at different water contents.

Specific model for argillite

Mechanical properties

Hydrostatic tests (up to 60 MPa confining pressure) have shown that the bulk modulus can be considered pressure independent, and that the initial anisotropy is small (figure 6.2). We will therefore consider a linear elastic behaviour.

The Cam-clay model described in section 3.3 is not really adapted to model rocks mechanical behaviour. Figure 6.2 shows that during isotropic compression, the stress-strain behaviour remains elastic even for very high pressures. The preconsolidation pressure will therefore be set to a high value (at least superior to 60 MPa). The expression of the plastic potential needs to be modified to obtain a critical state conform to experiments (figure 6.3), i.e. to enable strain hardening for mean pressures well below the preconsolidation pressure. Fig. 6.2: Hydrostatic compression test [START_REF] Shao | A coupled elastoplastic damage model for semi-brittle materials and extension to unsaturated conditions[END_REF] Fig. 6.3: Initial yield surface, compressibility/dilatancy boundary and failure surface [START_REF] Shao | A coupled elastoplastic damage model for semi-brittle materials and extension to unsaturated conditions[END_REF] The yield surface is

f p = q2 -M 2 p * (p * c -p * ) (6.1)
The new plastic potential is chosen of the form:

g p = q2 -M 2 p * (η d p * c -p * ) (6.2)
Deriving equation 6.2 gives the equation of the critical state line (in the double effective stress space):

q = M p * 2 η d -1 (6.3)
For the sake of simplicity, and since no experimental data are available to calibrate the dependence of the yield surface on suction, the preconsolidation pressure is taken independent of suction in the total stress space and we simply have:

p * c = p 0 + S l s (6.4)
Elastic, plastic and damage parameters chosen in this study are summarised in table 6.1.

Elasticity Plasticity Damage K ν M λ 0 κ = 1/K η d p 0 C 0 C 1 C 2 GPa GPa -1 GPa -

Sl

Van Genuchten model Experimental data [START_REF] Hoxha | Saturated and unsaturated behaviour modelling of Meuse-Haute/Marne argillite[END_REF] Fig. 6.4: Water retention curve α vg (MPa -1 ) n vg m vg S r 0.04 1.5 0.55 0 Table 6.2: Retention parameters Mechanical tests are carried on for 3 ranges of water content. Table 6.3 gives the saturation degree (given by [START_REF] Chiarelli | Modélisation élastoplastique couplée à l'endommagement anisotrope induit pour des argilites[END_REF]) and the suction value (calculated with Van Genuchten equation) corresponding to each water content range.

w (%) S l s (MPa) 4-5 0.61 32 5-6 0.71 23 7-8 0.93 7 Table 6.3: Hydraulic state for three water content ranges

Simulation of mechanical laboratory tests

Tests at natural water content

Triaxial compression tests, proportional compression tests and lateral extension tests are simulated. All the experimental data comes from the PhD thesis of [START_REF] Chiarelli | Étude expérimentale et modélisation du comportement mécanique de l'argilite de l'est[END_REF]. The tests are executed on samples at their natural water content, which is in the range w = 5-6%. Suction is therefore taken equal to 23 MPa.

Triaxial compression tests:

Triaxial compression tests with unloading-reloading cycles are simulated for four confining pressures (0 MPa, 5 MPa, 10 MPa and 20MPa). The comparison of computed stress-strain curves (a) and volumetric strains (b) with experimental data are represented on figures 6.5 to 6.8.

The agreement between computed and experimental stress-strain curves is quite good.

The plot of the volumetric strains versus axial strains shows that the volumetric strains are systematically underestimated. This could be improved by modifiying the plastic potential in order to increase the ratio εv / εs . The transition from volumetric compressibility to dilatancy (shown in figure 6.6.b) can not be represented by the plasticity model since it doesn't allow the stress path to cross the critical state line.

-0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 q (MPa) .9 shows the stress-strain curves obtained for the different confining pressures.

ε a (%) ε r (%)
The corresponding effective stress paths can be seen on figure 6.10. It can be seen that, since the initial elastic domain is small and the onset of plastic yielding is similar for every cases, the main parameter controlling the influence of initial confining pressure on the stress-strain behaviour is the position of the critical state line. The evolution of damage with the total deviatoric stress is given in figure 6.11. The confining pressure determines for which value of the deviatoric stress damage initiates. Afterwards, damage evolution is driven by the increase in deviatoric stress, independently of the value of the confining pressure.

The resulting degradation of elastic properties with radial strains is given in figure 6.12.

The experimental Young's (E L ) and radial (E R ) moduli have been evaluated from the slopes of unloading-reloading cycles [START_REF] Chiarelli | Modeling of elastoplastic damage behavior of a claystone[END_REF].

Since we consider a constant Poisson's ratio, the degradation of Young's and radial moduli are the same in our model: The computed Young's modulus degradation is in the same value range as the measured one and follows the same trend. However, due to the important scatter in experimental results, it can not really serve as a quantitative validation. Concerning the radial modulus variation, the computed values are smaller than the experimental data. The shape of the curve remains however similar. It appears that radial and axial stiffnesses evolve differently during the triaxial test. Our damage model should therefore be refined, either by considering an anisotropic damage variable, or, as do [START_REF] Shao | A coupled elastoplastic damage model for semi-brittle materials and extension to unsaturated conditions[END_REF], by considering that hydrostatic and deviatoric stresses are affected differently by damage.

R = E L E L 0 = E R E R 0 .

Proportional compression tests:

The simulation of three proportional loading tests are represented on figure 6.13. In these tests, axial and lateral stresses are simultaneously increased with a constant ratio (σ 3 = kσ 1 ). q (MPa) Lateral extension tests: Figure 6.14 shows the simulation of lateral extension tests for two initial confinements of 30 MPa and 60 MPa. The lateral stress σ 1 is then decreased while keeping the axial stress σ 3 constant. Experimental data and numerical simulations results are in good agreement at the begining of the test, and start to diverge when the lateral stress reaches one third of its initial value. The effective stress path is given is figure 6.15. At the end of the test, the stress state lies very close to the critical state line, which is most likely the cause of the high plastic strains observed. Figure 6.16 gives the evolution of damage as a function of the lateral stress. Good agreement with this test is as important as with triaxial compression test, since the stress path followed during tunnel excavation is closer to the one followed during a lateral extension test than to the one due to a triaxial compression test.

ε a (%) ε r (%) numerical results experimental data c.
Effect of water content on mechanical behaviour Figure 6.17 shows the stress-strain curves simulated for three water contents for a confining pressure of 20 MPa. Although numerical results do not agree perfectly with experimental data, the effect of suction on the mechanical behaviour is similar. Suction acts on the material by adding the component S l s to the confining pressure. q (MPa)

ε a (%) ε r (%)
w=4-5 % w=5-6 % w=7-8 % w=4-5 % w=5-6 % w=7-8 % Fig. 6.17: Stress-strain curves for different water contents (confining pressure: 20 MPa)

Effect of damage hardening parameter

During a triaxial test, the rate of damage evolution is given by ∆d

= 1 - C 2 3 ∆q C 1 (6.5)
C 1 is therefore the parameter controlling damage evolution rate. The lower C 1 , the higher the increment of damage for a given stress increment.

The simulation of the triaxial compression test for a confinement of 20 MPa has been repeated for three values of C 1 (100,150,200). The stress-strain curves are represented in figure 6.18, the evolution of damage with deviatoric stress in figure 6.19 and the degradation of the bulk modulus in figure 6.20. q (MPa) Choosing a low value for C 1 enables the simulation of a quicker rigidity degradation with radial strains. However, when damage reaches 50 %, the behaviour switches from an apparent hardening behaviour to an apparent softening behaviour, which under increasing deviatoric stress leads to failure. In the case of linear elasticity, with the principle of strain equivalence, the final apparent elastic moduli just before failure are equal to 50 % of their initial values. With the hypothesis of equivalent elastic energy, for 50 % of damage, the ratio R would be equal to R = 0.5 2 = 25%. However, the coupling between damage and plasticity would necessitate another strategy.

ε a (%) ε r (%) C 1 =200 MPa C 1 =150 MPa C 1 =100 MPa
C 1 =200 MPa C 1 =150 MPa C 1 =100 MPa data E R /E R0 data E L /E L0
At this point, the limited amount of experimental data and their scattering do not allow to conclude whether or not the moduli degradation obtained with the principle of strain equivalence is satisfactory. It should also be highlighted that the damage hardening parameter C 1 has an impact on both the rigidities degradation and the stress-strain behaviour. A compromise has to be made when choosing C 1 to correctly reproduce both.

Boom Clay

Boom Clay has been selected as a possible host rock for deep radioactive waste disposal in Belgium. It has therefore been extensively studied either from experiments on undisturbed natural samples, or on samples prepared by compaction from Boom clay powder. Many experimental data on saturated natural Boom clay are available in the literature [START_REF] Horseman | Geotechnical characterization of boom clay in relation to the disposal of radioactive waste[END_REF][START_REF] Baldi | Developments in modelling of thermohydro-geomechanical behaviour of Boom clay and clay-based buffer materials[END_REF][START_REF] Coll | Endommagement des roches argileuses et perméabilité induite au voisinage d'ouvrages souterrains[END_REF][START_REF] Sultan | Yielding and plastic behaviour of Boom clay[END_REF]. However, concerning the unsaturated behaviour, most of the studies have been made on compacted [START_REF] Bernier | Suction-controlled experiments on Boom clay[END_REF][START_REF] Romero | Characterisation and thermo-hydro-mechanical behaviour of unsaturated Boom clay: an experimental study[END_REF] or remoulded [START_REF] Al-Mukhtar | The fabric of a clay soil under controlled mechanical and hydraulic stress states[END_REF] samples, and only a few on undisturbed samples [START_REF] Cui | Suction effects in deep Boom Clay block samples[END_REF][START_REF] Vecchia | Some remarks on the hydromechanical constitutive modelling of natural and compacted Boom clay[END_REF]. Moreover, these mechanical tests at different suctions are limited to oedometer and isotropic compression tests. By comparing experiments on natural and compacted samples, Della [START_REF] Vecchia | Some remarks on the hydromechanical constitutive modelling of natural and compacted Boom clay[END_REF] concluded that the same constitutive framework seems to be applicable to natural Boom clay and to the material compacted from the clay powder. However, the mechanical parameters have to be adapted to reflect the different microstructure.

Boom clay is considered as an overconsolidated plastic clay. Although no data are available on rigidity degradation, the existence of an Excavation damaged Zone around excavations has been revealed. Excavation induced fractures are observed around galleries [START_REF] Bastiaens | Experiments and conclusions on fracturing, self-healing and self-sealing processes in clays[END_REF][START_REF] Bernier | Twenty-five years' geotechnical observation and testing in the Tertiary Boom Clay formation[END_REF][START_REF] Van Marcke | Excavation induced fractures in a plastic clay formation: Observations at the HADES URF[END_REF]. Damage has also been assessed by seismic [START_REF] Bastiaens | Experiments and conclusions on fracturing, self-healing and self-sealing processes in clays[END_REF] and acoustic [START_REF] Lavrov | Acoustic emission in host-rock material for radioactive waste disposal: comparison between clay and rock salt[END_REF] measurements. Boom clay can exhibit both ductile and brittle behaviours [START_REF] Dehandschutter | Brittle fractures and ductile shear bands in argillaceous sediments: inferences from Oligocene Boom Clay (Belgium)[END_REF].

The transition between the failure modes depends strongly on the confining pressure and is also influenced by the water content [START_REF] Al-Shayea | The combined effect of clay and moisture content on the behavior of remolded unsaturated soils[END_REF] and by the overconsolidation ratio.

Since all the available data comes from samples of diverse provenance, and that Boom clay parameters vary depending on the coring depth and the mineralogy, we will choose parameters in the range of the ones found in the literature and simulate typical experiments to show that our model is able to reproduce the main trends. However, we do not expect perfect agreement with experimental data coming from different laboratories with a single set of parameters.

Model parameters

Mechanical properties

Elastic, plastic and damage parameters chosen in this study are summarised in table 6.4. The pressure exponent is taken equal to 0.6. Due to the modifications of volumetric compression laws (equations 2.24 and 3.60) to incorporated hyperelasticity with an exponent non equal to 1, the parameters κ and λ take values different from the one usually found in the literature. It also has an impact on the values of the parameters of the Loading-Collapse curve. Water retention properties are chosen to obtain an average of the drying and wetting curves given by Della [START_REF] Vecchia | Some remarks on the hydromechanical constitutive modelling of natural and compacted Boom clay[END_REF]. Van Genuchten parameters are given in table 6.5.

Elasticity

The water retention curve as well as the suction stress component (S l s) are represented in figure 6.21.

α vg (MPa -1 ) n vg m vg S r 0.28 2.3 0.21 0

Simulation of mechanical laboratory tests

In the following examples, the preconsolidation pressure is taken equal to 5.5 MPa.

Isotropic compression test: Figure 6.22 shows the simulation of an isotropic compression test with comparison to experimental results from [START_REF] Baldi | Developments in modelling of thermohydro-geomechanical behaviour of Boom clay and clay-based buffer materials[END_REF]. The tests consist in drained isotropic cycles starting from 2 MPa of effective isotropic stress up to a maximum of 8 MPa. Back pressure value is 2 MPa. Two different measures of volumetric strains are made, one by calculating three times the vertical strain (ε v = 3ε 1 ), the other one obtained from water displaced from the specimen (ε v = ∆V w ).

The different volumetric strains given by the two types of measures makes it difficult to draw conclusions about the quality of our simulations. Oedometer tests: Figure 6.23 shows the simulation of an oedometer test in saturated conditions with comparison to the experimental data of [START_REF] Horseman | Geotechnical characterization of boom clay in relation to the disposal of radioactive waste[END_REF]. The agreement is quite good, although the slope of unloading-reloading curves shows that the loss of rigidity is underestimated. Figure 6.24 shows the simulation of a test which experimental results come from Della [START_REF] Vecchia | Some remarks on the hydromechanical constitutive modelling of natural and compacted Boom clay[END_REF]. A soil sample, initially at a suction of 3 MPa, is compressed under undrained conditions up to a vertical stress of 3 MPa, and later put in contact with water. The sample is then compressed up to 8 MPa, under drained conditions, and finally unloaded. We can see that numerical results give no volume change during the undrained loading part, contrary to the experimental results. It may be due to water leakage during the experiment, otherwise it means that we could improve our model by considering soil particles compressibility (or a higher water compressibility). The main behaviour trends are otherwise well represented, the swelling occuring during unloading being present, although underestimated.

Drained triaxial compression tests: Triaxial drained compression tests with unloadingreloading cycles are simulated for three confining pressures (2 MPa,3 MPa,4 MPa). The comparison of computed stress-strain curves (a) and volumetric strains (b) with experimental data are represented on figures 6.25 to 6.27.

The agreement between computed and experimental stress-strain curves is quite good. However, the model does not capture well the smooth transition between elastic and plastic behaviour. This behaviour could be improved by using more advanced versions of the Cam-Clay model, such as bounding surface plasticity [START_REF] Dafalias | Bounding Surface Plasticity. I: Mathematical Foundation and Hypoplasticity[END_REF] or continuous hyperplasticity [START_REF] Puzrin | Fundamentals of kinematic hardening hyperplasticity[END_REF].

The plot of the volumetric strains versus axial strains shows a satisfactory agreement between experimental and numerical results. Since damage tends to predict large volumetric strains compared to experimental data, the non-associativity parameter ζ has been chosen equal to 1.5 to increase the deviatoric strains predicted by the plasticity model. Figure 6.28 shows the stress-strain curves obtained for the different confining pressures.

The corresponding effective stress paths can be seen on figure 6.29. The evolution of damage with the total deviatoric stress is given in figure 6.30. It can be seen that failure is mainly due to the plasticity phenomenon, since the effective stress path reaches the critical state before damage attains the critical 50 % value. When the confining pressure increases, it is possible to get more plastic strains before failure, the behaviour is more ductile. For low confining pressures, although the damage value is smaller, the combination of damage with the critical state model leads to less plastic strains before failure, i.e. a more brittle behaviour. The resulting degradation of elastic properties with radial strains is given in figure 6.31. The predicted rigidity loss is low, around 10 % at its maximum. Due the the lack of experimental data, it is hard to say whether the principle of strain equivalence is acceptable in this case. However, the slopes of the unloading-reloading curves in figures 6.25 to 6.27 seem to be in good agreement with experiments. It is interesting to note that equation 2.38 gives a maximum possible loss of rigidity, for 50 % of damage, equal to Fig. 6.31: Evolution of apparent bulk modulus with axial strains

1 -R = 1 -(1 -d) (1-n) = 1 -0.5 0.4 = 24, 2%.

Suction effects on mechanical behaviour

Although no experimental data have been found in the literature about Boom clay, [START_REF] Al-Shayea | The combined effect of clay and moisture content on the behavior of remolded unsaturated soils[END_REF] showed that materials with high clay content exhibit a ductile/brittle behaviour transition when their water content decreases (see figure 6.32). Ductile behaviour is characterised by the ability to sustain large plastic strains during plastic hardening. Brittle behaviour is characterised by abrupt failure at a well-defined peak strength with strong softening. Figure 6.32 also shows higher shear strength for low water contents.

Simulations of triaxial compression tests are made for different constant suctions (0 MPa, 2 MPa and 6 MPa). The shape of the LC curve corresponding to the chosen parameters is given in figure 6.33.

The stress-strain curves obtained for different suctions are given in figure 6.34. The corresponding effective stress paths can be seen on figure 6.35. The evolution of damage with deviatoric stress is given in figure 6.36.

The ductile/brittle transition linked to the increase of suction is well reproduced. At low suction, large plastic strains occur and damage remains low. At higher suctions, shear strength is higher and the peak stress is attained for lower plastic strains and higher values of the damage parameter. This is due, partly to the shift in the stress path due to the increase in effective stress, and also to the increase of the elastic domain with suction. CHAPTER 7

APPLICATION TO FULLY-COUPLED HYDRO-MECHANICAL PROBLEMS

This chapter presents illustrative applications of the constitutive modelling framework to fully-coupled hydro-mechanical problems. First, damage induced by hydric loading is simulated, assuming that it is due to differential shrinkage in low permeable media. Then, plastic and damage zones induced by excavation are studied, first at different suctions under drained conditions, then under undrained conditions.

Ce chapitre présente des applications illustrant les capacités du modèle proposé. Premièrement, l'endommagement induit par un chargement hydrique est simulé, en supposant que celui-ci est dû à la contraction hétérogène des milieux peu perméables. Ensuite, les zones plastiques et endommagée dues à l'excavation sont étudiées, premièrement à différentes succions en considérant des conditions drainées, puis en conditions non drainées. Accurate excavation modelling is of a great importance in contexts such as tunnelling, mining, or radioactive waste disposal. Coupled hydro-mechanical phenomena have to be taken into account during the different phases of the excavation process. Due to the stress redistribution and the resulting convergence during the initial excavation stage, micro-and macro-fracturing is likely to appear in the surrounding rock. Depending on the type of rock and its permeability, either drained or undrained behaviours can be considered. The following tunnel ventilation dehydrates the rock around the opening. This desaturation tends to strengthen the rock, but, since strains are constrained, can also cause induced fractures and micro-cracks. Although beyond the scope of this thesis, other processes such as heating and long term behaviour (creep, self-sealing, self-healing) have to be taken into account in the context of radioactive waste disposal. A summary of the different processes to consider for different rock types is presented in the paper of [START_REF] Tsang | Geohydromechanical processes in the Excavation Damaged Zone in crystalline rock, rock salt, and indurated and plastic clays-in the context of radioactive waste disposal[END_REF]. Since a validation of our model behaviour under mechanical loading has been made in the previous chapter, the first section of this chapter will deal with purely hydric loading.

Damage induced by hydric loading

Under natural environmental conditions, evaporation tends to remove water from geomaterials, such as soil or concrete, and thus makes them shrink. If deformation is constrained, cracks may therefore appear.

Understanding how cracks and damage are created in drying geomaterials is of interest in geotechnical engineering applications, in agricultural science and in concrete technology. Although many research works have proposed explanations for the formation of these cracks, no consensus has yet been reached. The work presented here deals with damage, i.e. micro-cracks. However, the mechanisms of damage initiation are similar to the ones leading to macro-cracks.

Desiccating soils with constrained deformations leads to tensile total stresses. The most common approach is therefore to assume a tensile strength determining the onset of cracking. However, the constitutive model presented previously is based on the assumption that the damage criterion is expressed in terms of the double effective stress. This double effective stress, when suction increases, always remains positive, even under constrained conditions, and desiccation results in an increase of the compressive stress applied on the solid skeleton. A possible alternative could be to consider a tensile strength dependent on suction and becoming positive when suction increases. [START_REF] Péron | Desiccation cracking of soils[END_REF] proposed the following form for the tensile strength as a function of suction and the tensile strength in saturated conditions:

σ * t = σ sat t + k 2 1 -exp(- k 1 s k 2 ) .
Another process advanced to explain desiccation cracks is the air invasion of the largest pores. The air-water interface membrane then pushes particles away from the invaded pore, which creates a crack [START_REF] Shin | Desiccation cracks in saturated fine-grained soils: particlelevel phenomena and effective-stress analysis[END_REF]. Others highlight the role of heterogeneities in cracks initiation [START_REF] Lagier | Numerical strategies for prediction of drying cracks in heterogeneous materials: Comparison upon experimental results[END_REF]. He observed that the onset of microcracking occurs during the loading stage and is strongly controlled by humidification or desiccation rate. This microcracking seems therefore to be linked to moisture gradients induced by fast humidification/desiccation when the loading rate doesn't allow pressures to dissipate. Damage due to humidification has also been observed by [START_REF] Valès | Modes de déformation et d'endommagement de roches argileuses profondes sous sollicitations hydro-mécaniques[END_REF].

In the following, we will study how moisture gradients affect the stress state in an homogeneous porous media and may lead to damage initiation. Studies considering damage due to differential shrinkage usually find that microcracks are localised at the surface in contact with the atmosphere [START_REF] De Sa | Analysis of microcracking induced by differential drying shrinkage[END_REF][START_REF] Jia | Elastoplastic damage modeling of desaturation and resaturation in argillites[END_REF].

Material parameters

The mechanical parameters are the ones of the Callovo-Oxfordian argillite studied in chapter 6.1. A water pressure will be applied on surfaces BC and CD, the normal total stress remaining constant. The air pressure takes a constant null value in the whole specimen.

Intrinsic permeability K int (φ) = K w0 φ 3 (1 -φ) 2 (1 -φ 0 ) 2 φ 3 0 K w0 = 1.10 -15 m.s -1 φ 0 = 15% Relative permeability k r (S l ) = S l -S r 1 -S r m m = 3

Desiccation

The initial suction is set equal to 2 MPa and increased up to 20 MPa at the sample surface at a rate of 0.2 MPa/min (which corresponds approximatively to 0.1% RH/min).

Figure 7.5 shows the water pressure field at four time steps. Time t=1h30 corresponds to the moment at which the total increment of suction has been applied. It can be noticed that, due the high loading rate combined with the low permeability of the material, fluid pressures have not yet had time to equilibrate at that time. The likelihood of induced damage is therefore at its maximum for the given suction increase. We can even notice the appearance of overpressures at the sample core, certainly due to the compression exerted by the external ring which already started to desaturate. The spatial evolution of the total stress components at this time step along line AB are given in figure 7.6. It is verified that this pressure field induces compression at the sample core and traction at its surface. Water evaporates from the soil sample during the following time steps, until pressure equilibrium is reached around t=10h. The damage pattern at the end of the test is given in figure 7.7. Most of the damage is concentrated in the few millimetres near the surface of the sample with its highest values far from the sample free extremities. Damage is maximum at the center of the top surface, for which strain constraints are higher. However, the sample core also suffers non negligible damage, which is due to the overpressures appearing at the beginning of the test.

In the following, we will study point B to highlight the principal mechanisms leading to damage in our model. The total and effective stress paths followed at point B are represented in figure 7.9. Although the total mean stress is tensile, the double effective mean stress remains positive due to compressive suction effects. Since we considered all the dissipative processes to be driven by the double effective stress, damage can not be attributed to tensile stresses in the solid matrix within our modelling framework. However, figure 7.9 shows that the increase in deviatoric stress makes the stress path reach the damage criterion, and therefore induces damage. The damage pattern at the end of the test is given in figure 7.11. As for the desiccation test, damage is localised near the surface, far from the sample extremities. However, in this case, the sample core remains intact. decreases, the solid skeleton remaining in compression. As the applied suction decreases, the deviatoric stress increases due to the heterogeneous pressure field, reaches a maximum value, then decreases when pore pressures homogenise. Damage is not initiated from the beginning of loading, but its evolution afterwards is much more abrupt than for desiccation. The total and effective stress paths followed at point B are represented in figure 7.13. In this case, even total stresses are compressive. The simulated damage is again linked with the deviatoric stress increase due to differential shrinkage. These examples show that, within a constitutive framework based on the concept of an effective stress with a component related to suction, this double effective stress always remains positive, even when shrinkage is constrained and tensile total stresses develop. Indeed, the compressive effect of suction is more important than the effect of total tensile stresses due to constrained shrinkage. It is therefore, under this assumption, impossible to model damage initiation by considering it caused by traction of the solid skeleton. However, although this interpretation may need further discussion and experimental evidence, considering that damage is due to deviatoric stresses in the solid skeleton enables the reproduction of damage due to inhomogeneous deformation during desiccation and humidification. This outcome is of interest since none of the existing approaches is capable of modelling damage due to humidification. This approach, although not demonstrated here, is also sensitive to the rate at which suction increments are applied. Indeed, if loading is slow enough, pressures homogenise during suction application and no deviatoric stress is created.

Excavation modelling

Works about excavation, ventilation and heating of deep galleries can be found in [START_REF] Barnichon | Contribution of the bounding surface plasticity to the simulation of gallery excavation in plastic clays[END_REF]; [START_REF] Barnichon | Observations and predictions of hydromechanical coupling effects in the Boom clay[END_REF]; [START_REF] Hoxha | Poroplastic modelling of the excavation and ventilation of a deep cavity[END_REF]; [START_REF] Gatmiri | Effect of excavation on the thermo-hydro-mechanical behaviour of a geological barrier[END_REF]; [START_REF] Arson | Thermo-hydro-mechanical modeling of damage in unsaturated porous media: Theoretical framework and numerical study of the EDZ[END_REF] and [START_REF] Garitte | Analysis of hydro-mechanical processes in a ventilated tunnel in an argillaceous rock on the basis of different modelling approaches[END_REF] 7.2.1 Impact of suction on the extent and properties of the excavation damaged zone

The first study consists in simulation of excavation in Boom clay at different suctions.

Although the problem is simplified, assuming drained conditions and a simple model for Boom clay (linear elasticity), this example illustrates well how the saturation state around the tunnel affects the development of damage and plastic zones. This work has been published in the proceedings of the fifth Biot conference on Poromechanics [START_REF] Le Pense | Simulation of the Unsaturated Excavation Damage Zone around a Tunnel Using a Fully Coupled Damage-Plasticity Model[END_REF].

The specific equations chosen to represent clay behaviour are summarised in Table 7.4. Elasticity is assumed to be linear. Plasticity is modelled using the Barcelona Basic Model expressed in terms of the double effective stress.

Excavation Modelling

We simulate the excavation of a tunnel of radius R in Boom clay (note that in this conceptual study, the numerical value of R does not need to be specified). Excavation is simulated at different constant suctions, under the assumption that fluid transfers are fast enough to consider drained conditions.

Initial conditions and material parameters are taken of the same order of magnitude as the ones found in the literature for Boom clay at the HADES underground laboratory at Mol, Belgium [START_REF] Romero | Characterisation and thermo-hydro-mechanical behaviour of unsaturated Boom clay: an experimental study[END_REF][START_REF] Vecchia | Some remarks on the hydromechanical constitutive modelling of natural and compacted Boom clay[END_REF].

Plasticity

Barcelona Basic Model [START_REF] Alonso | A constitutive model for partially saturated soils[END_REF])

f p ( σ * , p * 0 , s) = q * 2 -M 2 p * (p * c (p * 0 , s) -p * ) g p ( σ * , p * 0 , s) = ζ q * 2 -M 2 p * (p * c (p * 0 , s) -p * ) p * c = p r p * 0 p r λ 0 -κ λ(s)-κ + S l s λ(s) = λ 0 [(1 -r) exp(-βs) + r] ṗ * 0 = p * 0 λ 0 -κ εp v Damage f d ( σ * , d) = q * -C 2 p * -C 0 -C 1 d ḋ = 1 C 1 - C 2 3 I d + 3 σ * d 2q * : σ *

Initial and boundary conditions

The initial stress state is supposed to be isotropic : σ r = σ θ = σ z = 4.5 MPa. The preconsolidation pressure is taken equal to p * 0 = 6 MPa. Under the assumption that the initial stress state is isotropic, the problem can be modelled by an axisymmetric analysis (see Figure 7.14). Displacements along e z direction are set equal to zero on surfaces AD and BC.

On surface CD, the stress state is supposed not to be disturbed by the excavation, ∆σ r = 0.

On surface AB, radial stress is supposed to be equal to zero after excavation, i.e. during the simulation, a stress difference ∆σ r = -4.5 MPa is applied.

Suction is kept constant at each node of the mesh during simulation.

Results and discussion

Stress path at the tunnel wall To understand better how our model handles the transition between plastic and brittle behaviour, let see the stress paths followed during excavation at tunnel wall. For two different values of suction (s=0 and s=1Mpa), figures 7.15 and 7.16 show the double effective stress path (thick solid line), the total stress path (thick dashed line) and plasticity (thin solid line) and damage (thin dashed line) criteria. These stress paths illustrate the main mechanism allowing for transition from ductile to brittle behaviour, which is the competition between the two dissipative phenomena. Indeed, for a saturated soil, the stress path reaches the plastic yield surface before reaching the damage criterion. The occurrence of plasticity modifies the stress path and prevents the development of high damage. For high suctions, the elastic domain is much larger, due to the dependence of the yield surface on suction, and the stress path will reach the damage criterion first, which leads to a more brittle behaviour.

Both the double effective stress path ( Ã * B * C * ) and the total stress path (ABC) are represented in Figure 7.16. Before damage is initiated (AB) the difference between the two stresses is only due to suction, which increases the confining pressure on the soil skeleton. After damage initiation (BC), the stress acting on soil skeleton increases due to the appearance of microcracks, while total stress tends to decrease.

Stress spatial evolution Figures 7.17 For the saturated case, plasticity is the main dissipative phenomenon and leads to a decompression of tangential stress around the opening.

In the unsaturated case, there is also a decompression of total stress but the value of double effective stress is much higher.

Development of plastic and damage zones

The spatial evolution of damage and plasticity (represented by the norm of the plastic strain matrix |ε p | = √ ε p : ε p ) are given in figures 7.19 and 7.20. The amplitude of plastic strains and the extent of the plastic zone both decrease with suction.

Damage at the tunnel wall increases when suction is higher, which is consistent with the model formulation, which was aimed to capture the shift from ductile to brittle behaviour associated with suction changes.

Damage develops further away from the excavation as suction is lower and as the plastic zone is wider. This observation is assumed to be due to the choice of constitutive models used for plastic and damage criteria, as explained from the stress paths illustrated in Figures 7.15 and 7.16. Indeed, Figure 7.15 shows that the occurrence of plasticity substantially modifies the stress path which contributes to the development of damage in the plastic zone.

This study allowed to show the effect of suction on the extent of damage and plastic zones created by excavation. When considering a drained behaviour, numerical results predict a wide plastic zone and high values of plastic strains when the ambient suction applied is low. For low suctions, damage remains low, but develops in the whole plastic zone, as opposed to states of high suction, in which damage is higher but localised in a narrower zone.

Simulation of excavation under undrained conditions

A more realistic excavation simulation is presented in this section. The excavation is assumed to be realised fast enough to take place under undrained conditions.

Geometry and boundary conditions

The initial stress state is supposed to be isotropic : σ r = σ θ = σ z = 4.5 MPa. The preconsolidation pressure is taken equal to p * 0 = 5.5 MPa. The initial pore water pressure is taken equal to p l = 2 MPa.

Under the assumption that the initial stress state is isotropic, the problem can be considered unidimensional and be modelled by an axisymmetric analysis (see figure 7.21). The excavation radius is taken equal to R = 1 m, and the spatial variable considered is the dimensionless radius, r/R. The problem is meshed from the tunnel wall, r/R = 1, to a distance r/R = 20, distance at which the stress state is not affected by the excavation (see figure 7.22). Displacements along e z direction are set equal to zero on surfaces AD and BC.

On surface CD, the stress state is supposed not to be disturbed by the excavation, ∆σ r = 0 and p w = 2 MPa.

On surface AB, the radial stress is supposed to be equal to zero after excavation, i.e. during the simulation, a stress difference ∆σ r = -4.5 MPa is applied. The stress release is supposed to be linear with time during excavation, which is supposed to last 15 days. Water flux is prevented across this section, Q w = 0.

The air pressure takes a constant null value in the whole specimen.

Material parameters

The mechanical parameters are the ones of Boom clay studied in chapter 6.2. They are summarised in table 7.7.

Fluid transfers are taken into account in this study. The permeability functions chosen as well as their parameters values are given in table 7.9.

The initial porosity is taken equal to 39 %. (1φ 0 ) 2 φ 3 0 K w0 = 1.10 -13 m.s -1 φ 0 = 39%

Elasticity

Relative permeability k r (S l ) = S l -S r 1 -S r m m = 3 Table 7.9: Hydraulic functions and parameters

Results

Figure 7.22 shows the spatial distribution of the different components of the total stress tensor as well as the water pore pressure at the end of excavation. It can first be noticed that the stress state is not significantly affected by the excavation from a distance to the tunnel r/R = 15 and thus that the choice of the outer boundary, at r/R = 20, is acceptable to consider that the response is similar to the one of an infinite medium. The stress release induces a zone in which water pressure is strongly reduced, which corresponds to the plastic zone. Negative pore pressures are generated close to the tunnel wall, what has been also observed at HADES research laboratory [START_REF] Bernier | Twenty-five years' geotechnical observation and testing in the Tertiary Boom Clay formation[END_REF].

The spatial evolution of radial displacements is given in figure 7.23. A radial convergence of 9.3 cm is obtained. Although the excavation radius is different, it is worth noting that this value lies in the same range as the radial settlements measured in Boom clay (9 cm for a 4 m diameter gallery with lining [START_REF] Bernier | Twenty-five years' geotechnical observation and testing in the Tertiary Boom Clay formation[END_REF], 14.5 cm for a 4 m radius gallery with lining [START_REF] Barnichon | Observations and predictions of hydromechanical coupling effects in the Boom clay[END_REF]). The simulated damage zone extends further than the plastic zone. The plastic strain evolution curve has a regular shape. Damage, however, takes a constant value from r/R = 2 to r/R = 4. To explain this surprising result, the time evolution of damage, plastic strains, and water pressure at the tunnel wall is represented in figure 7.25. There are three different periods during which different mechanisms take place. Damage is initiated first and water pressure remains unchanged as long as the behaviour remains elastic. When plasticity is activated, water pore pressure starts decreasing and damage stops evolving. Once water pressure reaches negative values, damage starts increasing again. The stress path, represented on figure 7.26, shows that this behaviour is due to the overconsolidated state of the clay, which implies that, once the stress path reaches the yield surface, only the extension of the plastic yield surface with suction can allow further increase of damage. Although these results should be validated against experimental data for other values of the tunnel radius, the values obtained for pore pressures and displacements are realistic. This example illustrates well how our modelling framework allows to simulate both damage and plastic strains. A major improvement would be to incorporate permeability dependence on damage, and study how it affects the overall behaviour.

CONCLUSIONS AND PERSPECTIVES

Conclusions

A constitutive modelling framework allowing for damage-plasticity couplings in unsaturated porous media has been proposed. This framework is based on the assumption of a double effective stress, accounting for damage and suction effects, which controls the material mechanical behaviour. The relationship between suction and the degree of saturation has proved to have a great influence on the value of the double effective stress. Water retention properties have therefore been studied in details and a proposition has been made to incorporate porosity into the water retention curve expression, in order to meet thermodynamical consistency requirements. Although a simple expression has been chosen for the double effective stress, it has been proved to be able to reproduce the main features of damaged unsaturated soils.

An hyperelastic formulation has been chosen to represent pressure dependent elasticity. Energy conservation in the elastic domain is therefore ensured. The two main hypotheses available to incorporate damage into elasticity and plasticity equations have been compared. It has been shown that the principle of equivalent elastic energy predicts a higher rigidity degradation than the principle of strain equivalence for the same value of the damage parameter. The pressure exponent in the hyperelastic formulation also proved to have a strong influence on rigidity degradation with damage. However, too few experimental damage and rigidity measurements are available to choose one hypothesis against the other.

The principle of strain equivalence has therefore been chosen for its ability to provide a straightforward way of coupling damage and plasticity. Damage and suction effects are taken into account by replacing the total stress by the double effective stress into elasticity and plasticity equations, which means that damage and plasticity criteria and evolution laws are expressed in terms of the double effective stress. This allows for a direct dependence of damage and plasticity criteria on suction and damage in the total stress space. Any plasticity model could be accommodated into this framework. The Cam-clay model has been used as an example, with an extra-dependence of the plastic yield criterion on suction, in order to account for the high sensitivity of clay on water content.

A specific local algorithm has been developed, to determine the stress increment associated to an increment in strain and pore pressures. The model formulation enabled the use of a fully explicit model, in which automatic substepping procedures have been incorporated. The double effective stress increment, as well as the elastic and plastic strains partition, are deduced from the plasticity algorithm. It then serves as an input to the damage algorithm, to determine the damage increment. This algorithm has the advantage of been flexible enough to accommodate modifications of the different model components.

This constitutive and numerical framework has been validated against laboratory mechanical tests for two materials, the Callovo-Oxfordian argillite and Boom clay. Provided that the plasticity model is adapted to fit the specific material behaviour, a good agreement is reached between experimental data and simulations.The effect of confining pressure on the mechanical behaviour is well reproduced, and the ductile/brittle transition associated to suction increase is captured.

In order to be able to simulate fully coupled hydro-mechanical problems, transfer laws depending on the saturation state have been introduced. The balance and mass conservation equations, their weak formulation and the spatial and time discretisations, achieved by means of the Finite Element Method, have been presented. The global resolution method is based on the iterative Newton-Raphson algorithm. This complete numerical framework allowed us to demonstrate the capabilities of the constitutive framework, by simulating fully-coupled hydro-mechanical problems. The first illustration concerns damage induced by desiccation and humidification. Since the solid skeleton always remains in compression when submitted to hydric loading, effective stress based formulations do not allow for tensile damage to be captured. However, important deviatoric stresses result from high moisture gradients experienced in low permeable materials, which may induce damage. What is interesting in this approach is its ability to capture, not only desiccation-induced damage, but also damage due to humidification, which has been observed but not yet explained. The second application deals with the creation of the excavation damaged zone around deep galleries. A first illustrative example allowed us to study the influence of suction on damage and plasticity zones extent, considering drained conditions. The second example presents the pore pressure field perturbation, and the damage and plasticity zones due to an undrained excavation.

Perspectives

Many improvements could be incorporated into the constitutive model. First, the definition of the double effective stress could benefit from a micromechanical analysis, to determine more accurately the stress resulting from damage and suction on the clay matrix. Since micro-cracks usually appear along preferential directions, depending on loading, an extension to an anisotropic damage tensor would certainly improve the model response. Evolution of the fluid related laws, such as water retention properties or permeabilities, with damage should also be accounted for to give accurate predictions of hydraulic phenomena. In the case of wetting-drying cycles, hysteresis of the water retention curve should not be neglected either. Concerning the plasticity framework, lots of advanced forms of the Cam-clay model are available in the literature, which could be used, provided that the resulting improvements are worth the increased number of required parameters.

Concerning the Finite Element implementation, Θ-Stock is already appropriate to simulate many hydro-mechanical non-linear problems. However, its efficiency could be greatly improved by using Newton-Raphson with the tangent rigidity matrix, instead of the elastic rigidity matrix, which would increase the convergence rate. A valuable improvement would also be to incorporate methods to prevent the appearance of spatial oscillations in the pore pressure field in coupled hydro-mechanical problems. This would allow the choice of small time steps without the need of mesh refinement.

Finally, only a few applications of the proposed model have been presented here. Further validation can be made by simulating full-scale real problems. The excavation and ventilation of deep galleries can already be simulated. However, for more realistic results, we think that the dependence of permeabilities on damage should be accounted for. What would be very interesting is to extend our modelling approach to long term effects, such as viscoplasticity, self-sealing and self-healing. Indeed, although the rock is damaged by excavation, clay-based geomaterials tend to recover part of their original properties after a certain time, due to chemical reactions with water.

APPENDIX A

DERIVATIVES FOR CAM-CLAY PLASTICITY

This appendix presents the derivatives of the plasticity functions used for Boom Clay modelling (based on CC and BBM). They are necessary to calculate the stiffness matrix in the Newton-Raphson method (chapter 3.4), and to calculate strains and stresses increments in the local stress-point algorithm (see chapter 5.6.3).

Introduction1M

  Definition of the double effective stress 1.1 Unsaturated soils . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Modelling approaches in the literature . . . . . . . . . . . . . . . . . 1.1.2 Constitutive stress based on energy considerations . . . . . . . . . . 1.1.3 Water retention properties and their impact on constitutive stress . 1.2 Principle of effective stress in Continuum Damage Mechanics . . . . . . . . 1.3 Introduction of a double effective stress . . . . . . . . . . . . . . . . . . . . 2 Elastic behaviour 2.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Pressure dependence of moduli in Cam-Clay . . . . . . . . . . . . . . . . . . 2.3 Pressure dependent hyperelastic models . . . . . . . . . . . . . . . . . . . . 2.4 Effect of damage on elastic behaviour . . . . . . . . . . . . . . . . . . . . . 2.4.1 Principle of strain equivalence . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Principle of equivalent elastic energy . . . . . . . . . . . . . . . . . . 2.4.3 Comparison of both assumptions . . . . . . . . . . . . . . . . . . . . 2.5 Unsaturated elastic behaviour . . . . . . . . . . . . . . . . . . . . . . . . . . vii LIST OF SYMBOLS This list of symbols aims at listing the main symbols which are used throughout the thesis. Slope of the critical state line m vg Parameter of the Van Genuchten equation n Pressure exponent in HAR hyperelastic model n vg Parameter of the Van Genuchten equation

Contents 1 . 1

 11 Unsaturated soils . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.1.1 Modelling approaches in the literature . . . . . . . . . . . . . . . 8 1.1.2 Constitutive stress based on energy considerations . . . . . . . . 11 1.1.3 Water retention properties and their impact on constitutive stress 13 1.2 Principle of effective stress in Continuum Damage Mechanics 18 1.3 Introduction of a double effective stress . . . . . . . . . . . . . . 19

Fig. 1 . 2 :

 12 Fig. 1.2: Repartition of the different phases

Fig

  Fig. 1.4: Influence of suction on volumetric compression and volume changes due to wetting and drying

Fig. 1 . 5 :

 15 Fig. 1.5: Evolution of χ parameter with the degree of saturation[START_REF] Jennings | Limitations to the use of effective stresses in partly saturated soils[END_REF]

Fig. 1. 6 :

 6 Fig. 1.6: Hysteresis in the movement of water meniscus during drying-wetting cycles (from[START_REF] Miller | Physical Theory for Capillary Flow Phenomena[END_REF]

  Fig. 1.7: Typical Desorption and absorption curves

  are given for different porosities in figure 1.8. The reference curve gives φ ref = 0.39, α vg ref = 1 3.44.10 6 = 2.91.10 -7 , n vg = 1.49 and m vg = 0.33. Porosity dependent water retention curves computed using equation 1.32 are represented in figure 1.9.

Fig

  Fig.1.8:[START_REF] Delahaye | Soil heterogeneity and preferential paths for gas migration[END_REF]

Fig. 1 .

 1 Fig. 1.11: (eq. 1.32)

Fig 2

 2 Fig. 1.13: Mine tailings

  Fig. 1.16: Principle of effective stress

Fig. 2 . 2 :

 22 Fig. 2.2: Deviatoric and volumetric strains obtained during an elastic cycle with pressure dependent hyperelasticity

Fig. 2 Fig. 2 .

 22 Fig. 2.3: Principle of strain equivalence

  Fig. 2.10: p = 1e5 Pa

  Fig. 2.13: Water Retention Curve

Fig. 2 .

 2 Fig. 2.14: swelling σ v = 0.1 MPa

  Le critère d'endommagement et sa loi d'évolution ainsi que le cadre élasto-plastique sont introduits. Les équations d'endommagement et de plasticité sont exprimées en fonction de la contrainte doublement effective. Nous montrons comment les critères d'endommagement et de plasticité évoluent avec la succion et l'endommagement dans l'espace des contraintes doublement effectives ainsi que dans l'espace des contraintes totales. Contents 3.1 Thermodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.1.1 Second law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.1.2 Helmholtz free energy . . . . . . . . . . . . . . . . . . . . . . . . 39 3.1.3 Clausius-Duhem Inequality -dissipation . . . . . . . . . . . . . . 39 3.1.4 Hypotheses on the form of Helmholtz free energy . . . . . . . . . 40 3.2 Damage criterion and evolution . . . . . . . . . . . . . . . . . . . 41 3.2.1 General framework . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2.2 Specific functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.3 Plasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3.1 Definition of the preconsolidation pressure and hardening law . . 48 3.3.2 Dependence of the preconsolidation pressure on suction . . . . . 50 3.3.3 Yield surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.4 Plastic flow rule . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.4 Derivation of the incremental law . . . . . . . . . . . . . . . . . . 55

Fig. 3 . 5 :

 35 Fig. 3.5: Shape of damage criterion, Drucker-Prager -Effective stress space

Fig. 3 . 6 :Fig

 36 Fig. 3.6: Effect of damage on damage criterion -Total stress space

Fig

  Fig. 3.8: Evolution of damage with deviatoric stress

  Fig. 3.10: Volumetric behaviour in Cam-Clay

  Fig. 3.11: Volumetric behaviour -mean pressure exponent n = 1

  Fig. 3.12: Shape of the yield surface in p -s plane

Fig

  Fig. 3.13: Effect of suction on yield point and consolidation curve for n = 1

Fig. 3 .

 3 Fig. 3.14: Original Cam-Clay yield surface

  Fig. 3.15: Modified Cam-Clay yield surface

Fig

  Fig. 3.16: Effect of suction on yield surface -Total stress space

  Fig. 3.18: Plastic potential

  states that the stress state has to remain on the yield surface. the dependence of the preconsolidation pressure on suction and on the hardening parameter, p 0 , (equation 3.61)

σ

  = σ * + p g I -S l (s)sI (3.100) σ = σ * + ṗg I -S l (s) ṡ + ∂S l ∂s ṡs I (3.101)

Contents 4 . 1

 41 Water phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.1.1 Liquid water flow . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.1.2 Water vapour flow . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.1.3 Water mass conservation equation . . . . . . . . . . . . . . . . . 64 4.2 Air phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

Incorporating equations 4. 5

 5 and 4.23 into equation 4.3, the total water flux becomes Q m = ρ w (D P ∇s -K w ∇z) (4.25) in which D P = D P l + D Pvap (4.26) Deriving equation 4.2, the variation of the average water density can be expressed as ρm = (ρ w S l + ρ vap (1 -S l )) φ + φ(ρ wρ vap ) Ṡl + φS l ρw + φ(1 -S l ) ρvap (4.27)

  31)Incorporating equations 4.28, 4.29, 4.30 and 4.31 into equation 4.27, the increment of water density becomes: ρm = -(ρ w S l + ρ vap (1 -S l ))I : ε + φ(ρ wρ vap )g 2 ṡ + φS l β p ṗl + φ(1 -S l )A ṡ (4.32) Incorporation of equations 4.25 and 4.32 into equation 4.1 gives the mass conservation equation:

  Fig. 5.1: Two-dimensional analysis types (image from Zienkiewicz et al. (2005))

Fig. 5 . 2 :

 52 Fig. 5.2: Correspondence between the original and the reference configuration for a 4-node isoparametric element

Fig. 5 .

 5 Fig. 5.3: Approximation in the time domain

  Fig. 5.4: Numerical solving of the non-linear load-displacement equation

2 .

 2 Fig. 5.5: Explicit algorithm for plasticity

Fig. 6 .

 6 Fig. 6.1: Example of a stress-strain curve from a triaxial compression test with unloading-reloading cycles (from Chiarelli et al. (2003))

  Fig. 6.5: Uniaxial compression test

Fig. 6 .

 6 Fig.6.9: Stress-strain curves for different confining pressures

Fig. 6 .

 6 Fig. 6.10: Effective stress paths

Fig. 6 .

 6 Fig. 6.13: Proportional compression tests (a. k = 2, b. k = 5, c. k = 10)

Fig. 6 .

 6 Fig. 6.14: Lateral extension tests (a. p c = 30MPa, b. p c = 60MPa)

Fig. 6 .Fig. 6 .

 66 Fig. 6.18: Stress-strain curves for different values of C 1

Fig. 6 .

 6 Fig. 6.20: Bulk modulus degradation for different values of C 1

Fig. 6 .

 6 Fig. 6.22: Isotropic compression test (test TBOOM1 from Baldi et al. (1991))

Fig. 6 .

 6 Fig. 6.23: Oedometer test (from Horseman et al. (1987))

Fig. 6 .

 6 Fig. 6.24: Oedometer test (from Della Vecchia et al. (2011))

Fig. 6 .Fig. 6 .

 66 Fig. 6.25: Confinement 2 MPa (Test TBOOM2 from Baldi et al. (1991))

Fig. 6 .

 6 Fig. 6.28: Stress-strain curves for different confining pressures

  Fig. 6.29: Effective stress paths

Fig. 6 .

 6 Fig. 6.32: Stress-strain behaviour at different moisture contents (from Al-Shayea (2001))

Fig. 6 .

 6 Fig. 6.33: Loading-Collapse curve

Contents 7 . 1

 71 Damage induced by hydric loading . . . . . . . . . . . . . . . . . 112 7.1.1 Material parameters . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.1.2 Geometry and boundary conditions . . . . . . . . . . . . . . . . . 114 7.1.3 Desiccation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 7.1.4 Humidification . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 7.2 Excavation modelling . . . . . . . . . . . . . . . . . . . . . . . . . 120 7.2.1 Impact of suction on the extent and properties of the excavation damaged zone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 7.2.2 Simulation of excavation under undrained conditions . . . . . . . 125

  Fig. 7.1: Processes leading to the creation of the Excavation Damaged Zone

2 :

 2 Fig. 7.2: Microcracking due to humidification from 65 to 90 % RH (from Wang (2012))

7. 1 . 2

 12 Fig. 7.4: Geometry of the problem

5 :

 5 Fig. 7.5: Water pressure field evolution with time

6 :

 6 Fig. 7.6: Spatial distribution of total stresses along line AB at t=1h30

  Fig. 7.9: Stress paths at point B during desiccation

Figure 7 .

 7 Figure7.10 shows the water pressure field at four time steps. As for the desiccation test, the rate at which suction is applied results in a very heterogeneous water pressure distribution

Fig. 7 .

 7 Fig. 7.11: Damage pattern after humidification

  Fig. 7.12: Time evolution of double effective mean and shear stresses and damage at point B

  Fig. 7.13: Stress paths at point B during humidification

  Fig. 7.14: Geometry of the problem.

Fig

  Fig. 7.15: Total and double effective stress path for saturated state

  and 7.18 below show the evolution of these stresses around the tunnel, for two states of suction: s=0 and s=1MPa (R denotes the radius of the tunnel).

  Fig. 7.17: Total and double effective stress for saturated state

FigFig. 7 .

 7 Fig. 7.19: Plastic strain distribution at various constant suctions.

  Fig. 7.21: Geometry of the problem.

  Fig. 7.22: Spatial distribution of total stresses and pore liquid pressure at the end of excavation

  Fig. 7.23: Spatial distribution of radial displacementsThe spatial distribution of damage and plasticity (represented by the norm of the plastic strain matrix) is given in figure7.24.

  Fig. 7.24: Plastic strain and damage distribution at the end of excavation

  Fig. 7.25: Plastic strain, damage and liquid pressure time evolution at the tunnel wall

Fig. 7 .

 7 Fig. 7.26: Total and effective stress paths at the tunnel wall

  κ (λ(s)κ) 2 λ ′ (s) + S l (s) + ∂S l ∂s s with λ ′ (s) = -λ 0 [(1r)β exp(-βs)] (α vg s) nvg ) mvg+1 = g 2

Table 1 .

 1

	Lakeland Sand	0.12	0.35	3.8 0.41
	Mine tailings	0.10	0.018	3.4 18.6
	Yolo light clay	0.39	0.35	2.41 0.24
	Beit Netofa clay	0	3.6.10 -3	0.75 0.39

-1 ) n vg m vg 2: Van Genuchten parameters for different materials

[START_REF] Leong | Review of Soil-Water Characteristic Curve Equations[END_REF]

Table 2 . 2

 22 on Boom clay. The vertical stress is kept constant (σ v = 0.1 MPa) while suction is decreased from 230 to 0 MPa. Mechanical and retention properties are given intable 2.2.

	Elasticity				Retention
	n	p r	κ	ν	S r	α vg	n vg	m vg
		kPa				kPa -1
	0.6	10	0.007	0.35	0	0.17.10 -3	2	0.4

: Elasticity and retention parameters Retention properties are chosen to fit the experimental retention curve represented in figure 2.13.

  1111 D e 1122 D e 1112 /2 D e 1133 D e 2211 D e 2222 D e 2212 /2 D e 2233 D e 1211 D e 1222 D e 1212 /2 D e 1233 D e 3311 D e 3322 D e 3312 /2 D e 3333 D e 1111 D e 1122 D e 1112 /2 D e 1133 D e 2211 D e 2222 D e 2212 /2 D e 2233 D e 1211 D e 1222 D e 1212 /2 D e 1233 D e 3311 D e 3322 D e 3312 /2 D e 3333

						
						
						 	(5.40)
						
	such as	  	σxx	  		
		   	σyy σxy	   	=	   
		 	σzz	 

  C wu C ww + θ∆tK ww C wa + θ∆tK wa C au C aw + θ∆tK aw C aa + θ∆tK aa

				92)
	in which
	K W A =			(5.93)
	6. Convergency is reached if the two following criteria are respected:
	||Ψ n || ||Ψ 0 ||	< T OL	(5.94)
	||∆dof n || ||dof n-1 ||	< T OL	(5.95)
	7. Steps 2 to 6 are repeated until the two convergence criteria are respected.

Table 6 .

 6

	1	MPa	MPa	MPa

1: Material mechanical parameters

Hydraulic properties

Van Genuchten parameters are chosen to fit the experimental water retention curve given by

[START_REF] Hoxha | Saturated and unsaturated behaviour modelling of Meuse-Haute/Marne argillite[END_REF]

(see figure

6

.4). They are summarised in table

6.2.

Table 6 .

 6 Fig. 6.21: Water retention curve and S l s component of effective stress for Boom clay

	5: Retention parameters

  They are summarised in table 7.1. Dissipation is supposed to be due to damage only and plasticity is ignored.Fluid transfers are taken into account in this study. The permeability functions chosen as well as their parameters values are given in table 7.3. The initial porosity is taken equal to 15 %.

	Elasticity	Damage
	K	ν	C 0	C 1	C 2
	GPa		MPa	MPa
	5	0.19	0	200	0.2
	Table 7.1: Material mechanical parameters
		α vg (MPa -1 ) n vg m vg S r
		0.04	1.5 0.55 0
	Table 7.2: Retention parameters

Table 7

 7

.3: Hydraulic functions and parameters

  The goal of the numerical work presented below is not to simulate a real experiment but to illustrate qualitatively the capabilities of the model.

	Elasticity	Linear elasticity	ψ e 0 =	K 2	(ε e v ) 2 + Gε e d : ε e d
	Retention properties	Van Genuchten (1980)	S l (s) =	1 + (αs) n 1	m

Table 7 .

 7 4: Specific functions for claysBoom clay parameters Retention and mechanical parameters are given in Table7.5 and 7.6.

					α vg		n vg	m vg
					0.28	2.3	0.21
				Table 7.5: Retention curve parameters
	Elasticity	Plasticity						Damage
	K	G	M λ 0	κ	p r	r	β	ζ	C 0	C 1	C 2
	MPa	MPa				MPa		MPa -1		MPa	MPa
	250	115 1	0.16	0.017	5	0.74	1	0.4	0	10	0.8
				Table 7.6: Material mechanical parameters

  Table 7.7: Material mechanical parameters α vg (MPa -1 ) n vg m vg S r

					Plasticity			LC curve		Damage
	n	p r	κ	ν	M	λ 0	ζ		p r	β	r	C 0	C 1	C 2
		MPa							MPa MPa -1		MPa	MPa
	0.6 0.1	0.0075 0.15 0.8 0.00755 1.5	0.1	0.5	0.998 0	25	0.3
						0.28		2.3 0.21 0
					Table 7.8: Retention parameters
	Intrinsic permeability K int (φ) = K w0	φ 3 (1 -φ) 2

Remerciements

ABSTRACT

CHAPTER 5 NUMERICAL IMPLEMENTATION

This chapter presents the discretisation in space and time of the weak formulation of the initial boundary value hydromechanical problem. The global iterative resolution algorithm used in θ-Stock as well as the specific local algorithm developed for the damage-plasticity model are also described.

Ce chapitre présente la discrétisation spatiale et temporelle de la formulation faible du problème hydromécanique aux conditions limites. Le schéma de résolution global itératif utilisé dans θ-Stock, ainsi que l'algorithme local développé pour le modèle d'endommagementplasticité sont également décrits.

A generalised displacement vector is defined:

The strain vector is related to displacements through the strain operator B.

ε(x, t) = B u(t) (5.45)