a) Boundary flow condition (b) Idealisation of flow at the interface

This research was funded by: FRIA -F.R.S. -FNRS Fonds pour la formation à la Recherche dans l'Industrie et dans l'Agriculture Rue d'Egmont, 5 B-1000, Bruxelles iii

Remerciement

Cette thèse est le résultat d'un travail de très longue haleine. En préambule, je veux adresser tous mes remerciements aux personnes avec lesquelles j'ai pu échanger et qui m'ont aidée pour l'aboutissement de cette thèse.

En premier lieu, je tiens à remercier mes directeurs de thèse, Panagiotis Kotronis et Frédéric Collin, pour leur confiance qu'ils m'ont accordée en acceptant d'encadrer ce travail doctoral, pour leurs multiples conseils et pour toutes les heures qu'ils ont consacrées à diriger cette recherche. J'aimerais également leurs dire à quel point j'ai apprécié leur grande disponibilité et leur respect sans faille des délais serrés de relecture des documents que je leur ai adressés. Enfin, j'ai été extrêmement sensible à leurs qualités humaines d'écoute et de compréhension tout au long de ce travail doctoral. J'exprime aussi toute ma gratitude envers Robert Charlier pour m'avoir laisser rejoindre son équipe ainsi que pour le suivi de mon travail.

Je remercie également Zhenyu YIN, pour son encadrement, ses conseils, sa disponibilité et son aide précieuse lors de mon séjour à Nantes. Mes remerciements vont aussi à Benjamin Cerfontaine, qui n'a jamais hésiter à m'aider, et répondre à mes questions techniques et théoriques.

Je remercie aussi Le Fonds de la recherche Scientifique FRIA-F.R.S.-FNRS ainsi que l'École Centrale de Nantes qui ont financé conjointement mes travaux de recherche et m'ont permis ainsi de poursuivre mes études doctorales dans des meilleures conditions.

J'aimerais exprimer mes remerciements à tous les chercheurs et doctorants, trop nombreux pour citer, pour les moments qu'on a partagés ensemble lors des pauses café, des sorties et des voyages. J'ai une pensée particulière à mes collègues de bureau présent et passés : Chaofa, Roxana, Borana et Shilton, avec qui je garde mes meilleurs souvenirs à Nantes, Georgia qui m'a beaucoup aidée dès mon arrivée à Liège, et Gilles qui s'est marqué d'une grande sympathie et d'une grande aide.

Mes sincères remerciements vont aussi à ma famille : mon oncle Mustapha ainsi que sa petite famille. Ils ont été pour moi une solide référence durant mon séjour en France et en Belgique. Je remercie pareillement ma grande famille, mes oncles, mes tantes ainsi que mes cousins, sans oublier mon frère et ma soeur.

Cette thèse n'aurait pu voir le jour sans l'aide de mes parents, sans qui je ne serais pas là, je leurs remercie pour leurs encouragements, leurs prières et leurs conseils, je ne pourrais jamais les en remercier assez. Merci Papa! Merci Mama! Enfin, je remercie mon cher Mohamed pour son soutien quotidien indéfectible et son enthousiasme contagieux à l'égard de mes travaux comme de la v vi vie en général. Notre couple a grandi en même temps que mon projet scientifique, le premier servant de socle solide à l'épanouissement du second.

Abstract

Nowadays, offshore wind energy industry is developing exponentially, due to the significant contribution of the North Sea wind turbines energy production to the total consumed energy in Europe. Given that the EU's energy target is to increase the share of renewable energy by 2020, there is a great potential of the offshore wind energy applications towards this direction. In order to achieve this, the support of scientific research is crucial.

Monopiles have been by far the most common support structure for offshore wind turbine, recently becoming applicable also for complex site conditions. Monopiles foundations consist of hollow, thin-walled piles, and they are generally subjected to lateral loading.

The standard design methods are not always suitable for the optimized design of this type of foundation. Especially design of the large monopile under cyclic loading is still demanding the support of scientific research. Therefore, the main objective of this research work is to study the soil-pile interaction problem of offshore wind monopiles under both monotonic and cyclic lateral loading.

The work focuses particularly on monopiles embedded in clayey soils. Special attention is given to the behaviour of natural clays, by proposing a novel constitutive law for natural clay under monotonic and cyclic loading. Furthermore, the constitutive law is employed for numerical analyses of the laterally loaded monopiles, by the mean of finite element modelling. vii

Résumé

La contribution des parcs éoliens offshore en termes d'énergie renouvelable ne cesse de croître. L'électricité produite notamment par les éoliennes en mer du Nord représente de ce fait une part importante de l'énergie consommée en Europe. Pourtant, les objectifs de la Commission Européenne en termes de transition énergétique à l'horizon de 2020 sont loin d'être réalisés sans l'optimisation des techniques existantes et le développement de technologies de plus en plus innovantes.

Les monopieux sont les fondations les plus populaires dans l'industrie des éoliennes offshores. Ce sont des pieux creux, à parois mince. Ils sont généralement soumis à différentes sollicitations latérales de nature monotone et cyclique.

Les méthodes standard de dimensionnement des fondations ne sont pas toujours adaptées pour le calcul en service de ce type de fondation. Cela est particulièrement vrai pour les monopieux larges, sous chargement cyclique. L'objectif général de ce projet est d'étudier les interactions sol-structure des monopieux du secteur éolien offshore sous sollicitations latérales monotone et cycliques.

Notre travail se focalise particulièrement sur les monopieux fondés dans les sols argileux. Une attention particulière est donnée au comportement de l'argile naturelle. En effet, une nouvelle loi de comportement pour l'argile naturelle est proposée afin de prédire différentes propriétés des sols argileux sous chargement monotone et cyclique. Ensuite, la loi de comportement est utilisée pour la modélisation en éléments finis du problème du monopieux sous chargement latéral. Given the huge potential to produce energy from renewable resources across the continent, such as solar, wind, geothermal, hydroelectric, and other forms, the scope of this research is set up, which is concerned with energy produced from wind, specifically from offshore wind turbines.

In fact, offshore wind turbines has the potential to contribute with 3.6% to 4.3% of the total energy production in Europe, with a combined installation capacity of 40 GW by 2020, and even more in the coming years. More specifically, Belgium is privileged of having seven high quality areas in the North Sea devoted to establish wind farms, which contributes to around 10% of the total energy consumption. 30km off the northern coast lies the C-Power project, the first wind farm construction in Belgium that consists of more than 50 wind turbines till now. When it comes to feasibility study of wind turbines, the expenditure on the foundations reaches up to 30% of the overall cost of wind farms (see Figure 1.1).

Briefly, there are several factors that affect the design of offshore foundations, such as the distance from the shore, geotechnical conditions of the seabed and lastly the sea depth, which therefore influence the overall cost of the wind turbine. Thus, it has been of a great interest to compromise between foundation's technical and economic aspects in order to increase the economic viability of the wind farm projects, and therefore encourage investing in wind turbines, creating a competition between private companies achieving the most cost-efficient wind turbines, and therefore increasing their turnovers.

Wind turbines foundations and Monopiles

An offshore wind turbine is a light structure that consists of a turbine installed at the top of a support structure resting on the foundation. The depth of the seabed and the geotechnical conditions of the site are the main criteria for choosing the type of foundation. However, the monopile type of foundation remain the most widespread solution in Europe since they occupy 80% of the market for offshore wind turbine foundations (EWEA (2017)).

Monopile foundations are tubular steel piles with thick walls. They are subjected to horizontal forces and significant reversal moment under the action of the wind on the turbine and the marine current on the supper-structure.

The first monopile foundations used in wind turbines returns back to the 80's. However, their design were based on classical standards and guidelines, which are sometimes limited when complex loading are involved. Fortunately, computational power has increased for the last decades, which allows to perform detailed analyses. In particular, numerical modelling of offshore structures by the mean of finite element modelling are already common, but still considered as a challenging, related to the soil behaviour, especially when it comes with cyclic loading. Soils are non-homogeneous media, having a strongly non-linear behaviour which is impossible to inspect entirely. Moreover, producing the cyclic behaviour of soils often provides the need of sophisticated constitutive models.

Natural clays and constitutive modelling

Fortunately, geotechnical conditions in the North Sea are remarkably uniform . Most of the deposit encountered are very dense layers. Moreover, soil composition varies between sandy, clayey and gravel layers, as shown in Figure 1.2. This PhD thesis deals especifically with clayey soil. In particular, this research work aims at modelling the behaviour of natural clays under both monotonic and cyclic loading. The developed model will be employed for the finite element analysis of the OWT Monopiles, under both monotonic and cyclic lateral loading.

Finite element modelling of offshore Monopiles

Simulation of the behaviour of monopiles under lateral loading, by the mean of finite element modelling, requires the control of two main keys: the soil-pile interface and soil matrix. On one hand, the interface between soil and pile is the seat of complex mechanisms: shearing, strain localization, dilatancy, in addition to the eventual opening of gap, and the evolution of suction around the pile in addition to the water flow inside and through the interface. On the other hand, an accurate modelling of the soil behaviour is not an easy task, especially when it comes with natural and in-situ clays, which are almost anisotropic and characterized by high sensitivity, without forgetting the cyclic features of soils, that represent an additional complexity to the soil behaviour.

These two keys, namely soil and interface behaviours, will be deeply analysed along the thesis. This challenging task will be achieved by the mean of finite element code. LAGAMINE have been employed in several geotechnical problem, especially those characterized by hydro-mechanical and thermo-hydromechanical problems such as the framework of nuclear waste, the geothermal problems, and suction caisson for offshore geotechnics. The main objective of this research work is to study the soil-pile interaction problems of the offshore monopiles, under lateral solicitations. In particular, this work focuses on the numerical modelling of the monopiles embedded in clayey soil. More especially, the aim is to:

Model the behaviour of natural clays under monotonic and cyclic loading using elastoplastic model. The concept of bounding surface plasticity is introduced to the baseline model SCLAY-1S; Elaborate the soil reaction curves using a fully hydromechanical approach and finite element method, which allow studying the effect of loading rate and drainage conditions on the response of laterally loaded piles;

1. Analyse the effect of clay behaviour on the response of laterally loaded monopiles. More specifically, study the effect of soil anisotropy and clay structure on the p-y curves;

Study the behaviour of rigid and slender monopiles under cyclic lateral loading using 3D Finite Element Analysis (FEA).

Research outline

The manuscript consists of 8 chapters covering the literature review and the developments achieved during the PhD research. Those chapters are organized in a logical sequence, starting with the laboratory test till the numerical application.

Chapter 2 presents the general approach considered in the finite element modelling throughout the research work, in particular, the basis of the finite element formulation required to solve coupled hydro-mechanical problems. Special care is given to the formulation of the elastoplastic constitutive relation for clayey soils.

Chapter 3 presents an overview of the behaviour of natural clays under monotonic loading and its features, such as plastic anisotropy, destructuration and dilatancy, which are described and discussed from experimental studies. Furthermore, a state of the art of constitutive laws predicting these features is provided and analysed.

Chapter 4 presents the formulation of a novel constitutive law for natural clays under monotonic loading, named SCLAY-1SB model. The proposed model introduces the concept of bounding surface plasticity to the constitutive law of SCLAY-1S. The model is formulated under both axisymmetrical and multi-stress conditions. The model parameters are described and have been calibrated on different clay using different laboratory tests.

Chapter 5 concludes the first part of this work by dissecting the offshore application. The response of a laterally loaded monopile is deeply studied. Monotonic loading is investigated through 2D and 3D finite element models. The 2D FEA aims to analyse the effect of clay behaviour on the response of laterally loaded piles. Furthermore, this chapter presents a 3D modelling of the centrifuge test on a laterally loaded piles using the SCLAY-1SB model for the soil matrix. Theis chapter ends up by presenting a 3D FEA of a rigid monopile.

Chapter 6 deals with the intrinsic cyclic features of natural clays, aiming to extend the proposed model to the cyclic loading. Three key aspects are discussed and introduced in the model. First, a moving projection center is proposed with an adequate evolution rule. Second, the strain accumulation effects due to the cyclic loading are taken into account. Lastly, a non-associate flow rule is adopted to enhance the prediction of the undrained cyclic stress path. The calibration method of the new added model parameters are explained and applied to undrained triaxial test of Georgia Kaolin clay. Finally, the model is validated on in-situ tests of cyclic self-boring pressuremeter carried out on intact and in-situ clay.

Chapter 7 concludes this work by presenting two different finite elements analysis, using the final version of the proposed model. First, a 2D FEA of a pile section under cyclic lateral loading is presented and compared to the response of the lateral expansion cavity during the cyclic pressuremeter tests. Second, a 3D finite element applications are developed, dealing with the behaviour of rigid monopiles under cyclic lateral loading using in the both cases the final version of the SCALY-1SB model.

Introduction

The numerical modelling of offshore monopiles is a complex soil-structure interaction problem and requires therefore the definition of a rigorous framework. This chapter aims to present the general framework adopted in this research work and constitutes a comprehensive basis for the developments of the next chapters.

The basics of the finite element code of LAGAMINE are first presented. More particularly, the first two sections provide the main concepts of the hydromechanical finite elements employed in the numerical modelling, where both solid and contact elements are presented. The description gathers the principal kinematics assumptions and the equilibrium equations for both hydraulic and mechanical behaviours. In the third part, attention is paid to the constitutive modelling of geomaterials. The theory of elastoplasticity is described and the adopted method for integrating elastoplastic laws is presented. Finally, emphasis is placed on the modelling of clayey soils; the concept of the critical state is notably discussed and the family of Cam-Clay models is detailed.

This chapter constitutes a short summary and highlights several points useful for this PhD. The interested reader should however refer to [START_REF] Charlier | Approche unifiée de quelques problèmes non linéaires de mécanique des milieux continus par la méthode des éléments finis (grandes déformations des métaux et des sols[END_REF], [START_REF] Collin | Couplages thermo-hydro-mecaniques dans les sols et les roches tendres partiellement satures[END_REF], [START_REF] Barnichon | Finite element modelling in structural and petroleum geology[END_REF], and [START_REF] Cerfontaine | The cyclic behaviour of sand, from the Prevost model to offshore geotechnics[END_REF] for an extensive description of the LAGAMINE code.

2 Hydro-mechanics of porous media

Classic porous media

It is commonly assumed that geomaterials, such as clayey soils, exhibit a porous structure which is usually considered as an assembly of grains forming the solid skeleton and fluids (liquid and/or gas) filling the porous space between the grains. Therefore, an elementary volume dΩ of a granular porous medium can be defined as the juxtaposition of two constituents:

• solid skeleton, which is an arrangement of grains. The volume of the skeleton is denoted dΩ s and its specific mass ρ s . dΩ s is equal to the volume of the grains and usually considered incompressible. Consequently, the mechanical behaviour of the granular porous material is entirely due to the grains arrangement.

• void volume, which is filled with fluid(s) of a specific mass ρ f , in liquid and /or gas phases. When the void volume is fully filled by liquid phase, the soil is referred as saturated. Otherwise, the soil is partially saturated. The present work deals with the former case. The soil is therefore here-after assumed fully saturated by water in its liquid phase, which makes sense in the particular case of offshore engineering.

= +

Porous Medium Soil Skeleton Interstitial Fluid At the macroscopic scale, homogenisation of the porous medium allows to define an equivalent continuous porous medium, superposition of the fluid phase and the skeleton, by considering that at any geometrical point a solid grain and the fluid coincide. The open porosity n is defined by:

n = Ω v Ω = Ω v Ω s + Ω v ∈ [0, 1] (2.1)
where the subscript v refers to the void medium and s to the solid grains. Hence, the total specific mass of the porous medium ρ is expressed as:

ρ = (1 -n) .ρ s + n.ρ f (2.2)
where the subscript f refers to the fluid phase.

Kinematics of porous media

Kinematics of continuous solids

The infinitesimal description of the deformation remains valid as long as the final and the initial configurations are very similar. In this case, strains, usually referred as Cauchy strains, are defined as a linear function of the displacements.

However, when the final configuration is very different from the initial one, the infinitesimal description does no longer hold as the non-linearity of kinematics becomes important. A large strain description of the kinematics is therefore required, which must account for rigid body motion (large translation and/or rotation) and for large deformation. In this case, an appropriate referential frame must be chosen for the description of kinematics.

Generally, the non-linearity of kinematics is described either by the Eulerian or the Lagrangian formulation. In the former, the equilibrium is formulated with respect to the deformable configuration but the basis is considered 14CHAPTER 2. HYDROMECHANICAL MODELLING OF GEOMATERIALS fixed. In the latter, the equilibrium is formulated with respect to the reference configuration but the basis evolves with the deformed body. In other words, the Eulerian mesh is fixed in the space and material points move through it, while the Lagrangian mesh is linked to the material points and moves with it (Figure 2.2). Figure 2.2: Eulerian and Lagrangian configurations, from [START_REF] Cerfontaine | The cyclic behaviour of sand, from the Prevost model to offshore geotechnics[END_REF] In the finite element code LAGAMINE, equilibrium is expressed by an updated Lagrangian formulation in the current configuration. This means that the reference configuration is Lagrangian and evolves at every step (Figure 2.3). Let consider a deformable body B. At time t, ϕ t : B → R n is the mapping that associates a point X of the reference configuration B onto its current position x = ϕ t (X).

The velocity of material point x in the current configuration reads:

v = dx dt (2.3)
Consequently, the gradient of the velocity L in the current configuration is given by: L = ∂v ∂x (2.4) which can be decomposed into a symmetrical and an anti-symmetrical part:

L = 1 2 L + L T + 1 2 L -L T (2.5)
The Cauchy strain rate is given as the symmetrical part of the velocity gradient L:

ε = 1 2 L + L T (2.6)
and the spin rate is the anti-symmetrical one:

ω = 1 2 L -L T (2.7)
The large deformation/ displacement process needs an objective stress computation, which must be independent of the rigid body motion. Therefore, the Jaumann objective stress rate tensor σ is adopted in the following:

σ = σ -ω.σ -σ.ω T (2.8)
Finally, the sign convention of soil mechanics is used, i.e., compression stress and compaction strains are counted positive.

Kinematics of fluid phases

The kinematics of continuous solids could be applied to the solid skeleton of a porous medium, since the skeleton is considered continuous. Nevertheless, it cannot be adopted for the fluid phase. On the opposite of continuous solids, porous media are open systems, which means that any elementary volume dΩ might exchange fluid with its environment. Moreover, the mechanical behaviour of the solid skeleton might reduce or increase the void volume dΩ v and thus influence the fluid behaviour. Hence, a different kinematics configuration must be adopted for the fluid phases. The code of LAGAMINE adopts the Eulerian description of the fluid. According to [START_REF] Coussy | Mécanique des milieux poreux[END_REF], the Eulerian description of fluid flow defines a control volume, without identifying individual fluid particles. Pressure, velocity and all other flow properties are described as fields within the control volume.

Let consider an elementary surface da, whose normal is n (the current configuration in Figure 2.3). At time t and per unit of time, a fluid mass J m f flows though da such as:

J m f = f w .nda (2.9) in which the fluid mass vector f w is expressed as:

f w = ρ f .v d (2.10)
where v d is the relative velocity of the fluid with respect to the solid skeleton, known as Darcy's velocity. It is written:

v d = n. (v w -v s ) (2.11)
where v w is the average spatial velocity of the water phase, v s the spatial velocity of the solid phase and n refers to the soil porosity.

Effective stress

The stress-strain constitutive relations must distinguish an effective stress tensor, which represents only the stress acting on the solid skeleton and hence causing the soil deformations. For porous materials that are saturated with water (one fluid), Terzaghi's postulate is conventionally employed. Accordingly, the total stress tensor σ is expressed as:

σ = σ + u w δ (2.12)
where δ is the identity tensor. σ is the effective Cauchy stress tensor, defined under the soil mechanics convention, in which compression stress is positive and u w is the pore water pressure. It must be point out that the Equation 2.12 is only representative of saturated soils, where grains are considered incompressible.

Stress and strain representation

The Cauchy stress tensor is a symmetrical tensor composed of 6 components

σ =    σ 11 σ 12 σ 13 σ 22 σ 23 Sym. σ 33    (2.13)
Invariants of the stress tensor are usually employed to represent the stress state independently of the referential frame. In the field of geomechanics, the stress tensor is conventionally described by a triplet of stress invariants

[I σ , II s , III s ].
The first invariant is given by I σ = σ ii (2.14) which is the trace of the stress tensor. The deviatoric stress tensor is then defined as: 2.15) where δ refers to the identity tensor.

s = σ - I σ 3 .δ (
The second and the third invariants are calculated from the deviatoric stress tensor as:

       II s = 1 2 s ij .s ij III s = 1 3 s ik .s kl .s lj (2.16)
The use of this triplet of invariants allows a complete representation of the stress state. Indeed, the three invariants act as a new polar coordinate system, as illustrated in Figure 2.4. The first invariant allows defining the deviatoric plane or Π-plane, the second invariant provides the distance between the stress state and the hydrostatic axis and the third invariant the Lode angle β, which determines the direction with respect to the pure shear line, obtained as: Furthermore, the first and the second invariant of the stress tensor are usually reformulated into the mean and deviatoric stresses, respectively:

β = - 1 3 sin -1 3 √ 3 2 III s II 3 s (2.
   p = I σ 3 q = √ 3II s (2.18) 18CHAPTER 2.

HYDROMECHANICAL MODELLING OF GEOMATERIALS

Similarly, the Cauchy strain tensor is a symmetric tensor composed of 6 components ε ij composed of a volumetric and a deviatoric par:

ε = ε d + ε v δ (2.19)
where ε d is the deviatoric strain tensor and ε v is the volumetric strain, given as the trace of the strain tensor (ε v = ε ii ).

The norm of the deviatoric strain tensor, denoted ε q , is usually used, and is given by:

ε q = 3 2 ε dij .ε dij (2.20)

Equilibrium conditions

Mechanical equilibrium

A quasi-static analysis is considered throughout the present research work, either under monotonic or cyclic loading conditions. If the former case is evident, the later could be justified by the relative low frequency of the cycles [START_REF] Leblanc | Response of stiff piles to random two-way lateral loading[END_REF]). Thence, inertial forces are neglected, and the local equilibrium equation is given in this case by:

divσ + ρg = 0 (2.21)
where σ refers to the total Cauchy stress tensor, g is the vector of gravity acceleration, and ρ is the apparent specific mass of the soil, determined in function of porosity and the specific mass of solid grains and fluid (Equation 2.

2).

A weak form of this local equilibrium equation can be obtained by considering an admissible virtual velocity field δv (i.e., a velocity field that respects the solid continuity and the boundary conditions). The result is known as the principle of virtual power, and states that, for any admissible δv, the solid is in equilibrium if the internal virtual power is equal to the external one. It reads in the current configuration:

ϕ(B) σ : ε.δvdV δ ẆI = ϕ(B) ρ.g.δvdV + ϕ(∂B) q.δvdΓ δ ẆE (2.22)
where ε is the virtual strain tensor, and q is the vector of surface forces. Equation 2.22 requires to define the behaviour of the solid skeleton and the fluid flow. For the mechanical behaviour, a constitutive relation has to be adopted between the effective stress and strain rates in the general incremental form:

σ = f ( ε, σ , κ) (2.23)
where κ represents a set of internal variables. Further details on the mechanical constitutive modelling of geomaterials are developed in Section 4.

The equations related to the hydraulic part gather the constitutive law of fluid flows, and the mass conservation equation, or storage law, as explained hereinafter.

Hydraulic equilibrium

As mentioned previously, a saturated porous medium can be separated in a solid phase and a fluid phase. Therefore, equations related to the hydraulic problem are necessary to fully describe the boundary value problem. Hence, the constitutive law of the fluid flow and the mass conservation equation have to be stated in the strong and weak forms.

Flow law

Similarly to the constitutive law for the solid skeleton, the fluid mass flow can be related to the fluid pressure. In the case of laminar flow and isotropic permeability, the Darcy's law is conventionally employed to describe the flow behaviour in saturated porous materials. It states that the interaction force between the skeleton and the fluid is proportional to the relative velocity between the two phases:

v d = - k µ f . (∇u w + ρ f g∇z) , ∇z = [0 0 -1] T (2.24)
where k is the intrinsic isotropic permeability of the porous medium and µ f is the dynamic viscosity of the fluid (µ f = 10 -3 Pa.s for water at 20 • C).

Storage law

Let's consider an elementary volume dV 1 in its initial configuration at t 1 . The porous volume is characterized with porosity n 1 and is filled by a fluid of specific mass ρ 1 f . The fluid mass in the elementary volume dV 1 is equal then to ρ 1 f n 1 dV 1 . Similarly, at time t 2 , the elementary volume is assumed in a deformed configuration where it occupies a volume dV 2 , with a porosity n 2 and a fluid specific mass ρ 2 f . The fluid mass at the deformed configuration is equal to ρ 2 f n 2 dV 2 . It can also be expressed as ρ 1

f n 1 dV 1 + SdV 1 , where SdV 1 represents the fluid mass exchanged between t 1 and t 2 by the elementary volume dV 1 . Therefore, the change in fluid mass content S with respect to the initial configuration reads [START_REF] Biot | Variational lagrangian-thermodynamics of nonisothermal finite strain mechanics of porous solids and thermomolecular diffusion[END_REF]):

S = ρ 2 f n 2 dV 2 dV 1 -ρ 1 f n 1 (2.25) 20CHAPTER 2.

HYDROMECHANICAL MODELLING OF GEOMATERIALS

The volume ratio between the initial and deformed configuration represents the determinant of the Jacobian transformation, denoted J, i.e J = dV 2 /dV 1 .

The evolution of the specific mass of the fluid phase is considered to follow the linearised state equation given according to its bulk modulus K f . It reads:

ρ 2 f = ρ 1 f + ρ 1 f u 1 w -u 2 w K f (2.26)
Introducing the Equation 2.26 into Equation 2.25, in addition to some manipulations, the following storage law is obtained:

S = ρ 1 f (J -1) + n 1 u 2 w -u 1 w K f + ρ 1 f (J -1) (2.27)
Finally, in the case of infinitesimal displacement, J is related to the volumetric strain, i.e J = 1 + ε v and the porosity term is approximated, i.e. n 1 J -1 + n 1 . In this case, Equation 2.27 is simplified:

S = ρ 1 f n 1 u 2 w -u 1 w K f + ρ 1 f ε v (2.28)
The Equation 2.28 represents the hydro-mechanical coupling, since it stems from both the compressibility of the fluid and the deformation of the solid skeleton.

According to [START_REF] Barnichon | Finite element modelling in structural and petroleum geology[END_REF], the overall fluid mass conservation reads:

Ṡ + divf w = Q (2.29)
where the divergence operator is relative to Eulerian coordinates, f w is the flow of fluid mass vector and Q is the imposed volume flux.

Finally, the weak formulation can be found from the equilibrium Equation 2.29, by applying the principle of virtual power and considering an admissible virtual displacement of pore pressure δu w . The weak formulation of the problem in the current configuration is stated as follows:

ϕ(B) Ṡ.δu w -f w .∇ (δu w ) dV δ ẆI = ϕ(B) Q.δu w dV + ϕ(∂B f ) q.δu w dΓ δ ẆE (2.30)
where q is the imposed surface flux and ∂B f is the boundary where the flux q is applied.

Hydromechanical contact problem

The finite element modelling of offshore monopiles requires not only the definition of the hydromechanics of porous media (seen in section 2), but also the behaviour of the soil-structure (pile) interface. In the following, the principal concepts and constitutive relations adopted in the finite element code LAGAMINE to develop efficient contact finite elements are recalled. The presented concepts have been essentially summarized from the work of [START_REF] Barnichon | Finite element modelling in structural and petroleum geology[END_REF] and [START_REF] Cerfontaine | The cyclic behaviour of sand, from the Prevost model to offshore geotechnics[END_REF].

Stress and deformation vectors

Let's consider two deformable bodies B 1 and B 2 with boundaries ∂B 1 and ∂B 2 in contact through the boundary ∂B c (Figure 2.5).

Figure 2.5: Contact between deformable solids

At any point X of the contact surface, a local reference can be defined for each solid as indicated on Figure 2.5. The e 1 axis is normal to the contact, while e 2 and e 3 are tangent. In this local referential, the stress tensor in the contact surface of each body reduces to the contact stress vector t, which has three components:

t = p N [t T ] T T (2.31)
where p N is the normal pressure and t T = [t T 2 t T 3 ] the vector of shear stress, reduced to a scalar in the case of a 2D geometry. The contact stress vector t is associated with the deformation vector g of three components measuring the relative displacements of both solids:

g = g N [g T ] T T (2.32)
where g T = [g T 2 g T 3 ] the vector of relative tangential displacements and g N the gap function, which measures the distance between the two bodies with respect to the normal direction. The gap function is defined as a distance that relates every point X 2 of the contact boundary ∂B 2 to a point X 1 related on ∂B 1 , such that:

g N = X 2 -X 1 .e 1 1
(2.33)

The point X 1 in Equation 2.33 is the closest projection point of the so-called slave point X 2 onto the master surface ∂B 1 (see Figure 2.5).

It is noteworthy that the gap function is one of the principal features of the contact problem. The definition given above allows capturing the interpenetration of the two bodies, if g N < 0, as well as the gap between them, if g N > 0.

Ideal mechanical contact constraint

The ideal mechanical contact constraint prevents the penetration of body B 1 into B 2 . This condition is written as g N ≥ 0. If contact occurs the gap is null (g N = 0), while the associated normal component (p N ) of the effective Cauchy stress vector along the interface must be non-zero.

The ideal mechanical contact constraint is also named condition of Hertz-Signorini-Moreau and is stated as:

g N ≥ 0, p N ≥ 0, p N .g N = 0
(2.34)

If a gap exists between the two solids, i.e. g N > 0, the contact pressure is null. If contact happens, i.e. g N = 0, it involves a positive contact pressure.

If friction is allowed at the contact zones between the two bodies, two cases have to be distinguished: the stick and slip states. The stick state corresponds to the state when two points in contact are not allowed to move in the tangential direction. Hence, no relative tangential motion occurs during the displacements of the bodies. The slip state is the opposite. The transition from the stick to the slip state is ruled by a constitutive law, which defines a threshold of admissible shear stress:

f (t, κ) < 0 (2.35)
where κ is a set of internal variables. If f = 0, the slip condition stands and a plastic relative tangential displacement occurs, i.e. g T ≥ 0. Otherwise, both bodies are stuck and no relative motion occurs, i.e. g T = 0. A formalism identical to Equation 2.34 can be employed to write these conditions:

g T ≥ 0, f ≤ 0, g T .f = 0 (2.36)

Mechanical equilibrium

Both bodies, B 1 and B 2 , are subjected to equilibrium conditions, identical to those defined in Section 2. Furthermore, a new condition holds, where the contact occurs between the two bodies, through the boundary ∂B c . Therefore, a mechanical contact problem is simply the classic problem of solid mechanics with an additional constraint related to the contact boundary. Accordingly, each body B α , α = 1, 2 verifies the classic equilibrium equation stated in Section 2, i.e.:

divσ

+ ρ α g = 0 (2.37)
The solution of the problem consists in finding the field of velocity δv, for all x ∈ B α , that verifies Equation 2.37 and subjected to the constraint inequality g N ≥ 0 over ∂B α c . Henceforth, applying the principle of virtual power in the current configuration including the contribution of both solids is rewritten as: 2.38) where q α is the applied surface forces on the boundary ∂B α .

α=1,2 Bα σ.ε.δv T dV ≥ α=1,2 Bα ρ α g.δvdV + ∂Bα q α .δvdΓ (
The inequality in the weak formulation finds its origin in the contact constrain inequality (g N ≥ 0). If the contact zone has been previously established, the equation stands instead of the inequality:

α=1,2 B α σ.ε.δv T dV δ ẆI = α=1,2 B α ρ α g.δvdV + ∂B α q.δvdΓ δ ẆE +δ Ẇ c E (2.39)
where δ Ẇ c E is the mechanical contact contribution to external virtual work. In the finite element code LAGAMINE, this component is calculated by the penalty method, which enforces the contact constraint while authorising an interpenetration. The contribution to external virtual power reads:

δ Ẇ c E = ∂Bc [K N .g N .δ ġN + K T .g T .δ ġT ] dΓ (2.40)
with K N and K T the penalty parameters, which control the contact constraint.

The exact solution can be recovered for (K N , K T ) → ∞. However, very high penalty parameters might lead to ill-conditioned systems, numerical troubles and loss of accuracy of the results [START_REF] Cerfontaine | 3d zero-thickness coupled interface finite element: formulation and application[END_REF]).

Hydraulic characterization of the interface

Ideal flow conditions

In the LAGAMINE code, the void between both solids is modelled as an independent equivalent medium, B 3 , in which a fluid flow holds. Moreover, the flow

Hydraulic equilibrium equations

Similarly to the mechanical formulation, both solids verify the local hydraulic equilibrium (Equation 2.29), in addition to the effect of the flow at the contact surface. Therefore, the weak form of virtual power principle, including the fluxes in both bodies as well as in the interface is written in the current configuration as:

α=1,2 Bα Ṡ.δ uw -f w .∇ (δu w ) dV δ ẆI,f +δ Ẇ c I,f = α=1,2 Bα Q.δu w dV + ∂Bα q.δu w dΓ δ ẆE,f +δ Ẇ c E,f (2.41)
where δu w is an admissible field of pore water pressure rate, δ Ẇ c I,f , and

δ Ẇ c E,f
are the contribution of the interface to the internal and external virtual powers, respectively. The former stands for the fluid flow inside the gap B 3 , noted f l w :

δ Ẇ c I,f = B 3 Ṡ.δu w -f l w .∇ (δu w ) dV (2.42)
while the later represents the effect of the fluxes coming from the solid volumes B 1 and B 2 , namely f T 1 w , and f T 2 w :

δ Ẇ c E,f = B 3 f T 1 w .δu w dV + ∂B 3 f T 2 w .δu W dΓ (2.43)

Constitutive relations of contact elements

Mechanical law

The mechanical constitutive law for the interface element relates the effective stress rate to the relative velocity by means of the local compliance tensor. The stiffness parameters act as penalty coefficients and are considered rate independent:

ṫ = C st : ġ (2.44)
where the two components of the local effective stresses are:

• ṫT : the rate of shear stress.

• ṗ N : the rate of normal effective stress.

The two components of the relative velocity are:

• ġT : the relative tangential velocity.

• ġ N : the normal relative velocity.

The compliance matrix depends on the state of the interface (without contact, stick or slip). When there is no contact, the compliance tensor and the stresses are null. The stick and slip states are modelled within the framework of elastoplasticity. For the sake of simplicity, the Mohr-Coulomb criterion and perfect plasticity is assumed hereafter. The yield surface reads:

f ≡ |t T | -µ • p N -c ≤ 0 (2.45)
where µ and c are the friction coefficient and the cohesion of the interface, respectively.

In the general case, the total relative velocity is split into its elastic (stick) and plastic (slip) components:

ġ = ġst + ġsl (2.46)
Therefore the effective stress rate inside the interface is calculated as:

ṫ = C st . ġ -ġsl (2.47)
For the stick state, the compliance tensor becomes (using the penalty parameters K N and K T ):

C st = K N 0 0 K T (2.48) 26CHAPTER 2.

HYDROMECHANICAL MODELLING OF GEOMATERIALS

The slip component is defined as:

ġsl = λ. ∂h ∂t (2.49)
where λ is the plastic multiplier and h stands for the plastic potential. An associative rule is considered, i.e. h = f .

When the slip state is reached, the consistency condition must be enforced, i.e.:

ḟ = ∂f ∂t . ṫ + ∂f ∂µ μ + ∂f ∂c ċ = 0 (2.50)
where μ and ċ allow hardening of the yield surface. They are calculated as:

μ = λ.h µ , ċ = λ.h c (2.51)
As already mentioned, a perfect plastic behaviour is considered in this work, i.e. μ = 0, and ċ = 0.

Introducing Equations 2.49, 2.50 and 2.51 into Equation 2.47 leads to the explicit computation of the final stress state:

ṫ = C st . ġ -λ. ∂f ∂t (2.52)
Finally, the local elastoplastic compliance tensor is computed:

C st = 0 K N .µ. t T |t T | 0 K N (2.53)

Hydraulic laws

Similarly to the mechanical laws, hydraulic constitutive laws relate the fluid flux f w to the gradient of pressure ∇u w and the gap function variation ġ.

f w = h (∇u w , ġ) (2.54)
Three different fluxes are distinguished: (i) the storage flow f s w which comes from the mechanical opening/closing of the interface and / or the variation of the specific mass of the fluid filling the interface; (ii) the longitudinal flow f l w which represents the flow through a preferential path in the interface, and (iii) the transversal fluxes f t w,α standing from the exchange between sides of the interface and the surrounding bodies. i Storage f s w The main hypothesis is that the interior medium of the interface is fully saturated with water. Consequently, a storage flux Ṡ is necessary to verify the mass conservation equation.

The storage flux per unit of surface is obtained from the development of the mass conservation equation applied on an elementary volume V = g N .L (see Figrue 2.6), for which the mass conservation equation reads:

Ṡ = d dt (ρ w .V ) (2.55)
The specific mass of the fluid is assumed pressure dependent and follows a linear relation:

ρ w = ρ w,0 . 1 + u w -u w,0 χ w (2.56)
where ρ w,0 and u w,0 are the reference specific fluid mass and pressure, respectively and χ w is the fluid compressibility.

Hence, Equation 2.55 can be re-written as:

Ṡ = d dt (ρ w .g N .L) (2.57) = ∂ρ w ∂u w uw .g N .L + ρ w . ġN .L + ρ w .g N . L (2.58) = ρ 0 δu w χ w .g N .L + ρ w . ġN .L (2.59)
Finally, the storage flux is straightforward:

f s w = Ṡ L = ρ w,0 χ w . uw .g N + ρ w . ġN (2.60) ii Longitudinal f l w
The generalized Darcy law is adopted to describe a laminar flow between the two sides of the interface:

f l w = -k l µ w (∇u w + ρ w .g.∇s) .ρ w (2.61)
where k l [m 2 ] is the permeability in the longitudinal direction, g the gravity acceleration and ∇s is the gradient of the curvilinear coordinate in the longitudinal direction.

iii Transversal f t w,α

According to [START_REF] Cerfontaine | The cyclic behaviour of sand, from the Prevost model to offshore geotechnics[END_REF] the transversal fluxes are considered proportional to a transverse transmitivity T t w and they only depend on the pressures between solids in contact and the interior medium (B 3 in Figure 2.6)

f t w,1 = T t w,1 . u 1 w -u w .ρ w (2.62) f t w,2 = T t w,2 . u w -u 2 w .ρ w (2.63)
where u 1 w is the pore pressure on the solid B 1 , u w the pore pressure inside and u 2 w on the body B 2 . The transversal flow differs if a fracture is opened or closed, and filled with water or alterred material. Details for each case are provided in [START_REF] Cerfontaine | The cyclic behaviour of sand, from the Prevost model to offshore geotechnics[END_REF].

Constitutive modelling of geomaterials

The constitutive model is a complementary set of equations to the balance equations and kinematic conditions. Together with the loading and boundary conditions they allow a complete formulation of a boundary value problem and thus to calculate the motion and deformations of a given body. Constitutive laws are an important field of engineering mechanics and require clear premise and definitions. The major concepts and theories useful in this research work are briefly introduced hereafter.

Throughout this section, all the stress quantities are regarded as effective stresses and a compression positive convention is used for both stress and strain measures.

Modelling approaches

Elastic models

In elastic models the current stress is function of the current strain, the behaviour does not depend on the stress or strain path. The material can recover its original configuration, i.e, the deformation is reversible and there is no energy dissipation.

The elasticity theory dates back to the end of 17 th century, when the English scientist Robert Hook discovered that the displacement of some materials was proportional to the applied force. This discovery was the foundation of the linear elasticity in terms of force and displacement. In a multidimensional space, the elastic behaviour is often described by the generalization of Hook's law, which reads:

σ = E : ε (2.64)
where E is a fourth order elastic stiffness tensor. If E depends on the current stress state, the elasticity is non-linear, and termed hypoelasticity, whereas if the elastic stiffness tensor depends on the direction of loading, the elasticity is called anisotropic.

Elastoplastic models

The basic hypothesis is the decomposition of the total strain or its rate into elastic (recoverable) and plastic (irrecoverable) components. This is done in the literature in two ways;

• Using the additive decomposition of an appropriately defined finite strain tensor ε into elastic and plastic components [START_REF] Green | A general theory of an elastic-plastic continuum[END_REF])

ε = ε e + ε p (2.65)
• Using the multiplicative decomposition of the deformation gradient into elastic and plastic components [START_REF] Lubarda | Constitutive analysis of large elasto-plastic deformation based on the multiplicative decomposition of deformation gradient[END_REF])

= e p (2.66)
The first method is adopted by the majority of researchers in the field of geomechanics, while the multiplicative decomposition is motivated by the deformation process in crystal plasticity. Throughout this dissertation, the additional decomposition of total strain increment into elastic and plastic parts is considered within the framework of elastoplasticity. The main ingredients of an elastoplastic constitutive law are listed hereinafter:

• Elastic properties, the way in which the soil elastic deformations are defined.

• Yield criterion, is the boundary of elastic and plastic deformation.

The yield criterion is a function of the stress state and a set of internal variables, denoted κ, which varies with plastic strains to reproduce the magnitude of yielding.

• Flow rule, which determines the direction of plastic deformation. In the framework of elastoplasticity, it is assumed that plastic deformation derives from a plastic potential, which defines the mode of plastic deformation at yielding.

It is also often assumed that the plastic potential surface and the yield surface are identical or at least have the same direction. In this case, the material is said to obey the postulate of normality (regarding the stress state) or follows an associate flow rule (regarding the strains). Otherwise, the flow rule is called non-associated.

• Hardening rule describes the way in which the state parameter κ depends on stress and plastic strain states. There are two principal families of hardening rules, as shown in Figure 2.7: isotropic hardening, when the size of the yield surface varies with plastic strains, and kinematic hardening rule, when the yield surface moves with plastic deformations. In the general case, a combination of isotropic and kinematic hardening rules could be employed in the constitutive law.

Other approaches

Apart from the elastoplastic theory, which is the adopted modeling approach in the present research work, various theories of constitutive modelling exist and are commonly employed in geomechanics. Among them, the hypo-plasticity and hyper-plasticity are often adopted.

In Hyper-plasticity, originally suggested by [START_REF] Ziegler | An introduction to Thermomechanics[END_REF] and based on thermo-mechanical principles, the constitutive behaviour of a dissipative material is completely defined by two potential functions. The first function is either the Gibbs free energy or the Helmholtz free energy, whereas the second potential is the dissipation function. The interpretation of this framework in terms of conventional plasticity demonstrates that the classical yield and plastic potential surfaces, and hardening rules stem from these two scalar potential functions [START_REF] Ziegler | An introduction to Thermomechanics[END_REF]).

Following [START_REF] Lanier | A general formulation of hypoplasticity[END_REF], the Hypo-plasticity omits the assumption that the plastic strain is derived from a plastic potential, and does not a priori distinguish between elastic and plastic strains. Plastic deformation can therefore occur at the very beginning of loading process.

Numerical integration of elastoplastic laws

Finite element codes such as LAGAMINE employ an incremental formulation of the constitutive law. Henceforth, the incremental form of the general equations associated with an elastoplastic constitutive model is presented hereinafter.

Incremental formulation of elastoplasticity

The present research work deals with the rate-independent elastoplastic constitutive laws. The time is therefore fictitious and the material time derivative or rate, denoted hereinafter by a superposed dot, indicates simply the incremental form.

1. Decomposition of total strain increment :

ε = εe + εp (2.67) 2. Stress-strain relation σ = E : ε -εp (2.68)
where E is the fourth-order elastic stiffness tensor.

Flow rule

εp = λ • ∂g ∂σ (2.69)
where λ is the plastic multiplier. In the case of elastic loading, λ = 0, which leads to εp = 0

Yield surface:

The domain of the stress states where no plastic deformation occurs is termed elastic domain. The boundary of this region in the stress space is called yield surface, which is a convex surface in a general stress space.

Based on the yield surface definition, the elastic domain is defined in the stress space as follows:

E = f {(σ, κ) | f (σ, κ) ≤ 0} (2.70)
where f is the yield function.

In the strain space, the elastic domain can be defined as follows:

E = f (σ, κ) εp = 0 ⇔ λ = 0 (2.71)
Both Equations 2.70 and 2.71 can be gathered in one condition, called Karush-Kuhn-Tucker condition:

λ • f (σ, κ) = 0, where: λ ≥ 0; f (σ, κ) ≤ 0 (2.72) 32CHAPTER 2.

HYDROMECHANICAL MODELLING OF GEOMATERIALS

This equation provides the definition of the elastic and plastic domains.

In the former, f (σ, κ) < 0, and λ = 0, while in the later λ > 0 and f (σ, κ) = 0.

Hardening rule

κ = λ • h (σ, εp ) (2.73)
h is function of effective stress and plastic strain invariants.

Consistency condition

The consistency condition ensures that the stress state always lies on the yield surface, if plastic strains occur.

ḟ = ∂f ∂σ • σ + ∂f ∂κ • κ = 0 (2.74)
This set of equations represents the general framework of elastoplasticity for an arbitrary yield surface f , plastic potential g, and internal variables κ. Moreover, relationships between these components may exist. For example, combination of Equations 2.73 and 2.74 leads to the relationship between the plastic multiplier and the yield function, as fallows:

λ = 1 K p ∂f ∂σ • σ (2.75)
where K p is the plastic modulus:

K p = - ∂f ∂κ • h (σ, κ) (2.76)
On the other side, introducing Equations 2.67, 2.68 and, 2.69 yields:

λ = ∂f ∂σ • E : ε K p + ∂f ∂σ • E : ∂g ∂σ (2.77)
Finally, the expression of the elastoplastic tensor can be deduced:

σ = E -λ ∂f ∂σ • E • ∂g ∂σ T K p + ∂f ∂σ • E : ∂g ∂σ : ε (2.78)

Discretised formulation

Solution of Equation 2.78 over the time interval t ∈ [0, T ] requires discretisation of the time into N steps such that:

t 0 = 0, ∆t = T N , and
t n+1 = t n + ∆t (2.79)
For an initial value problem, the discretization of the constitutive equations is based on the Euler integration scheme, which considers that the variable X at the time step t n is known and assumes that its variation ∆X is linear between t n and t n+1 . Therefore, the constitutive Equation 2.78 is discretised for a given finite increment of deformation ∆ε corresponding to a time step ∆t:

     ∆σ = F (σ n+θ , κ n+θ ) ∆ε σ n+θ = (1 -θ) .σ n + θ.σ n+1 κ n+θ = (1 -θ) .κ n + θ.κ n+1 (2.80)
which is termed the generalised midpoint rule. The integration is termed explicit (forward Euler scheme) if all the variables required to compute F are known at the beginning of the step (θ = 0). On the other hand, if variables required to compute F have to be computed in n + 1 (θ = 1) the integration is called implicit (backward Euler scheme). If θ = 1 /2 the integration scheme is termed midpoint rule.

Application of the midpoint rule to an arbitrary elastoplastic model that is represented by the Equations 2.68 and 2.73, leads to:

∆σ = E : ∆ε -E : ∂g ∂σ n+θ ∆λ (2.81) ∆κ = ∆λ • h (σ n+θ , κ n+θ ) (2.82)
where ∆λ is the discrete plastic multiplier.

The yield criterion has to be verified at the end of the step:

f (σ n+1 , κ n+1 ) ≤ 0 (2.83)
where f (σ n+1 , κ n+1 ) denotes the yield criterion at step n + 1. Thence, if the stress state is elastic, no plastic deformation occurs, i.e., ∆λ = 0. Otherwise, the evolution of the plastic strains is established such that:

ε p n+1 = ε p n + ∆λ • ∂g ∂σ n+θ (2.84)
Therefore, the integration of an elastoplastic constitutive law consists in finding solution of the set of equations 2.81, 2.82, and 2.84. The method generally employed is based on a two-step algorithm, termed return mapping, which is described in the following.

Return mapping algorithm

Starting with a strain increment ∆ε, a trial stress is calculated assuming that the strain increment is elastic. This step is termed the elastic predictor :

σ tr = σ n + E : ∆ε (2.
σ n+1 = σ n + ∆σ = σ n + E • ∆ε -E • ∆λ ∂g σ n+θ = σ tr -E : ∂g ∂σ n+θ • ∆λ
Therefore, integrating the constitutive law for a given increment of deformation ∆ε reduces to finding ∆λ that verifies Equations 2.81, 2.82, and 2.83. Most of the time, this set of equations is highly non-linear and an iterative procedure must be performed.

The specific procedure adopted in the present work is named the cutting plane algorithm [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF], [START_REF] Wilkins | Calculation of elastic-plastic flow[END_REF]) and it is developed hereinafter for an arbitrary elastoplastic model.

Let's denote

σ (i) = σ (i) n+1 , κ (i) = κ (i)
n+1 and f (i) = f σ (i) , κ (i) , g (i) = g σ (i) , κ (i) the stress state, internal variable, yield criterion, and the plastic potential surface at the iteration (i).

The stress and internal variables are discretized between the iteration (i+1) and (i) by the forward Euler method (explicit):

σ (i+1) = σ (i) + δσ (i)
(2.90)

κ (i+1) = κ (i) + δκ (i) (2.91)
The increments of stress and internal variable at the iteration (i) are expressed in function of the discrete plastic multiplier as follows:

δσ (i) = δλ (i) .G (i) (2.92) δκ (i) = δλ (i) .h (i) (2.93)
where

G (i) = -E : ∂g (i) ∂σ (i)
(2.94)

h (i) = h σ (i) , κ (i) (2.95)
Finally, the yield surface is linearised around the stress state σ (i) n+1 :

f (i+1) = f (i) + ∂f (i) ∂σ (i) : δσ (i) + ∂f (i) ∂κ (i) : δκ (i) (2.96) = f (i) + δλ (i) • ∂f (i) ∂σ (i) : G (i) + δλ (i) • ∂f (i) ∂κ (i) : h (i) (2.97)
Thence the correction δλ (i) that must be applied in order to recover f (i+1) = 0 is easily computed :

δλ (i) = -f (i) ∂f (i) ∂σ (i) : G (i) + ∂f (i) ∂κ (i) : h (i) (2.98)
and iterative variable are updated using Equations 2.90 and 2.91. If f (i+1) is greater than a given tolerance, the procedure is started again until convergence is reached.

σ n σ n+1 = σ tr = σ (0) (0)
σ (1) σ (2) σ n+1 f n+1 = 0 Cuts (Elastic predictor) Figure 2.8: Geometrical interpretation of the cutting plane algorithm in the stress space, from [START_REF] Ortiz | An analysis of a new class of integration algorithms for elastoplastic constitutive relations[END_REF] The return path followed by the cutting plane algorithm is schematically illustrated in Figure 2.8. At every iteration, the stress and internal variables are projected through planes defined by Equations 2.92 and 2.93. Thence, it may be said that a series of cuts are made in order to minimize the value of the plastic mupliplier (Equation 2.98) hence the name cutting plane algorithm. The optimized value of δλ corresponds to a cut plane that is tangent to the yield surface at the end of the step.

Critical state theory

Critical State Soil Mechanics (CSSM) was developed for the application of the plasticity theory to soil mechanics. The first critical state models were the series of Cam-Clay constitutive laws, developed at the University of Cambridge by Roscoe and his co-workers. The formulation of the Original Cam-Clay model as an elastoplastic constitutive law was presented by [START_REF] Roscoe | On the yielding of soils[END_REF]. Afterwards, [START_REF] Roscoe | On the generalized stress-strain behaviour of wet clay[END_REF] proposed the Modified Cam-Clay model (MCC). CSSM has been widely used since then resulting to the development of various constitutive models. The purpose of all of these models is to achieve a better agreement between the predicted and the observed soil behaviour. In this section, a brief description of the CSSM theory is presented.

Critical state concept

Let us analyse some results of drained monotonic triaxial tests on London clay carried out by [START_REF] Gens | Stress-strain and strength characteristics of a low plasticity clay[END_REF]. Samples of the same soil had been consolidated and sheared at different initial mean effective stresses. Results are presented in Figure 2.9, in terms of stress path in p -q plane, and evolutions of the stress ratio (η = q p ) and volumetric strain function of axial deformation for different values of over consolidation ratio (OCR). The tests show a strong dependency on the mean effective stress. The different behaviours are depicted hereafter.

• OCR < 2: the stress-strain relationship is continuously hardening, the stress ratio η and volumetric strains are increasing. This behaviour means that the elastic locus (or the yield surface) is expanding, and the material is contracting.

• OCR > 2 : the stress-strain relationship shows a smoothed peak on the stress ratio η, followed by a softened behaviour. The volumetric deformation is decreasing, and reaches negative value for the highest OCR. The behaviour is thus dilative.

• All the η -ε a curves tend to the same asymptotic plateau.

These observations are not specific for the aforementioned clay, but have been observed in almost all fine-grained soils (e.g., [START_REF] Gens | Stress-strain and strength characteristics of a low plasticity clay[END_REF], Burland and[START_REF] Burland | On the compressibility and shear strength of natural clays[END_REF][START_REF] Biarez | Elementary mechanics of soil behaviour: saturated remoulded soils[END_REF]). The first test dates back to the late 50s, when Professor Roscoe developed a simple shear apparatus [START_REF] Roscoe | On the yielding of soils[END_REF]), followed by a much more extensive data of triaxial tests from research under the direction of Professor Schofield, which led to the publication of the critical state theory for soils [START_REF] Schofield | Critical state soil mechanics[END_REF]).

The critical state concepts states that, if a soil sample is allowed to change its volume during a shearing test, it will either dilate or contract depending on its initial state of density (initial values of p, q and e). The volumetric yielding process will continue until the soil sample reaches a critical void ratio, after which the soil deforms at a constant void ratio, constant mean effective stress and constant deviatoric stress. This final state is known as the critical state. This behaviour may occur in clays and sands as well, for any drainage and loading conditions.

The critical state starts when clay particles reach a final flow structure that offers the minimum shearing resistance. It was demonstrated that the friction angle mobilised at the critical state is constant. Moreover, the soil behaviour at this state is only function of the critical void ratio and lies on the critical state line in the ln p -e plane. This leads to the following equations:

q = M p (2.99) e = Γ -λ ln p (2.100)
where M is the slope of the critical state line in p -q plane. M is function of the friction angle of the material, Γ is the reference value of the void ratio, and λ the slope of the critical state line in ln p -e plane. As shown in Figure 2.10, the critical state line in ln p -e plane is parallel to the isotropic consolidation line.

Nowadays, the critical state concept is fundamental, especially for finegrained soils. This theory gave birth to numerous constitutive models for clay and sand, such as the well-known models of Original and Modified Cam-Clays described briefly hereafter.

Original Cam-Clay

The original Cam-Clay model was suggested by Roscoe and his co-workers [START_REF] Roscoe | On the generalized stress-strain behaviour of wet clay[END_REF], as a result of the publication of the critical state theory. The formulation of the original Cam-Clay model is based on the following assumptions:

• Elasticity is assumed non-linear, with a pressure-dependent bulk modulus; • The state of the soil is presented with the mean effective and deviatoric stresses (p, q), in addition to the void ratio. The effect of the third invariant of the stress tensor is not included;

• The flow rule is associated;

• The final state of the material during shearing is described with the critical state theory, i.e, when the stress ratio reaches the state η = M , one has εp v = 0 and the stress state does not evolve further;

• The change in size of the yield surface is related to the change of the void ratio (or volumetric deformation);

• Negative pressures are not admissible. Cohesion is hence not included and only the behaviour of reconstituted samples is captured (i.e, the bonding effect is not included);

• Soil is considered isotropic.

Based on these assumptions, the constitutive equations of the Original Cam-Clay are listed as follows:

1. Elastic properties:

It is assumed that recoverable changes in volume vary with the mean effective stress, according to the expression :

εe v = κ 1 + e ṗ p (2.101)
where κ is the slope of the swelling line in ln p -e plane and e is the void ratio of the soil. Equation 2.101 implies a linear relationship between the void ratio and the logarithm of p . Therefore, the bulk modulus is:

K = 1 + e κ p (2.102)
However, the elastic shear strains are negligible in the original Cam-Clay, i.e: εe q = 0 (2.103)

Plastic potential and yield surfaces

The plastic potential is related with the work dissipation. In the triaxial plane, an increment of the work input per unit volume is given as [START_REF] Taylor | Fundamentals of soil mechanics[END_REF]:

Ẇ = q εq + p εv (2.104)
The elastic portion of this total increment is recoverable, whereas the plastic portion, Ẇ p , is available for dissipation, then

Ẇ p = q εp q + p εp v = p η εp q + εp v (2.105)
Within the context of the critical state theory, the plastic work increment Ẇ p is assumed entirely dissipated by friction [START_REF] Schofield | Critical state soil mechanics[END_REF]. Hence, at the critical state, the rate of plastic volumetric deformation is null, and the stress ratio η is equal to M . These conditions lead to:

Ẇ p = M p εp q (2.106)
where M is the slope of the critical state line in p -q plane. Depending on the compression or extension character of loading, M may be positive or negative. The increments of volumetric and deviatoric plastic strains are given in the case of an associated flow rule as:

εp v = λ∂f ∂p , εp q = λ∂f ∂q (2.107)
where the plastic multiplier λ verifies the Equation 2.75, which is written in the triaxial plane as:

λ = 1 K p ∂f ∂p ṗ + ∂f ∂q q (2.108)
where K p is the plastic modulus.

On the other side, combination of Equations 2.105 and 2.106 gives rise to the following equation:

η + εp v εp q = M (2.109)
and given that the dilatancy function is defined as the ratio of the increments of plastic volumetric and deviatoric strains, it yields from Equation 2.109:

Ψ = εp v εp q = M -η (2.110)
According to Equation 2.71, the differential form of the yield surface can be obtained by setting λ = 0, which on view of Equations 2.108 and 2.110 yields:

λ = 0 ⇐⇒ ∂f ∂p ṗ + ∂f ∂q q = 0 ⇐⇒ Ψ ṗ + q = 0 ⇐⇒ 1 M p [(M -η) ṗ + q] = 0 ⇐⇒ 1 M p - q p ṗ + M ṗ + q = 0 ⇐⇒ - q M p 2 ṗ + ṗ p + q M p = 0 ⇐⇒ ∂ ∂p q M p + ln p ṗ + ∂ ∂q q M p + ln p q = 0
Thence, the yield function can be rewritten in the p -q plane as:

f ≡ q M p + ln p + A = 0 (2.

111)

The constant A can be determined using any pair of conjugate values (p in , q in ) on f = 0. Hence, by replacement yielding conditions of the normal consolidation (i.e., p in = p m ,q in = 0, where p m is the preconsolidation pressure) in the Equation 2.111, the final equation of the yield function is obtained: In multi-axial space, the yield surface is expressed as:

f ≡ q + M p ln p p m = 0 (2.
f ≡ 3 2 s : s -M p ln p m p = 0 (2.113)
where s is the deviatoric stress tensor (see Section 2.2).

Hardening rule:

The yield surface is assumed to expand/contract with a constant shape, and the size of the yield surface is assumed to be related to the change in volume only, according to the following equation: .114) where e is the void ratio of the soil, λ and κ are the slope of the normal compression and swelling lines in the ln p -e plane respectively.

ṗm = λ -κ 1 + e εp v ( 2

42CHAPTER 2. HYDROMECHANICAL MODELLING OF GEOMATERIALS

Modified Cam-Clay

One of the major limitations of the original Cam-Clay model is the singular point of the yield surface at its intersection with the hydrostatic axis (the stress state (p m , 0), see Figure 2.11). At this yielding state, the normality condition provides deviatoric plastic strain during an isotropic consolidation, which contradicts with the isotropy assumption.

To avoid this singularity, [START_REF] Roscoe | On the generalized stress-strain behaviour of wet clay[END_REF] proposed a novel yield dissipation equation by including a dependency on the volumetric plastic strain:

Ẇ p = p εp v + qδε p d = p (M εp d ) 2 + ( εp v ) 2 (2.115)
Similarly to the yield surface of the original Cam-Clay, assuming an associated flow rule and combining the work dissipation equation (Equation 2.115) with Equations 2.108, and 2.107 leads to the expression of the yield surface in p -q plane:

f = q 2 + M 2 p (p -p m ) (2.116)
The hardening rule of the MCC is identical to the one of the original Cam-Clay (Equation 2.114).

Limits of the MCC model

The publication of the critical state theory has been a turning point in geomechanics, as it gave birth to the one of the simplest and most used constitutive laws, the MCC model. Indeed, the model is able to capture, by the mean of one state parameter, dilatancy as well as contractancy of clayey soils, and it was widely employed in various geotechnical applications [START_REF] Wroth | Interpretation of in situ soil tests[END_REF]). Nevertheless, the MMC model presents several restrictions regarding the behaviour of clays. First and foremost, the early applications of the MCC to overconsolidated clays have revealed some major drawbacks of the model. Large elastic domain, sudden transition from elastic to plastic region, as well as inadequate predictions of peak failure stress and dilatancy.

Furthermore and as already mentioned, the MCC is restricted to isotropic and reconstituted soils. However, it is well known that the response of natural clays is very different from the reconstituted and isotropic ones. Natural clays are structured materials, their behaviour is dominated by the nature of the fabric and the strength of the inter-particle forces. Moreover, the fabric of natural clays is known to be always anisotropic, which significantly affects the clay behaviour.

In addition, the MCC is inappropriate to capture the cyclic behavior of clays, since a purely elastic behaviour is predicted at each stress reversal.

Therefore, the main objective of the present research work is the development of a constitutive law, based on the MCC model, able to capture the main different features of natural clays, under monotonic and cyclic loading.

Conclusion

The present dissertation deals with the numerical modelling of offshore foundations using the finite element code LAGAMINE. The approach consists in representing the soil (clay) as a saturated porous material, where grains and fluid (water) are assumed to be in thermodynamic equilibrium. Quasi-static analysis is considered, either for monotonic or cyclic loadings.

In the first part of the chapter,the mechanical equilibrium equation is expressed in the usual differential local form and the weak form is written in the current configuration, which corresponds to the updated Lagrangian configuration. The hydraulic equilibrium is based on the equation of mass conservation, written in the local and global forms in function of Eulerian variables. The flow of the fluid mass is described by the Darcy's law. The flow is hence considered laminar, and the soil permeability is assumed isotropic. Finally, the stresses acting on the solid skeleton are the effective Terzaghi stresses. This assumption makes sense in offshore geotechnics, where soil is always saturated.

In the second part, the general finite element formulation of contact problems, as implemented in LAGAMINE, has been recalled. The adopted hydromechanical contact element belongs to the zero-thickness family, which is capable to take into account:

• sliding, captured by the Mohr Coulomb criterion,

• fluid flow parallel to the interface, which is described by the generalized Darcy law;

• fluid flow from the solids through the interface;

• storage flow, which ensures the fluid mass conservation and gives a hydromechanical coupling term.

Special attention has been then paid to geomechanical constitutive modelling. The general framework of elastoplasticity has been presented and specific details on the numerical integration schemes for elastoplastic models have been provided.

The present work focuses on offshore piles founded in clayey soils. Hence, the basic concepts that encounter the constitutive modelling of clays are finally introduced. Details are given on the critical state theory, the formulation of the Original and Modified Cam-Clay models and their restrictions. The soil beneath a structure is henceforth subjected to various stress states as schematically presented in Figure 3.1, which shows loading conditions along a potential failure surface beneath an offshore structure. Along this failure surface, the loading conditions approach either triaxial compression and extension or Direct Shear Stress (DSS) experiments. Therefore, a large range of offshore foundations can be studied in laboratory through the aforementioned experiments.

Figure 3.1: Simplified stress conditions along a potential failure surface in the soil beneath an offshore structure from [START_REF] Andersen | Bearing capacity under cyclic loading-offshore, along the coast, and on land. the 21st bjerrum lecture presented in oslo, 23 november 2007[END_REF] One of the first laboratory study that respects the in-situ conditions dates back to 1930's when Casagrande in Harvard University developed, under the direction of Terzaghi, the triaxial device for soils Terzaghi (1936a). Thereafter, experiments in the field of geomechanics had taken a turning point, which resulted in a deep understanding of the soil behaviour and development of different design models and constitutive laws.

The monotonic behaviour of clays depends strongly on various factors such as the initial mean effective stress, the orientation of principal axes and the loading history. Moreover, properties such as strength, permeability and compressibility, depend on the size of soil particles, their shape and arrangement (e.g., [START_REF] Biarez | Elementary mechanics of soil behaviour: saturated remoulded soils[END_REF], [START_REF] Hicher | Comportement mécanique des argiles saturées sur divers chemins de sollicitations monotones et cycliques application à une modélisation élastoplastique et viscoplastique[END_REF], [START_REF] Locat | The compressibility and sensitivity of an artificially sedimented clay soil: The grande-baleine marine clay, québec, canada[END_REF]). Natural clays are known to exhibit a high degree of anisotropy due to their mode of deposition and the subsequent loading. Neglecting this feature may lead to inaccurate prediction of soil response [START_REF] Ahayan | L'effet de l'écrouissage rotationnel sur le comportement des sols argileux[END_REF]). Furthermore, natural clays may show a significant degree of sensitivity, which is quantified by the loss of strength when the naturally deposited soil is disturbed. Generally, the natural structure of intact soils renders their behaviour different from that of the reconstituted material. The sensitivity of clays has been the reason of many instability problems and landslides [START_REF] Bjerrum | Engineering geology of norwegian normally-consolidated marine clays as related to settlements of buildings[END_REF]).

On the other hand, naturally deposited clays may also exhibit certain degree of overconsolidation due to their past loading, such as erosion, excavation, as well as water table variation and cyclic loading. A distinct behaviour of the overconsolidated (OC) clays is related to their volume change when subjected to shearing. Typically, a heavily OC clay shows volume expansion, known as dilatancy, while a lightly OC clay exhibits volume contraction when subjected to shearing [START_REF] Biarez | Elementary mechanics of soil behaviour: saturated remoulded soils[END_REF].

For an accurate prediction of clay behaviour is therefore necessary to capture anisotropy, sensitivity and dilatancy. In this chapter, the main features of clays are firstly presented, notions such as structure, fabric, and bonding are explained. Afterwards, the principal properties of natural clays are defined and described. The essential characteristics of dilatancy, anisotropy, and sensitivity are developed, and modelling concepts are reviewed.

Main features of clays

Mineralogy

Before presenting the main properties of clay soils, the principal features of clay minerals are introduced.

Clay minerals belong to the phyllosilicate group, characterized by layered structures composed of polymeric sheets of silica tetrahedra attached with octahedral sheets.

As shown in Figure 3.2 (a), the tetrahedral sheet, named also silica sheet, is composed of silica tetrahedra (SiO 4 ) 4-, which are linked together by sharing three of their four oxygen ions. All silica sheets are oriented in the same direction and define hexagonal cavities. Likewise, Figure 3.2 (b) depicts the structure of the octahedral sheet, called also alumina sheet, which is made up of aluminium or magnesium octahedral. The cation of the alumina sheet bonds with six oxygen atoms or hydroxyl groups. Octahedral sheets are all linked together in a triangular face, by sharing their six oxygen or hydroxyls.

The association of tetrahedral and octahedral sheets form a layer. Two types of layers are defined, namely the TO and the TOT layer. The TO layer, consists on one tetrahedral and one octahedral sheets. The TOT layer is com- posed of an octahedral sheet sandwiched between two tetrahedral sheets. In both cases strong primary bonds exist between sheets of a same layer [START_REF] Stçepkowska | Aspects of the clay/electrolyte/water system with special reference to the geotechnical properties of clays[END_REF], but the bonds between two successive layers are generally weaker.

The properties of clay minerals, particularly in presence of water, depend on the type of layers and the nature of the interlayering. On this basis, clay minerals are classified into different groups, as presented in Figure 3.3. The structure of Kaolinite is built on the TO layer, while Illite and Smectite are based on the TOT layer. Kaolinite is characterized by hydroxyl groups and oxygen ions in two consecutive layers, which are bonded through hydrogen bonds. The bonds maintain the interlayer closed and prevent swelling. In Illite, the layers are strongly locked together through ions bonds, that form the interlayer space. Conversely, the interlayers of Smectite are easily penetrated by water molecules and ions, leading to the expansion of the mineral upon wetting.

Finally, the layers stack together to form particles. The orientation and distribution of these particles in a specific soil mass define the fabric of the clay.

Structure: bonding and fabric

Even though clays are composed of discrete particles and particle groups, clays are often treated as continuum media in engineering analysis. However, prop-Figure 3.3: Classification of clay minerals after [START_REF] Mitchell | Fundamentals of Soil Behaviour[END_REF] erties such as strength, permeability and compressibility, depend on the size and shape of the particles, their arrangement and the force between them [START_REF] Mitchell | The fabric of natural clays and its relation to engineering properties[END_REF]). In other words, the soil properties depend on its structure. The term structure means the combination of bonding, the inter-particle forces and fabric, the arrangement of the component particles.

Bonding

The inter-particle bonding is defined as the combination of all inter-particle forces. They can be of electrostatic or electromagnetic nature. In general terms, bonding gathers all factors acting to keep the soil particles together [START_REF] Gens | Stress-strain and strength characteristics of a low plasticity clay[END_REF]).

Two different types of bonds can be distinguished in a clay mass: interatomic bonds within clay particles and inter-particle bonds between particles and particles groups. The former is generally stronger than the latter. Hence, soil mass can be considered as an assemblage of particles in which deformation processes are dominated by displacements between particles and not by deformations of the particles themselves.

The structure of clay results hence from a physico-chemical equilibrium between particles and particles groups. This equilibrium develops during the geological life of the soil as a result of electrostatic and magnetic interaction between crystals, ions and water concentration, temperature during deposition, in addition to osmotic pressure and organic content [START_REF] Biarez | Elementary mechanics of soil behaviour: saturated remoulded soils[END_REF]).

External variation of these factors or the development of chemical reactions within the sediments can cause substantial changes of both fabric and bonding of the clay during time. The soil structure is henceforth considered as the result of all the processes that the soil has undergone during its geological history.

Furthermore, the soil structure is considered evolving with loading. According to [START_REF] Leroueil | Stress-strain rate relation for the compressibility of sensitive natural clays[END_REF], any soil, at any state, can be considered to have a structure. The authors distinguished an intact state of structure, a remoulded state and a destructured one. The first state occurs in natural deposits, and is characterised by maximum bonding forces and strength. The second state is obtained when sufficient mechanical energy is spent to reduce the strength to a minimum. The third state corresponds to the state of reconstituted soils, which are subjected to large strains that destroy totally its original structure. [START_REF] Burland | On the compressibility and shear strength of natural clays[END_REF] defined the properties of a reconstituted soil as intrinsic, referring to the basic or inherent properties of the material.

In the present dissertation, the reconstituted state will be considered as a reference state to define the intrinsic parameters of natural clays.

Fabric

The term fabric refers to the size, shape and arrangement of the solid particles, the organic inclusions and the associated voids of a soil mass. [START_REF] Mitchell | The fabric of natural clays and its relation to engineering properties[END_REF] defined fabric as the appearance of pattern produced by the particles shape and described the geometrical interrelationships among soil particles with respect to the local and general degree of orientation of particles.

Knowledge of the history of deposition and soil sedimentation is a goodstarting point for understanding how soil fabrics are formed and changed throughout the history of deposition. Accordingly, [START_REF] Burland | On the compressibility and shear strength of natural clays[END_REF] explained that slow deposition in still water leads to an open fabric, while rapid deposition gives rise to more oriented fabric, which is consequently more compact. Hence, different types of fabric could be reorganized depending on the mode of deposition and sedimentation. Sides and [START_REF] Sides | The microstructure of dispersed and flocculated samples of kaolinite, illite, and montmorillonite[END_REF] provided a classification of fundamental fabrics (see Figure 3.4), based on the shape, arrangement and bonding network of the soil skeleton. They identified two main groups, termed the flocculated and the dispersed fabrics. In the former, clay particles tend to move toward each other and become attached, while in the later particles tend to move away due to repulsive inter-particle forces. For the flocculated fabric, various subgroups are distinguished as depicted in Figure 3.4, such as (i) the Cardhouse fabric, when a single particle arrangement is detected, (ii) the Bookhouse fabric, when particles arranged in groups, called domains, are parallel, and (iii) the Honeycomb fabric, characterized by non-uniform domains occurring during compression. For Dispersed fabrics, they defined a Turbostatic structure, when domains are not oriented into the same direction. Stacks is when the fabric is characterized by a highly oriented domains.

The fabric significantly influences the mechanical response of clays. According to [START_REF] Oda | Initial fabrics and their relations to mechanical properties of granular material[END_REF], soil specimens with different sample preparation methods have different fabric and, expectedly, respond differently to loading. It is however necessary to keep in mind the scale at which the fabric is of interest. For instance, the behaviour of an homogeneous soft clay is expected to be greatly influenced by the particle arrangement at a micro-scale, whereas a stiff and intact clay where some macro-cracks are present, is mainly controlled by properties of the crack than the nature of its fabric [START_REF] Han | Étude expérimentale du comportement mécanique d'une argile fortement surconsolidée sous chargements monotones et cycliques[END_REF]).

Overconsolidated clays and modelling

Naturally deposited clays may exhibit a certain degree of overconsolidation due to their past loading history such as cyclic loading, repeated compaction, excavation and backfilling, as well as water table variation. Other sources of overconsolidation include phenomena such as desiccation, secondary compression, and cementation. Traditionally, the overconsolidation degree is expressed by a single scalar, the OverConsolidation Ratio OCR = σp /σv, where σ p is the preconsolidation pressure (maximum past vertical effective stress), measured from one-dimensional consolidation tests and σ v is the current vertical effective stress. An alternative definition that uses the mean effective stress instead of vertical stresses is also usually employed. OCR is a basic element of normalized soil behaviour as described by the critical state soils mechanics. A distinct feature of an OverConsolidated (OC) clay is its dilatancy behaviour, which depends closely on the degree of overconsolidation. Typically, a heavily OC clay shows volume expansion while a lightly OC clay exhibits volume contraction when subjected to shearing. Furthermore, OC clays exhibit a higher degree of non-linearity compared to the normally consolidated soils. According to [START_REF] Jardine | The measurement of soil stiffness in the triaxial apparatus[END_REF], the majority of OC soils are highly non-linear even at axial strain levels as small as ε a = 0.005%.

There have been numerous efforts to capture the behaviour of OC clays in the early 1960's, most of which based upon the original or the Modified Cam-Clay models. Representative works include the bounding surface or the subloading surface models that extend the theory of elastoplasticity to consider inelasticity within the yield surface. There exist also attempts to model the behaviour of OC clays using the classical theory of elastoplasticity, by replacing the yield surface of the MCC in the supercritical region (dry side) by an appropriate surface.

The main features of the behaviour of OC clays are reviewed in the following sections. Several experimental studies are firstly presented and then details are given about several modelling aspects.

Main features

Dilatancy

Dilatancy was first put forward by [START_REF] Reylonds | On the dilatancy of media composed of rigid particles in contact with experimental observations[END_REF], when he demonstrated that the volume change of a regular packing of rigid spheres under shear depends on the initial arrangement of the spheres and that contraction would occur in a loose packing, whereas dilatancy would occur in a dense packing. Consequently, the shear strength of the sphere assembly depends not only on the friction between the spheres, as assumed up to that time, but also on the volume change (see Figure 3.5). Following this pioneer work, several experimental studies on remoulded soils, particularly sands, improved our understanding of the phenomenon. Results of sand drained triaxial tests clearly showed that the initial density is an important parameter having a great influence on the dilatancy rate. In undrained triaxial tests, the dilatant behaviour of dense sand leads to a pore pressure decrease while a pore pressure increase is observed in loose samples.

For clay materials, the results of drained and undrained tests showed that the volume change (or pore pressure evolution) during shearing mainly depends on the clay loading history, characterized by the OCR. The clay behaviour is relatively similar to sand in this respect; contraction or dilatancy in sand are directly connected to the relative density, while in clays these volumetric changes are more precisely related to the OCR [START_REF] Hattab | Dilating behaviour of overconsolidated clay[END_REF].

More specifically, drained and undrained triaxial tests published by [START_REF] Henkel | The effect of overconsolidation on the behaviour of clays during shear[END_REF] showed that the volume change for clayed materials mainly depends on the OCR. Normally consolidated or slightly overconsolidated clay samples showed a contractant behaviour, while highly overconsolidated samples presented a dilatant behaviour. Similar results were obtained by [START_REF] Zervoyannis | Etude synthétique des propriétés mécaniques des argiles saturées et des sables sur chemin oedométrique et triaxial de révolution[END_REF], who showed that dilatancy phenomena appeared generally for an OCR greater than 2 or 2.5, depending on the clay mineralogy. [START_REF] Shimizu | Effect of overconsolidation on dilatancy of a cohesive soil[END_REF] performed triaxial tests at constant mean effective stress on Fujinanori clay with different overconsolidation ratios. The author showed that the stress plane is divided in several regions of contractant and dilatant behaviours according to the stress ratio and the OCR (see Figure 3.6). These first experimental studies confirmed that the dilatancy of OC clays is mainly dependent on the degree of overconsolidation.

Non-linearity at small strain levels

Most overconsolidated soils are highly non-linear even at small strain levels [START_REF] Jardine | The measurement of soil stiffness in the triaxial apparatus[END_REF]. Heavily OC clays showed a highly non linear, but very stiff, initial behaviour. A marked loss of stiffness was generally obtained at an early stage of almost 0.1% axial strain. Following the same authors, the stiffness ratio Eu /su (where E u is the undrained Young Modulus and s u the undrained shear strength of the material) was strongly dependent on the OCR for the intact and reconstituted samples. As depicted in Figure 3.7, lightly overconsolidated samples produced the lowest values of the ratio, while heavily overconsolidated and remoulded samples the highest.

Inadequacy of the Modified Cam-Clay model

The critical state concept and the Cam-Clay family of models are considered as a turning point in the field of geomaterials, since they attempt to simulate volume change and friction effects using one simple isotropic hardening rule. The Modified Cam-Clay (MCC) model was the first constitutive law for soils capable of reproducing the contractant and dilatant behaviour. From the early computational applications it became however evident that the MCC model presents some limitations/drawbacks.

Constitutive laws for overconsolidated clays

A simplified way to model the behaviour of OC clays is to use the classical theory of elastoplasticity by changing the yield surface of the MCC in the supercritical region (dry side of the yield surface) by an appropriate surface. Another approach is the bounding surface or the subloading surface models that extend the theory of elastoplasticity to include inelasticity within the yield surface. They are therefore capable of describing the gradual transition from elastic to plastic region, the behaviour of overconsoslidated soils as well as the cyclic features of soil behaviour. These different concepts are explained hereafter:

• Appropriate yield surface A first succesful attempt to capture realistic values of the undrained shear strength was proposed by [START_REF] Zienkiewicz | Finite element studies of soils and porous media[END_REF]. The authors replaced the Cam-Clay surface in the dry side with the Hvorslev failure surface, which imposes a more accurate limit to admissible stress states above the critical state line providing thus a more realistic description of the peak shear stress value (see Figure 3.9). [START_REF] Zienkiewicz | Finite element studies of soils and porous media[END_REF] Following a similar idea, Whittle and co-workers [START_REF] Whittle | Evaluation of a constitutive model for overconsolidated clays[END_REF]) proposed a more appropriate yield surface that reduces the elastic domain in the supercritical region. More recently, [START_REF] Collin | Couplages thermo-hydro-mecaniques dans les sols et les roches tendres partiellement satures[END_REF] and [START_REF] Jin | Unified modeling of the monotonic and cyclic behaviors of sand and clay[END_REF] proposed a yield surface that combines the Cam-Clay surface in the wet side with the Mohr Coulomb line in the dry side.

Deviatoric stress q (kPa)

• Small-strain stiffness

The aforementioned modified constitutive laws are able to capture more accurately the peak strength of OC clays, they are however unable to reproduce non-linearities at small stress levels. To improve this, several authors modified the formulation of the elastic strain components. For example, [START_REF] Hueckel | Plastic-flow of granula and rocklike materials with variable elasticity moduli[END_REF] proposed an elastoplastic coupling formulation in which the elastic moduli are related to the plastic deformation. [START_REF] Hueckel | Some hysteresis effects of the behaviour of geologic media[END_REF], [START_REF] Hardin | Shear modulus and damping in soils: design equations and curves[END_REF], [START_REF] Biarez | Elementary mechanics of soil behaviour: saturated remoulded soils[END_REF] among others, introduced a more elaborate piecewise elastic formulation, which describes the small strain non-linearity of clays and sands and models the hysteretic unload-reload behaviour.

• Extended elastoplastic theories Evaluation of the foregoing models shows that while they capture some aspects of the observed behaviour of overconsolidated clays, no single model describes fully the observed aspects summarized above [START_REF] Hicher | Comportement mécanique des argiles saturées sur divers chemins de sollicitations monotones et cycliques application à une modélisation élastoplastique et viscoplastique[END_REF]. This challenge has been achieved when Mroz proposed a two surfaces model with kinematic hardening [START_REF] Mroz | An anisotropic hardening model for soils and its application to cyclic loading[END_REF] that is able to predict plastic strains for stress state within the yield surface, while respecting the general framework of elastoplasticity. Subsequently, various approches have been introduced following the Mroz concept, such as multisurface models [START_REF] Prevost | Plasticity theory for soil stress-strain behavior[END_REF], bounding surface models Dafalias (1986b), and subloading surface models [START_REF] Hashiguchi | Elastoplastic constitutive equation of soils with the subloading surface and the rotational hardening[END_REF]. These models describe well the gradual transition from elastic to plastic region, as well as the soil behaviour under cyclic loading.

Bounding Surface plasticity for monotonic behaviour

The general aspects of bounding surface plasticity associated with rate independent plasticity are presented in this section. The material state is defined in terms of the effective stress tensor σ and an inelastic internal variable κ, which considers proper measures of inelastic deformations. The bounding surface is defined analytically in the stress space by:

F (σ, κ n ) = 0 (3.1)
where a bar over a stress quantity (i.e. (•)) indicates an image point on the bounding surface. The actual stress point lies always within or on the bounding surface. To each σ a unique image stress is assigned by a properly defined mapping rule (which becomes the identity mapping if σ is on the surface).

The radial mapping rule Dafalias (1986b), [START_REF] Hashiguchi | Elastoplastic constitutive equation of soils with the subloading surface and the rotational hardening[END_REF], shown schematically in Figure 3.10 is widely used in the framework of Bounding Surface Plasticity. Using σ c as the projection center, σ is obtained by radially projecting σ on the bounding surface, hence the name radial mapping rule. Analytically, this projection is expressed by:

σ = σ c + b (σ -σ c ) (3.2)
The definition of the projection center is specific to each constitutive law. In the case of monotonic loading, σ c is generally considered fixed at the origin of the stress.

The variable b (b ≥ 0) is determined in terms of the stress state and the internal variables of the model, by substituting σ from Equation 3.2 onto an explicit form of Equation 3.1 and solving for b. As shown in Figure 3.10, a consequence of assuming a radial mapping rule is that a loading surface homologous to the bounding surface with respect to σ c and passing through σ is indirectly defined. As for the case of rate independent associated plasticity, the direction and magnitude of the plastic strain rate are function of the plastic potential G derived at the image stress, and the rate of the plastic multiplier λ:

εp = λ ∂G ∂σ (3.3)
where .. are the Macauley brackets. If an associate flow rule is considered, the direction of the plastic potential will be identical to the direction of the bounding surface.

Respecting the general framework of the elastoplasticity, the rate of plastic multiplier is given as follows

λ = 1 K p ∂F ∂σ σ (3.4)
where K p is the plastic modulus at the current state of stress σ. For a stress state that lies on the bounding surface (i.e., σ = σ), the plastic modulus at the image stress K p = K p (σ) is defined by applying the consistency condition with respect to the image stress, i.e. Ḟ (σ, κ) = 0. For a stress state that lies within the bounding surface, an equivalent relationship is defined by [START_REF] Kaliakin | Theoretical aspects of the elastoplasticviscoplastic bounding surface model for cohesive soils[END_REF]:

1 K p ∂F ∂σ σ = 1 K p ∂F ∂σ σ (3.5)
Equation 3.5 is the crux of the bounding surface plasticity. It describes how plasticity occurs within the bounding surface and serves also to determine the rate equation of b.

As mentioned before, the plastic modulus K p at the image stress point is defined via the consistency condition (as in classical plasticity). Within the bounding surface, K p is a function of K p and the distance to the bounding surface. This relation is analytically expressed by:

K p = K p + H (3.6)
where H denotes a proper scalar hardening function of the state, named shape function. This function defines the shape of the response curve during plastic hardening (or softening) for points within the bounding surface and relates the plastic modulus K p to its bounding value K p . Generally, H is an increasing function of b (or the Euclidean distance δ between σ and σ, see Figure 3.10). Due to the Equation 3.6, a smooth transition from the purely elastic to purely plastic response is captured. When the stress state lies on the bounding surface, the current and image stresses are identical, and hence, the shape function is null (for b = 1 or δ = 0), and the bounding surface model is reduced to the classical elastoplastic model. Equation 3.6 embodies the meaning of the bounding surface concept. If δ < r or b > 1 and H is not approaching infinity, the concept allows for plastic deformations to occur for points either within or on the surface at a progressive rate that depends upon δ (or b). The closer to the bounding surface is the actual stress σ, the smaller is K p and the greater is the magnitude of plastic strain increment for a given stress increment. The stress state may eventually reach the bounding surface in the course of plastic loading and then it remains on the surface, if the loading continues. Thus, for states within the bounding surface the function H and its associate parameters and constants are intimately related to the material behaviour response. As such, they constitute important new ingredients of the bounding surface model with regard to classical elastoplasticity.

Table 3.1 summarizes shape functions of some constitutive models that employ the bounding surface concept. Dafalias (1986b), [START_REF] Jiang | An associative and non-associative anisotropic bounding surface model for clay[END_REF] and [START_REF] Huang | Simulation of yielding and stress-stain behavior of shanghai soft clay[END_REF] among others, proposed simple forms of H that introduce an elastic nucleus, homologous and smaller in size than the loading surface. The size of the elastic nucleus is controlled by the constant s and the magnitude of the plastic modulus is controlled by the constant H 0 . Jiang et al.

(2012) added an additional power constant W to capture more realistically the plastic strains within the bounding surface. It is also possible to assume an elastic nucleus limited to a point. This is the case of the shape function of [START_REF] Yin | Modeling time-dependent behavior of soft sensitive clay[END_REF], which do not consider a maximum value of the ratio b before starting yielding within the bounding surface and employs only one additional parameter H 0 controlling the magnitude of the shape function. There is also the possibility to capture an exact formulation of H without employing additional constants, by introducing state parameters that depend on b, δ, the current stress σ, the void ratio, and so on [START_REF] Ahayan | Loi de comportement à surface limite pour les argiles naturelles[END_REF]). [START_REF] Yao | Modified uh model: constitutive modeling of overconsolidated clays based on a parabolic hvorslev envelope[END_REF] for example, proposed a unified model, based on the bounding surface plasticity and the original Cam-Clay model, that captures the behaviour of overconsolidated soils with the same number of parameters that the classical Cam-Clay model. To do so, the authors have related the overconsolidated ratio (OCR), which is related to the parameter b of the mapping rule, to the Hvorslev surface. In the same spirit, [START_REF] Jocković | Bounding surface model for overconsolidated clays with new state parameter formulation of hardening rule[END_REF] proposed a modified form of the MCC, based on the bounding surface concept, and introduced a shape function relating the overconsolidation ratio to a well-defined state parameter Ψ. Table 3.1: Some Bounding surface models for the behaviour of overconsolidated clays ,

Shape function H Constants Bounding surface

Hardening Reference

H = H 0 p 3 m δ r-sδ H 0 , s Modified Cam-Clay Isotropic Dafalias (1986b) H = H 0 δ r-sδ W H 0 ,s,W Rotated MCC Isotropic, Anisotropic Jiang et al. ( 2012
)
H = H 0 δ 0 δ 0 -δ ψ -1 H 0 ,δ 0 ,ψ Rotated MCC Isotropic, Anisotropic, Destruc- turation Huang et al. (2011) H = H 0 p 3 m (b -1) H 0 Rotated MCC Isotropic, Anisotropic Yin et al. (2011)
H = 2 1+e λ-κ p m (bp) 2 1 -1 bω - Original Cam-Clay Isotropic Yao et al. ( 2012
)
H = 2 1+e λ-κ M 4 p m (bp) 2 (1 -ω) - Modified Cam-Clay Isotropic Jocković and Vuk- ićević (2017)

Anisotropy and modelling

Natural clays are commonly anisotropic due to their mode of deposition, i.e. their mechanical properties depend on the orientation of stress and strain tensors. The importance of soil anisotropy is well known in geotechnical engineering [START_REF] Ahayan | L'effet de l'écrouissage rotationnel sur le comportement des sols argileux[END_REF]). For instance, Bjerrum (1973b) and [START_REF] Jardine | Evaluating design parameters for multi-stage construction[END_REF], among others, drew attention to the principal stress axis rotation in soft clays under embankment loading, and highlighted the need to better understand the anisotropy of shear strength. In foundation design, [START_REF] Leroueil | Discussion on: Correlations between index tests and the properties of remoulded clays[END_REF] demonstrated that the failure mode strongly depends on soil anisotropy. The author demonstrated analytically that the bearing capacity of a strip foundation, in an elastic-perfectly plastic soil mass, may reduce or increase up to 40% depending on the strength directionality.

Hence, the effect of anisotropy on the clay behaviour may become critical. Several parameters such as the undrained strength, the consolidation coefficient and elastic deformability may exhibit a strong directional preference. This feature is the most important single cause of discrepancy between classical numerical modelling and field results. Soil anisotropy is mentioned by [START_REF] Biarez | Elementary mechanics of soil behaviour: saturated remoulded soils[END_REF] as an uttermost research need in geotechnical engineering. Neglecting soil anisotropy is a major shortcoming of most in-situ measurement techniques.

Inherent and Stress-Induced Anisotropy of clay

The deposit of soil particles occurs under gravity and hence it is directionaly dependent. The soil fabric is therefore anisotropic, which explains the differences in soil response depending on the direction of loading. For example, non-spherical particles tend to deposit with their long axis in the direction perpendicular to gravity. This mechanism leads to inherent or structured anisotropy. According to [START_REF] Gens | Stress-strain and strength characteristics of a low plasticity clay[END_REF], the inherent anisotropy arises from the structure of the soil as a consequence of geological processes, created by the application of an anisotropic stress tensor during sedimentation. [START_REF] Ladd | Stress-strain behaviour of anisotropically consolidated clays during undrained shear[END_REF] assumed that the inherent anisotropy is caused by the initial one-dimensional deposition, named K 0 consolidation. Anisotropy might also exist in some reconstituted soils if they have undergone an anisotropic plastic strain history. In this case, anisotropy is known as stress (or strain) induced anisotropy and results on the anisotropy of the current stress conditions [START_REF] Biarez | Elementary mechanics of soil behaviour: saturated remoulded soils[END_REF].

The notions of inherent and induced anisotropies were first introduced by [START_REF] Casagrande | Shear failure of anisotropic materials[END_REF] as the result of an extensive experimental investigation on Boston clay. This classification is conceptually elegant because it separates the effects of soil structure and pre-straining, the former leading to the inherent anisotropy, while the later causes the induced anisotropy. However, the afore-mentioned classification provides little help in the experimental quantification of soil anisotropy. Inherent and induced anisotropy are usually measured together, named combined anisotropy, due to the difficulty in measuring the two components separately and because engineers are usually interested in their combined effects.

In the following, the term anisotropy refers to the combined anisotropy and describes the directional dependency of the stress-strain relationship at a soil material point. The effect of anisotropy on clay behaviour and different aspects of the way to model this feature are discussed in the next two sections.

Effects

According to numerous experimental studies (e.g., [START_REF] Leroueil | Behavior of destructured natural clays[END_REF], [START_REF] Leroueil | Stress-strain rate relation for the compressibility of sensitive natural clays[END_REF] and [START_REF] Wood | True triaxial tests on boston blue clay[END_REF]), anisotropy may affect the undrained and drained behaviour of clays. This feature has a pronounced effect on the soil stress-strain response and strength. For instance, [START_REF] Bishop | The influence of the size and orientation of the sample on the apparent strength of the london clay at maldon, essex[END_REF] and [START_REF] Wesley | Influence of stress path and anisotropy on the behaviour of a soft alluvial clay[END_REF], among others, showed that the undrained shear strength of specimens trimmed along the consolidation direction is higher in compression than in extension. In addition, according to [START_REF] Stipho | Theoretical and experimental investigation of the behavior of anisotropically consolidated Kaolin[END_REF], the pore water pressure is significantly affected by the extent of initial anisotropic consolidation. [START_REF] Ladd | Stress-strain behaviour of anisotropically consolidated clays during undrained shear[END_REF] observed that the yield surface of normally consolidated clays is oriented along a line close to the K 0 -consolidation line.

In the following, the effect of anisotropy on the yield and plastic potential surfaces and the undrained shear strength of natural clays is discussed through past experimental studies. It is well known that the undrained shear strength is anisotropic, although its exact range or magnitude differs for various soils. Figure 3.11 shows the variation of the undrained shear strength for four soft Norwegian clayey soils. It is shown that when the angle θ between the major principal stress and the vertical direction is equal to zero, the undrained shear strength reaches maximum values. However, it continuously decreases as θ increases until it reaches its minimum value at θ = 90 • . [START_REF] Ladd | Stress-strain behavior of saturated clay and basic strength principles[END_REF], [START_REF] Ladd | Stress-strain behaviour of anisotropically consolidated clays during undrained shear[END_REF], [START_REF] Henkel | The influence of stress history on stress paths in undrained triaxial tests on clay[END_REF], [START_REF] Donaghe | Effects of anisotropic versus isotropic consolidation in consolidated-undrained triaxial compression tests of cohesive soils[END_REF], among others, studied the effect of K 0 consolidation on the clay behaviour compared to the isotropic consolidation. They observed that the K 0 consolidation had a significant effect on the undrained shear strength with differences up to 10 or 15%. [START_REF] Mayne | Stress anisotropy effects on clay strength[END_REF] reviewed published data for over 40 different normally consolidated and overconsolidated clays subject to isotropic or anisotropic consolidation before shearing. He found that the anisotropic undrained shear strengths for triaxial compression and extension were respectively about 87% and 60% of the respective isotropic strengths.

Effect on undrained shear strength

Effect on yield and plastic potential surfaces Tavenas and [START_REF] Tavenas | Effects of stresses and time on yielding of clays[END_REF] experimentally studied the behaviour of Champlain clay deposits and more specifically the effects of the initial K 0 consolidation on the shape of the yield surface. The results are shown in Figures 3.12 (a) and (b) in stress space. The shape of the experimental yield locus differs strongly from the theoretical shape of the yield locus of the isotropic Cam-Clay models. The locus is not centred on the hydrostatic axis, but rather on a line close to the K0-line. [START_REF] Wood | True triaxial tests on boston blue clay[END_REF] found also similar results from laboratory measurements.

Similar conclusions have been also reported recently from an experimental study on Shanghai clay [START_REF] Gao | Compression and extension yield of an anisotropically consolidated soil[END_REF]. A series of drained triaxial tests under constant stress ratio were conducted by the authors, their results are depicted in Figure 3.13 in terms of stress paths and ε v -ε d curves. The tests clearly show the anisotropic behaviour of Shanghai clay. Shear strains are not symmetric with respect to the loading direction, i.e. compression produces less shear deformations than extension. Furthermore, shear plastic strains are observed during the isotropic consolidation (the curve corresponding to η = 0), which means that the plastic potential surface and the yield surface (in the case of an associated flow rule) are not centred on the hydrostatic axis.

Modelling

A considerable effort has been put from the late 60's to develop constitutive models accounting for the anisotropic behaviour of natural clays. The first attempt dates back to 1970, when [START_REF] Prevost | Plasticity theory for soil stress-strain behavior[END_REF] proposed an anisotropic model for soils as an extension of the work of Mroz for metals [START_REF] Mroz | An anisotropic hardening model for soils and its application to cyclic loading[END_REF]). The Prevost's model inspired many authors [START_REF] Cerfontaine | Possibilities and limitations of the prevost model for the modelling of cohesionless soil cyclic behaviour[END_REF], see for example the model of [START_REF] Banerjee | A plasticity model for the mechanical behaviour of anisotropically consolidated clay[END_REF], with a kinematic hardening and an inclined yield surface to capture anisotropy.

The next section provides a review of some well-known theories and concepts to model the anisotropic behaviour of soils, particularly clays.

Elastic anisotropy

The concept of a soil as an anisotropic elastic medium has been considered by many researchers (e.g., [START_REF] Barden | Stresses and displacements in a cross-anisotropic soil[END_REF], [START_REF] Hooper | Elastic settlement of a circular raft in adhesive contact with a transversely isotropic medium[END_REF]). The stress-strain relationship is given by a 6x6 matrix of moduli (the elastic stiffness matrix).

However, given that the stiffness matrix is symmetrical [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF]), a total of 21 independent parameters are needed to describe a fully anisotropic elastic behaviour.

Nevertheless, soils often show limited forms of anisotropy. For example, a transversal isotropic material possesses an axis of symmetry, i.e., its properties are independent of rotation around the axis of symmetry. Transverse isotropy characterizes almost all granular soils and particularly lightly consolidated post-glacial clays, as a result of their mode of deposition [START_REF] Graham | Anisotropic elasticity of a natural clay[END_REF]). A soil deposited vertically and then subjected to equal horizontal stresses is expected to exhibit a vertical axis of symmetry and therefore be transversely isotropic.

According to [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF] and [START_REF] Graham | Anisotropic elasticity of a natural clay[END_REF], the elastic behaviour of a transversely isotropic material is fully described by five independent parameters. Under triaxial test conditions, with the axis of symmetry of the test corresponding to the vertical direction, i.e., on a vertically cut specimen, only three parameters are needed to define elasticity. Anisotropy is therefore in this case expressed by multiplying the stiffness terms measured in the vertical direction by an anisotropy coefficient, termed β. Following this concept, [START_REF] Graham | Anisotropic elasticity of a natural clay[END_REF] and [START_REF] Houlsby | Study of plasticity theories and their applicability to soils[END_REF] described the elasticity of a transverse isotropic material using the following matrix:

C e = Ev (1+ν vh )(1-2ν vh )           1 -ν vh βν vh βν vh βν vh β 2 (1 -ν vh ) β 2 ν vh βν vh β 2 ν vh β 2 (1 -ν vh ) β 2 (1 -2ν vh ) β (1 -2ν vh ) β (1 -2ν vh )          
(3.7) where E v and ν vh are the Young Modulus and Poisson's coefficient measured in the vertical and horizontal directions, respectively.

Hill's anisotropic theory

Hill's anisotropic theory or Hill's first elastoplastic criterion [START_REF] Hill | The mathematical theory of plasticity[END_REF]) provides a general description of materials with anisotropy and orthotropic symmetry (i.e., materials where each point possesses three mutually normal planes). The failure criterion is of the following form:

σ 11 X 2 + σ 22 Y 2 + σ 33 Z 2 - 1 X 2 + 1 Y 2 - 1 Z 2 σ 11 σ 22 + 1 Y 2 + 1 Z 2 - 1 x 2 σ 22 σ 33 + 1 Z 2 + 1 X 2 - 1 Y 2 σ 22 σ 33 + τ 12 S 2 + τ 23 T 2 + τ 13 R 2 = 1 (3.8)
where X, Y, Z are tensile yield stresses and R, T, S are the shear yield stresses with respect to the axes of anisotropy. Equation 3.8 is an extension of the Von Mises criterion from isotropic to orthotropic material. In the same spirit, [START_REF] Nova | A generalized failure condition for orthotropic solids[END_REF] and [START_REF] Yu | Constitutive modelling of swedish cohesice soils accounting for anisotropy[END_REF] extended the MCC yield function from an isotropic to a transversal isotropic material.

Equivalent shear stress

In equivalent shear stress models, the stress ratio η = q p is replaced by an equivalent stress ratio η = q p -η 0 , where η 0 is the stress ratio at the end of consolidation. Hence, the model allows rotation of the yield surface during consolidation, but no further rotation occurs during subsequent shearing. The equivalent shear stress models were first proposed by [START_REF] Sekiguchi | Induced anisotoropy and time dependency in clay[END_REF], based on experimental observations of the yielding behaviour of clays. It was also adopted afterwards to a series of MIT models (e.g., MIT-E1 [START_REF] Kavvadas | A constitutive model for clays based on non-associated anisotropic elasto-plasticity[END_REF], and MIT-E3 [START_REF] Whittle | Formulation of mit-e3 constitutive model for overconsolidated clays[END_REF]).

Similarly, [START_REF] Anandarajah | Bounding surface plasticity, iii: Application to anisotropic cohesive soils[END_REF] and [START_REF] Liang | Anisotropic plasticity model for undrained cyclic behavior of clays. ii: Verification[END_REF] introduced an induced stress invariant, based on the fabric tensor, to replace the corresponding stress invariant in the yield surface of the Modified Cam-Clay model in order to capture soil anisotropy. [START_REF] Banerjee | A plasticity model for the mechanical behaviour of anisotropically consolidated clay[END_REF] proposed an anisotropic model based on the fact that the yield surface is centred close to the K 0 -line. The yield function was obtained from experimental data without additional theoretical considerations and two hardening rules were considered, an isotropic and a kinematic one. The former was the same as the MCC, whereas the later was given as an empirical function accounting for anisotropy.

Yield function based on experimental results

Banerjee

Critical state models including anisotropy

In the aforementioned models, the assumptions related to anisotropy were directly postulated, rather than being related to the fundamental theory of plasticity. This gap was covered by Dafalias (1986a), who introduced an anisotropic theory based on the critical state concept. Details of this theory are explained in the following section.

Anisotropic critical state plasticity (ACSP)

Dafalias (1986a) proposed a plastic work dissipation equation with the following form in the triaxial space:

Ẇ p ≡ p εp v + q εp d = p ( εp v ) 2 + (M εp d ) 2 + 2α εp v εp d 1/2 (3.9)
with α being a non-dimensional variable accounting for the effect of anisotropy. α may be positive or negative depending on the compression or extension character of the triaxial loading. For α = 0 the work dissipation equation of the MCC model is recovered (Equation 2.106). It is noteworthy that the positive definition of the quadratic term under the square root in Equation 3.9 requires that |M | > |α| (given that M is negative in extension). Dividing Equation 3.9 by p εp d provides the following dilatancy equation:

εp v εp d = M 2 -η 2 2 (η -α) (3.10)
Assuming an associated flow rule, the equation of the yield surface in a triaxial plane can be derived from Equation 3.10 as follows:

f ≡ (q -pα) 2 + M 2 -α 2 p (p -p m ) = 0 (3.11)
where p m is the preconsolidation pressure. In the triaxial plane, the yield surface is an inclined ellipse, Figure 3.14. As mentioned before, when α = 0, the yield surface of the isotropic MCC model is retrieved. The yield surface can expand/contract and simultaneously rotate to model the magnitude of yielding under plastic deformations. The expansion and/or contraction of the yield surface is given by the variation of the preconsolidation pressure, as in the isotropic MCC model:

ṗm = λ -κ 1 + e 0 εp v (3.12)
whereas the rotation of the yield surface is controlled by a Rotational Hardening Rule (RHR), describing the rate equation of the internal variable α. It must be properly defined to correctly reproduce the evolution of the anisotropy under different loading conditions. A review of various RHR is provided hereafter.

Review of rotational hardening rules (RHRs)

The majority of the rotational hardening rules proposed in the literature are expressed as a function of the deviatoric and/or volumetric components of the plastic strain tensor (Dafalias (1986a), [START_REF] Whittle | Evaluation of a constitutive model for overconsolidated clays[END_REF], [START_REF] Ling | Anisotropic elastoplastic bounding surface model for cohesive soils[END_REF], among others). Generally, the RHR takes the following form:

α = A α b v (η) -α | εp v | + B α b d (η) -α | εp d | (3.13)
where A and B are functions of the stress state and the rotational hardening variable, or simply given as constants and α b v and α b d are instantaneous equilibrium values of α at the current stress ratio η. For example, if B = 0, a loading along a stress path of constant η would lead to α asymptotically approaching to α b v ; in the case of A = 0, α tends asymptotically to α b d . In the general case, when A and B are not null, α tends toward α b v for stress paths where plastic volumetric strains dominate (typically when η approaches zero) and towards α b d when shear plastic strains dominate (typically when approaching the critical state). Equation 3.13 does not represent all RHRs. Other contributions can be cited, such as the RHR proposed by [START_REF] Yang | From compression behavior to plastic anisotropy of reconstituted soft soils[END_REF] who considered a quadratic norm of deviatoric and volumetric strains and series of RHRs proposed by Dafalias [START_REF] Dafalias | Saniclay: simple anisotropic clay plasticity model[END_REF], [START_REF] Dafalias | Anatomy of rotational hardening in clay plasticity[END_REF], [START_REF] Dafalias | Rotational hardening with and without anisotropic fabric at critical state[END_REF]) who considered the plastic multiplicator ( λ) instead of components of the plastic strain tensor. Other models add an explicit dependency of the Lode Angle, precisely in the terms α v b and/or α d b , such as the model of [START_REF] Hashiguchi | Elastoplastic constitutive equation of soils with the subloading surface and the rotational hardening[END_REF], and [START_REF] Kobayashi | Numerical interpretation of a shape of yield surface obtained from stress probe tests[END_REF]. RHR can also consider the accumulation of strain due to cyclic loading, such as the model of [START_REF] Hueckel | Modeling of elastic anisotropy due to onedimensional plastic consolidation of clays[END_REF]. Some of these contributions and others propositions are summarized in Table 3.2. The list is not exhaustive, but includes the principal contributions based on the critical state concept.

It is noteworthy that the formulation of the RHR is not trivial, not simply a question of choice or accuracy with experimental results. The RHR must satisfy specific conditions. For example, the RHR must provide appropriate values of α under various loading conditions, such as loadings up to the critical state and should avoid to generate excessive rotation of the yield surface, since the work dissipation and the yield surface equations (Equation 3.11) imply |α| > |M | for real solutions to exist. The different conditions for the formula-Table 3.2: List of some well-known rotational hardening rules

RHR Constant Reference α = C p p 0 η x -α εp v C, x Dafalias (1986a) α = bp p 0 k -|α -ξ| k (η -α) εp v b,k,ξ Whittle (1993) α = b [η exp (-bε acc )] εp d b
Hueckel and Tutumluer (1994)

α = b |η -α| η -α |η -α| M (θ) -α εp d b Hashiguchi
and Chen (1998)

α = c β s + |β s | η -2 |β s | α 2 |η -α| εp d (β s = m e |η -α| (β 0 -tα)) c, m e ,t,β 0 Kobayashi et al. (2003) α = ω 3 4 η -α εp v + ω d 1 3 η -α | εp d | ω,ω d Wheeler et al. (2003b) α = b r (b l M -α) η |η| εp d b l , b r Zhang et al. (2007) α = (β s + |β s |) η -2 |β s | α εp d -ω d (2 |β s | α) εp v (β s = ω (|M -α| -η)) ω,ω d Yin et al. ( 2011
) α = λ C (α b -α) α b =        η αc M + m 1 -|η| M n if |η| ≤ |M | η αc M exp -µ |η| M -1 if |η| ≥ |M | C, α c ,m,n,µ Dafalias and Taiebat (2014) α = ω (α b -α) εp r where    α b = η + M 2 2ξ 1 -η M 2 ln 1 -η M ; εp r = εp v 2 + εp d 2 /2 ω,ξ Yang et al. ( 2014
)
tion of RHR are explained herein, based on three representative RHRs from the Table 3.2 .

The first RHR is proposed by Dafalias (1986a) and captures the correct stress ratio under K 0 consolidation. It reads:

α = A | εp v | (α b -α) ; α b = η x (3.14)
where A is a model parameter controlling the rate of α and x is a constant defining the equilibrium value of α under fixed stress-ratio loading.

The second RHR, proposed by Wheeler et al. (2003b), for the yield and plastic potential surfaces of equation 3.11), reads:

α = A α b v -α εp v + ω α b d -α | εp d | ; α b v = 3 4 η and α b d = 1 3 η (3.15)
where A and ω are model constants.

The third RH rule is the proposition of [START_REF] Dafalias | Rotational hardening with and without anisotropic fabric at critical state[END_REF] that separates the rate equation on the wet and dry sides (for η < M and η > M , respectively):

α = λ A (α b -α) ; α b =        η αc M + m 1 -|η| M n if |η| ≤ |M | η αc M exp -µ |η| M -1 if |η| ≥ |M | (3.16)
where λ is the rate of plastic multiplier, α c the value of the rotational hardening variable at the critical state, given as model parameter and A,m, n, and µ are constants.

In the following, advantages and limitations of the aforementioned RHRs are studied providing modifications and comparisons.

Uniqueness of CSL

The uniqueness (or lack of) critical state within the framework of rotational hardening was discussed by several authors (e.g., [START_REF] Karstunen | Discussion of "finite strain, anisotropic modified cam clay model with plastic spin. i: Theory" by george z. voyiadjis and chung r. song[END_REF], Wheeler et al. (2003b), and [START_REF] Dafalias | Anatomy of rotational hardening in clay plasticity[END_REF]). Wheeler and his co-workers highlighted the importance of shear strain rate in the evolution equation of α, and Dafalias discussed in details the question of the critical state value of α and the uniqueness of the CSL in ln p -e space. Conclusions and arguments presented by the afore-mentioned researchers are summarized and analysed hereinafter.

According to [START_REF] Dafalias | Anatomy of rotational hardening in clay plasticity[END_REF], the critical state is reached when the stress ratio η = M one has εp v = 0 (precisely ∂g ∂p = 0). Thence, the rate of p m is zero and the variable α reaches a critical state value α c and does not vary further. However, it must be asked how these conditions guarantee the uniqueness of the critical state in the ln p -e plane, independently of the shearing mode (triaxial compression, extension or anything in between) and loading history. This question is related to the uniqueness of the mean effective stress at the critical state, which defines the CSL in the ln p -e plane.

Let's examine this condition for the above-mentioned RHRs. It is clear that Equation 3.14 does not yield a unique CSL because of the term εp v . When the stress ratio reaches the state η = M , the volumetric plastic strain freezes the evolution of both α and p m (due to the zeroing of ∂g ∂p ). Such blocking of the evolution of α induces various values of the mean effective stress and therefore different CSLs in ln p -e space, as depicted in Figure 3.15 (a).

On the other hand, if the RHR is function of the shear plastic strain rate (Equation 3.15), or simply if the term ∂g ∂p is erased ( Equation3.16, the evolution of α is not stopped when εp v = 0. At the state η = M , α continues evolving until it reaches the equilibrium value in accordance with its rate equation. The p c value can be determined from the yield surface equation (by replacing q c = M p c ), which leads to:

p c = p 0 2 1 + α c M (3.17)
Thence, the value of the rotational hardening variable at the critical state α c is determined from the RHR. It follows from Equation 3.17 that the location of the CSL in the ln p -e plane depends on the ratio αc /M. The uniqueness of CSL is therefore related to the ratio αc /M that must remain the same for any loading process at any mean effective stress and shearing mode. In this respect, one must account for the fact that M , and eventually α c , depend on the Lode angle in multi-axial stress space. α c corresponds to the equilibrium value of α under constant stress ratio at the critical state. It is obtained by substituting η = M in the term α b of the RHR. Henceforth, the value α c in accordance with the chosen RHR is as follows:

• α c can not be defined for Equation 3.14

• α c = M /3, for Equation 3.15

• α c is given as model constant for Equation 3.16 Equation 3.15 provides therefore a unique CSL, since the ratio αc /M is assumed constant and does not depend on the shearing mode. α c is given as a model parameter in Equation 3.16. Thus, in order to have a non-dependency of αc /M on the shearing mode, α c should be chosen to have the same dependency on the Lode angle as M , so that the ratio αc /M provides a unique value of the mean effective stress at critical state. Otherwise, the condition of a unique CSL is violated.

Excessive rotation of the yield surface

The RHR must prevent excessive rotation of the yield and plastic potential surfaces, mainly because the work dissipation and yield surface equations require that |α| < |M | for real solutions to exist.

In order to test this condition, loading at constant stress ratio η > M should be considered. This condition leads to the following equilibrium values for the chosen RHR:

• α b = η /x for Equation 3.14;

• α b = η /3 for Equation 3.15;

• α b = α c for Equation 3.16.

It is clear from Equations 3.14 and 3.15 that the equilibrium value of α can exceed M for an appropriately high value of η. Thus, the inequality α < M can be. On the contrary, because of the negative exponential term in Equation 3.16, no excessive rotation occurs when |η| > M . Equation 3.16 is particularly interesting as it represents a common remedy for all RHRs regarding the pending question as to what a rotation hardening might be for |η| exceeding |M |. Following [START_REF] Dafalias | Rotational hardening with and without anisotropic fabric at critical state[END_REF] proposition, an exponential decay of α b can be applied on a RHR when the stress ratio exceeds the critical state slope M . In addition, one can write the RHR as a single equation by using the Macauley brackets .. . Despite the advantages of the Dafalias and Taiebat (2014)'s RHR formulation (Equation 3.16), its use (calibration) is not trivial and follows a complicated procedure as explained in [START_REF] Dafalias | Anatomy of rotational hardening in clay plasticity[END_REF]. Furthermore, the RHR needs the definition of five constants, compared to the Equation 3.15, which uses only two . For this reason, the RHR presented by the Equation 3.15 will be corrected in the following chapter to avoid excessive rotation of the yield surface. Most of the constitutive models have been developed based on the results of careful and comprehensive studies on remoulded and reconstituted natural soils or artificial materials such as kaolinite or illite (e.g., Taiebat et al. (2010a)). These studies have the advantage of repeatability and sampling uniformity. However, engineering practice typically encounters natural soils, which often behave differently, mainly due to the state of their structure. For instance, insitu soils develop tensile strength, and since their structure is often fragile they can develop unexpected instabilities (collapse). A key difference is that many natural soils present an undisturbed shear strength higher than a remoulded strength [START_REF] Skempton | The sensitivity of clays[END_REF], [START_REF] Nova | A constitutive model for bonded geomaterials subject to mechanical and/or chemical degradation[END_REF].

Sensitivity

The loss of strength that results from the disturbance of natural clays is generally quantified by the sensitivity of the initial structure. Terzaghi (1936b) defined sensitivity as the ratio of the intact to remoulded strength at the same water content. Sensitivity is generally considered as the parameter embodying the differences between natural and remoulded clays structures Cotecchia and Chandler (1998). [START_REF] Schmertmann | Swell sensitivity[END_REF] defined swelling sensitivity as the ratio of the swelling indices of intact and remoulded samples, respectively. Whereas, [START_REF] Burland | A laboratory study of the strength of four stiff clays[END_REF] suggested to use the ratio of the normalized strength at the critical state as an indicator of soil sensitivity. [START_REF] Burland | A laboratory study of the strength of four stiff clays[END_REF] highlighted the role of cohesion as a significant parameter of bonding and suggested the ratio of the cohesion of intact and reconstituted (unstructured) samples to quantify the degradation of the structure.

Sensitivity varies from about 1 for heavily over-consolidated clays to over 100 for the so-called extra-sensitive or quick clays. Figure 3.16 illustrates the sensitivity of some clays as a function of their liquid limit and undrained shear strength. Clays of low sensitivity (1-2) are very rare (e.g. the insensitive clay of London [START_REF] Gens | Stress-strain and strength characteristics of a low plasticity clay[END_REF]), sensitivities of 2-4 are very common among normally consolidated clayey soils, sensitivities of 4 to 8 are also frequently encountered, while the extra-sensitive clays are frequent in parts of Scandinavia, Canada [START_REF] Locat | The compressibility and sensitivity of an artificially sedimented clay soil: The grande-baleine marine clay, québec, canada[END_REF], Japan and Alaska [START_REF] Ohtsubo | Marine quick clays from ariake bay area, japan[END_REF]. [START_REF] Leroueil | Discussion on: Correlations between index tests and the properties of remoulded clays[END_REF] Many factors can be taken into account when studying the sensitivity of natural soils, such as temperature, pH, and salinity [START_REF] Torrance | Towards a general model of quick clay development[END_REF]. The first work on the role of chemical factors was undertaken in Scandinavia [START_REF] Rosenqvist | Considerations on the sensitivity of norwegian quickclays[END_REF], when Rosenqvist observed that the Norwegian clays, which were subjected to disastrous landslides, all contained pore water of low salinity despite their marine origins. This finding prompted Rosenqvist to propose leaching as responsible for the development of the marine quick clays. A series of investigations carried out by [START_REF] Bjerrum | Geotechnical properties of norwegian marine clays[END_REF] proved moreover that a relationship exists between the Atterberg limits and the pore-water salinity. The author observed that the liquid limit of the Norwegian clays decreased with the porewater salinity until a low value for salinities less than 5g/l. A similar effect was observed for the plastic limit. As a result, the sensitivity of marine clays was proved to be inversely proportional to the pore-water salinity. [START_REF] Leroueil | The general and congruent effects of structure in natural soils and weak rocks[END_REF] demonstrated through an extensive experimental study that the chemical and geological characteristics of Scandinavian clays are common in natural geological materials. The sensitivity of natural soils have been recognised in soft and stiff clays, granular soils, as well as weak and weathered rocks. [START_REF] Leroueil | The general and congruent effects of structure in natural soils and weak rocks[END_REF] highlighted, through the same study, the need to consider the behaviour of structured materials in relation with the behaviour of the equivalent reconstituted ones.

The following section briefly discusses some of the most common experi-mental observations on the mechanical behaviour of intact soils.

Experimental observations on intact clays

As [START_REF] Leroueil | The general and congruent effects of structure in natural soils and weak rocks[END_REF] had shown, similar patterns of behaviour are observed in materials of different geological origins. Some of these patterns, are presented in Figures 3.17 The figures reveal the most common behaviours of intact natural soils. Three different effects of the soil structure are discussed hereafter:

( 1 ¢ 3 )/2£ p (-) + (-) (¤ 1 ¥ 3 )/2¦ p (-) - 0 0.4 0.2 0.6 0.8 0 0.2 0.4 ( § 1 ¨ 3 )/2© p (-) - + 0.6 1.0 Intact Remoulded ϕ' =30°( 1 3 )/2 p (-) Intact Remoulded ϕ' =23°Y ield points
• Effect on the Isotropic Consolidation Line (ICL) Evidence of structure is usually provided by comparison of oedometer or isotropic consolidation tests on natural and remoulded soils. As shown in Figure 3.17 (a), the consolidation line for the intact natural soil has a higher void ratio than for the remoulded soil (from which the structure has been partially removed). This behaviour have been observed in numerous studies, e.g. [START_REF] Casagrande | The structure of clay and its importance in foundation engineering[END_REF], [START_REF] Quigley | The fabric of anisotropically consolidated sensitive marine clay[END_REF], [START_REF] Mesri | Composition and compressibility of typical samples of mexico city clay[END_REF], [START_REF] Locat | The compressibility and sensitivity of an artificially sedimented clay soil: The grande-baleine marine clay, québec, canada[END_REF].

• Effect on the yield surface

The structure of an intact soil produces an increase of the preconsolidation pressure, the limit state curve and the entire yield curve, as shown in Figure 3.17 (b), compared to a remoulded soil. The yield surface of some intact natural clays is clearly determined, even in isotropic consolidation (e.g., [START_REF] Gens | Conceptual bases for a constitutive model for bonded soils and weak rocks[END_REF], [START_REF] Leroueil | The general and congruent effects of structure in natural soils and weak rocks[END_REF]).

• Effect on the tensile strength

The yielding behaviour of intact soils usually presents some tensile strength, as depicted in Figure 3.17 (d). The tensile strength of natural and undisturbed soils was also highlighted by [START_REF] Saxena | Static properties of lightly cemented sand[END_REF] and [START_REF] Clough | Cemented sands under static loading[END_REF]. The authors demonstrated that bonding of intact and cemented soils increases the peak strength, the initial stiffness and provides some tensile strength.

On the whole, the behaviour of intact and natural soils is related to their state of stucture, which is strongly affected by the decrease of bonding during yielding. The presence of bonding allows the intact natural soil to exist at higher strengh and void ratio than would otherwise be possible. Degradation of bonding, as the result of yielding, leads to a decrease of this apparent strength. This behaviour is known as desctructuration and is typically introduced into an elastoplastic model by assuming that the size of the intact soil is proportional to the amount of destructuration, which follows a softening hardening, (e.g. [START_REF] Gens | Conceptual bases for a constitutive model for bonded soils and weak rocks[END_REF], [START_REF] Rouainia | A kinematic hardening constitutive model for natural clays with loss of structure[END_REF], and Taiebat et al. (2010b)). Such an isotropic softening could be combined with the isotropic hardening due to consolidation. Modelling destructuration in natural soils is explained hereafter in the framework of elastoplasticity.

Modelling destructuration

Gens and Nova (1993) presented a general framework for incorporating destructuration within elastoplastic constitutive models. In addition to the yield surface of the intact natural material, a notional intrinsic yield surface is introduced to represent the yield surface of the equivalent reconstituted material, at the same void ratio and with the same fabric. The intrinsic yield surface is assumed to be of the same shape as the yield surface of the natural soil, but smaller in size. The size of the intrinsic yield surface is characterized by the parameter p mi , which is the preconsolidation pressure of the equivalent reconstituted sample. p mi is related to p m , the size of the real yield surface, by the mean of the parameter χ, which defines the amount of destructuration. χ may be considered as an indicator of the sensitivity of intact natural soils, according to [START_REF] Schmertmann | Swell sensitivity[END_REF]:

p m = (1 + χ) p mi (3.18)
The intrinsic yield surface is allowed to expand and/or contract in accordance with the isotropic hardening law for reconstituted soils:

ṗmi = λ i -κ 1 + e 0 εp v (3.19)
where λ i is the slope of the isotropic consolidation line of reconstituted sample, which is considered as an intrinsic parameter.

The destructuration is related to the plastic strain increment, and given by a well-defined softening law. In the following, a brief review of rules controlling the destructuration is presented.

Destructuration law

In recent years, there have been important developments in formulating appropriate softening laws for bonded soils. Table 3.4 summarizes some well-known destructuration laws applied within the framework of elastoplasticity. The various models differ in the precise form of the destructuration law and/or the form of the baseline model employed for the intrinsic material. Generally, the destructuration is assumed to start at the beginning of yielding until complete loss of structure. Thus, the majority of destructuration laws propose a dependency on volumetric and deviatoric plastic strains. The former leads to the damage of bonding under small values of stress-ratio, while the later gives rise to complete loss of structure for loading up to the critical state (e.g., [START_REF] Gens | Conceptual bases for a constitutive model for bonded soils and weak rocks[END_REF], [START_REF] Rouainia | A kinematic hardening constitutive model for natural clays with loss of structure[END_REF], Taiebat et al. (2010a), [START_REF] Hong | A two-surface plasticity model for stiff clay[END_REF] and [START_REF] Yang | From compression behavior to plastic anisotropy of reconstituted soft soils[END_REF]). On the other hand, some models assumed that the destructuration is only function of deviatoric deformation, which is the case of [START_REF] Oka | An elasto-viscoplastic model for clay considering destructuralization and prediction of compaction bands[END_REF] and [START_REF] Baudet | A constitutive model for structured clays[END_REF]. Furthermore, [START_REF] Baudet | A constitutive model for structured clays[END_REF] have tried to capture the possibility of having a residual state of structure. Instead of a complete loss of bonding, Baudet assumed that the destructuration parameter χ is able to decrease with plastic straining until reaching a final state χ f , that defines the residual state of structure.

Among these destructuration laws, the softening rule proposed by [START_REF] Gens | Conceptual bases for a constitutive model for bonded soils and weak rocks[END_REF] is selected in the present research, owing to its analytical simplicity and its success in capturing the destructuration of various clays within the framework of the critical state plasticity (e.g., Wheeler et al. (2003b), Karstunen and Koskinen (2004a)). The expression of the softening rule is given in Table 3.4, and implies two constants, ξ and ξ d .

Tensile strength of structured material

The initial structure of the intact natural clays is known also to reveal tensile strength and cohesion to the material. This cohesion is also affected by the destructuration, so it decreases with plastic straining and causes an additional softened behaviour to the structured material. Therefore, the destucturation law can be combined with an additional softening law that controls the degradation of tensile strength with damage of bonding. The combination of bonding and cohesion degradation is known as destructuration mechanism Taiebat et al. (2010b), and Karstunen and Koskinen (2004b). Degradation of the cohesion as the result of bonding damage is usually given as a decreasing Table 3.4: List of some well-known destructuration laws

Destructuration law

Constant Reference

χ = -χξ (| εp v | + ξ d | εp d |) ξ, ξ d Gens and Nova (1993) χ = -ξ (χ -1) ψ εp a εp a = (1 -A) ( εp v ) + A ( εp d ) ξ, ψ ,A Rouainia and Muir Wood (2000) χ = -χξ εp d ξ,
Oka and Kimoto ( 2005)

χ = -ξ (χ -χ f ) εp d ξ, χ f Baudet and Stallebrass (2004) χ = -ξ i (ξ -1) εp a εp a = (1 -A) ( εp v ) 2 + A ( εp d ) 2 ξ i , A Taiebat et al. (2010a) χ = -ξ d εp a where εp a = ( εp v ) 2 + εp d 2 /2 ξ d Yang et al. ( 2014
)
function of shear strain increments [START_REF] Gens | Conceptual bases for a constitutive model for bonded soils and weak rocks[END_REF], [START_REF] Nova | A constitutive model for bonded geomaterials subject to mechanical and/or chemical degradation[END_REF], and Wheeler et al. (2003a). Generally, it takes the following form [START_REF] Nova | A constitutive model for bonded geomaterials subject to mechanical and/or chemical degradation[END_REF]:

ṗt = -ξ t p t | εp d | (3.20)
where p t is the tensile strength of the structured material, and ξ t is a constant that controls its degradation.

Anisotropy and destructuration

The question that finally arises is the relation between anisotropy and bonding. Since these two properties have been usually treated separately, anisotropy was neglected in several models for structured soils, while for some others anisotropy was solely supposed related to the structure, i.e. as soon as the destructuration is complete, anisotropy disappears (e.g., [START_REF] Kobayashi | Numerical interpretation of a shape of yield surface obtained from stress probe tests[END_REF], [START_REF] Yin | Modeling time-dependent behavior of soft sensitive clay[END_REF]). However, this assumption does not fit with the experimental investigations on soft clays, which show that anisotropy of reconstituted samples can be as great as the anisotropy of the corresponding natural and structured samples [START_REF] Leroueil | The general and congruent effects of structure in natural soils and weak rocks[END_REF], Taiebat et al. (2010b).

Hence, anisotropy should be considered as an intrinsic feature related entirely to the fabric, i.e., the real and intrinsic yield surfaces are subjected to a unique rotational hardening variable α.

Conclusion

This chapter deals with the principal features of natural clays behaviour under monotonic loading and presented several experimental and numerical issues. The clay behaviour depends on the nature of their structure, a combination of the fabric, the component particles arrangement and bonding.

The first section discusses the behaviour of overconsolidated clays, which are characterized by two principal features: dilatancy and high non-linearity. The classical critical state models capture successfully the first, but are unable to reproduce the later. In this context, appropriate constitutive laws for overconsolidated clays have been reviewed, in particular bounding surface plasticity models within the framework of the critical state theory have been presented in detail.

In the second section, the effect of anisotropy on the clay behaviour has been discussed and related to the undrained and drained behaviour of cohesive soils. In the framework of the critical state theory, anisotropy is usually modelled via a Rotational Hardening Rule (RHR). The anisotropic critical state plasticity has been described and a short review on several RHR has been presented

The chapter ends with the behaviour of intact and sensitive clays, which loose rapidly their strength when disturbed. Several experiments have been presented to reveal the effect of structure on the clay behaviour. The destructuration mechanism, which gathers softening rules that describe damage of bonding, has been defined and some-well known softening rules from the literature have been presented.

Introduction

Since its publication, the Modified Cam-Clay (MCC) model has given birth to numerous constitutive models aiming to capture different features of clays (see Chapter 3). For instance, Wheeler et al. (2003b) proposed a constitutive law, named SCLAY-1 for reconstituted and soft clays by introducing a rotational hardening to the MCC in order to predict the initial and induced anisotropies. Later on, Wheeler et al. (2003a) proposed a versatile extension of SCLAY-1 model, termed SCLAY-1S law, which incorporates two distinct mechanisms for modelling anisotropy and destructuration. The former is the rotational hardening rule of SCLAY-1, whereas the later is governed by two distinct softening rules that control the state of the clay structure.

Although the SCLAY-1S model is able to predict anisotropy as well as destructuration for normally and lightly overconsolidated natural clays, the model does not capture the response of overconsolidated (OC) clays. This deficiency arises from the classical theory of elastoplasticity, in which plastic strains are avoided for stress states within the yield surface. Consequently, the SCLAY-1S model, as all classical elastoplastic laws, overestimates the elastic strains as well as the peak strength for OC clays.

The bounding surface, proposed initially by Dafalias (1986b), came to overlap this gap, by allowing plastic deformations for stress states within the bounding surface, at an early state of straining. The formulation of bounding surface plasticity belongs to the framework of elastoplasticity, and includes a dependency of the plastic modulus to the distance between the stress point and a so-called image point lying on the bounding surface.

This chapter presents an extension of SCLAY-1S model. The proposed model, named SCLAY-1SB, introduces the concept of bounding surface plasticity to capture the behaviour of overconsolidated clays, as well as the anisotropy and destructuration features of natural and sensitive clays.

The chapter opens with the formulation of the SCLAY-1SB, where equations of the initial model and bounding surface plasticity are recalled. The constitutive equations are stated with respect to triaxial and multi-stress conditions. The second part describes the calibration of constants and parameters of the model. Finally, experiments on three different natural clays are simulated to validate various aspects of the constitutive law. The tests include oedometer as well as drained and undrained triaxial tests.

The SCLAY-1SB model

In this chapter, stress quantities are considered as effective stresses and the compression positive convention is used for both stress and strain measures.

Triaxial formulation

The mean effective stress p = (σa + 2σr) /3, and deviatoric stress q = (σ a -σ r ), are used as stress measures, while the volumetric strain, ε v = ε a + 2ε r , and deviatoric strain ε q = 2 /3 (ε a -ε r ), are used as strain measures, where ε a and ε r being the axial and radial strains, respectively.

Elastic properties

The model adopts the additive decomposition of the total strain rate, into elastic (suscript e ) and plastic (subscript p ) parts. It reads:

ε = εe + εp (4.1)
The elasticity is supposed to be hypoelastic (non linear elasticity) and isotropic, for which the bulk and shear moduli are expressed as:

K = 1 + e κ p, G = 3K (1 -2ν) 2 (1 + ν) (4.2)
where e is a void ratio of the material, ν the Poisson's ratio, and κ the rebound line in the e -ln p plane.

Bounding surface: Intrinsic and Structured state

Following [START_REF] Gens | Conceptual bases for a constitutive model for bonded soils and weak rocks[END_REF], the concept of the intrinsic yield surface is adopted in the SCLAY-1S model. Accordingly, the intrinsic surface takes the same form as the yield surface of the baseline model (SCLAY-1S). As illustrated in Figure 4.1 (a), the yield surface of the SCLAY-1S is a rotated and distorted ellipse in triaxial stress space. The rotation and distortion is controlled by the internal variable α, which represents the degree of plastic anisotropy, while the size of the intrinsic surface is governed by the internal variable p mi . The yield surface shares the same shape as the intrinsic surface but is enlarged both towards the compressive and tensile regions of the stress space, through the internal variable p m and p t , respectively, see figure 4.1 (a).

The former corresponds to the preconsolidation pressure of the natural and structured clay, while the later is the tensile strength. The same concept is employed hereafter and combined with the bounding surface plasticity for the SCLAY-1SB model (see Figure 4.1 (b)). The bounding surface F = 0 of the structured soil is expressed as:

F ≡ (q -(p + p t ) α) 2 + M 2 -α 2 (p + p t ) (p -p m ) = 0 (4.3)
where a bar over stress quantities indicates an image point that lies on the bounding surface. Image and current stresses are properly linked by a well defined mapping rule. p m and p t are related to the amount of inter-particle bonds. p m is linked to p mi by the mean of the internal variable χ (see Equation 3.18), which defines the degree of destructuration.

Radial mapping rule

The radial mapping rule proposed by Dafalias (1986b) is employed herein to define the image stress. Using (p c , q c ) as the projection center, (p, q) is obtained by radially projecting the current stress (p, q) into the bounding surface. It yields:

p = p c + b (p -p c ) q = q c + b (q -q c ) (4.4)
where the variable b varies from 1 to ∞, with the two limits being attained when the current stress coincides with either the image stress (b = 1) or the projection center (b = ∞).

The projection center is assumed lying on the origin of the bounding surface, at any stress state (see Figure 4.1 (b)). In the presence of a tensile strength it reads:

q c = 0 ; p c = p t (4.5) qc = 0 ; ṗc = ṗt (4.6)
It is noteworthy that this assumption has sense only in the case of virgin loading, and is not suited for complex loading when cycles of unloading/reloading are considered. In this cases, elastic strains are often overestimated. This point will be reviewed in Chapter 6, when cyclic features of clay will be discussed.

Plastic modulus

The bounding surface theory requires the definition of a shape function, which relates the plastic moduli at image and current stresses, respectively. The former is determined by enforcing the consistency condition on the image stress, while the later is analytically expressed as follows:

K p = K p + H (4.7)
The shape function proposed by [START_REF] Yin | Modeling time-dependent behavior of soft sensitive clay[END_REF] is employed herein:

H = H 0 . M 4 1 + e λ -κ (bp m ) 3 . 1 - 1 b (4.8)
H 0 is a constant controlling the steepness of the response in the q -ε q plane.

Hardening rules

1. Isotropic hardening:

The hardening rule of the intrinsic preconsolidation pressure p mi is identical to the hardening rule of the MCC. It yields:

ṗmi = λ -κ 1 + e p mi εp v (4.9)
where κ and λ i are the slope of the rebound and compression isotropic consolidation line of a reconstituted sample, and e is the void ratio of the material.

Destructuration mechanism:

The destructuration mechanism is governed by the softening rules of [START_REF] Gens | Conceptual bases for a constitutive model for bonded soils and weak rocks[END_REF], which control the degradation of the tensile strength and the degree of destructuration of sensitive clays:

χ = -χξ | εp v | + ξ d εp q (4.10) ṗt = -ξ t p t εp q (4.11)
ξ, ξ d are constants controlling the degradation of degree of destructuration χ, and ξ t is the constant of the softening rule of the tensile strength.

Rotational hardening:

The reference model (SCLAY-1S) employs the RHR of Wheeler et al. (2003a), which is expressed by Equation 3.15 and has two constants A and ω d . Equation 3.15 has been studied together with other representative RHRs from the literature in Chapter 3. It has been pointed out that the RHR of Wheeler et al. (2003a) has some strong and weak points, summarized in the following:

• Its formulation is relatively simple and uses few constants, compared to others RHRs.

• It was calibrated and validated by many authors in the past. It has been implemented in various finite element codes and used in many numerical studies.

• It captures the uniqueness of the CSL in ln p -e plane. An enhancement of Equation 3.15 is presented herein. The proposed modification explores the formulation employed by [START_REF] Dafalias | Rotational hardening with and without anisotropic fabric at critical state[END_REF], which uses a negative exponential term in the RHR for stress ratio exceeding the critical state. Thus, no excessive rotation occurs during plastic straining. This formulation is particularly interesting, as it represents a common remedy for all RHRs regarding the formulation of the RHR for |η| exceeding |M |. Following [START_REF] Dafalias | Rotational hardening with and without anisotropic fabric at critical state[END_REF] proposition, an exponential decay of α b can be applied on a RHR when the stress ratio exceeds the critical state slope M . In addition, the RH rule can be written as a single equation by using the Macauley brackets .. . It yields:

α = A 4 3 η -α εp v + ω d 1 3 η η-M |η-M | M η exp µ 1 -η M -α εp q (4.
12) where µ > 1 is an additional constant, which controls the decay of the variable α when |η| > |M |.

For stress state ratio η < M , Equation 4.30 is identical to the RHR of Wheeler et al. (2003b). When the stress ratio exceeds the critical state slope, the term inside the Macauley brackets becomes positive and the negative exponential term is activated to avoid excessive rotation of the yield surface.

It is important to point out that the proposed correction guarantees the continuity of the RHR at the stress ratio η = M . At this state, α tends towards the equilibrium value α

b (η = M + ) = α b (η = M -) = M /3.

Generalized formulation

The three-dimensional formulation of the SCLAY-1SB model is provided in this section. Second-order tensors are hereafter denoted by boldface symbols.

Triaxial stress quantities are generalized to multiaxial space using a second order tensor. The stress tensor σ is decomposed into its hydrostatic and deviatoric components, pδ and s, as follows:

p = 1 3 tr (σ) (4.13) s = σ -pδ (4.14)
where tr is the trace operator and δ is an identity tensor.

Similarly, triaxial strain quantities are generalized to multiaxial space. The strain tensor ε is decomposed into its volumetric ε p v δ/3 and deviatoric parts ε d :

ε v = tr (ε) (4.15) ε d = ε - 1 3 ε v δ (4.16)
A systematic generalization of the triaxial constitutive relation is based on the observation that under triaxial conditions, the following holds true:

x = [( 3 /2) X : X] 1 /2 (4.17)
where X is any deviatoric stress or strain tensor and x is its triaxial counterpart, the symbol : denotes the trace of the product of adjacent tensors. On the basis of Equation 4.17, it yields:

q = [( 3 /2) s : s] 1 /2 (4.18) ε q = [( 3 /2) ε d : ε d ] 1 /2 (4.19) α = [( 3 /2) α : α] 1 /2 (4.20) (4.21)
where s and ε d are deviatoric stress and strain tensor. α is a deviatoric tensor controlling the RHR.

Elastic relations

The multiaxial generalization of the elastic relationship is given by:

ε = ṡ 2G + ṗ 3K δ (4.22)
where the bulk and shear moduli K and G are given by Equation 4.2.

90CHAPTER 4. THE SCLAY-1SB MODEL FOR MONOTONIC LOADING

Bounding surface

The equation of the bounding surface in the multiaxial space takes the following form:

F ≡ 3 2 (s -(p + p t ) α) : (s -(p + p t ) α) + M 2 -3 2 α : α (p + p t ) (p -p m ) = 0 (4.
23) M is considered to be Lode angle dependent. It is interpolated from its values M c in compression and M e in extension as proposed by [START_REF] Van-Eekelen | Isotropic yield surfaces in three dimensions for use in soil mechanics[END_REF] :

M = a (1 + b sin 3β) n (4.24)
Parameter n controls the yield surface convexity in the deviatoric plane, and is generally chosen equal to -0.229 (default value) to ensure the convexity of the yield surface (Figure 4.2). Coefficients a and b allow an independent choice of the compression and extension friction angles, ϕ c , and ϕ e :

a = M c (1 + b) n ; b = m 1 n -1 m 1 n + 1 , with m = M c M e
where the slopes of the critical state line in compression M c and in extension M e read:

       M c = 6 sin ϕ c 3 -sin ϕ c M e =
6 sin ϕ e 3 + sin ϕ e As a result, the shape of the critical state surface in the deviatoric plane is not circular, and it is possible to define a friction angle in compression and extension. Furthermore, in the framework of anisotropic critical state plasticity, the Lode angle β is calculated from the deviatoric stress-ratio tensor (η -α) (Dafalias (1986a):

sin3β = - 3 √ 3 2 J 3 J 2 (4.26)
where J 2 and J 3 are respectively the second and the third invariant of the deviatoric tensor η -α; with η = 1 p+pt s.

Mapping rule

The radial mapping rule is generalized simply by replacing the deviatoric quantity q in Equation 4.4 by the deviatoric tensor s. It reads:

s = s c + b (s -s c ) (4.27) p = p c + b (p -p c ) (4.28)
Similarly, the rate equation of the ratio b is expressed in the multi-axial space as follows:

ḃ = - Kp K P -b ∂F ∂κ κ + Kp Kp ∂F ∂p c ṗc + ∂F ∂sc ṡc ∂F ∂p (p -p c ) + ∂F ∂s (s -s c ) b + Kp K P -b ∂F ∂b (4.29)

Rate evolution equation of the internal variables

The rate equations of the internal variables p mi , χ and p t in the mutliaxial space are identical to those in the triaxial plane, and given by the Equations 4.9, 4.10 and 4.11, respectively. On the other hand, the counterpart of the RH variable α is a stress-ratio tensor, and the generalization of its rate equation takes the following tensorial form:

α = A 4 3 η -α εp v + ω d 1 3 η η-M |η-M | M η exp µ 1 -η M -α εp q ( 4
.30) where η denotes the stress ratio tensor (η = 1 p+pt s), while η is its counterpart scalar value in the triaxial plane.

Calibration method

The SCLAY-1SB model employs 7 internal variables and requires definition of 12 constants when all the mechanisms are activated (see Tables 4.1 and4.2). Determination of these parameters is based entirely on standard laboratory tests such as oedometer and triaxial tests, as explained hereafter.

Poisson's ratio ν

The Poisson's ratio is calibrated using K 0 -consolidation test results (see Figure 4.3). The test is carried out in the triaxial cell by imposing lateral strains equal to zero (ε q = 0) and a constant strain rate ratio, i.e.:

= εv εq = εa + 2 εr 2 /3 ( εa -εr ) = 3 2 (4.31)
If the unloading path is considered elastic, Equation 4.31 gives:

q ṗ = 2G K . 1 = 2 (1 -2ν) (1 + ν) (4.32)
Hence, comparison with the slope of the experimental data leads to determination of the Poisson's ratio. For the Bothkennar clay, this method leads to ν = 0.20, as presented in 3.2 Parameters κ and λ i κ and λ i are defined as the slopes of the swelling and the normal compression lines in ln p-e plane, respectively. Their values can therefore be obtained from isotropic compression or one-dimensional consolidation (or oedometer) tests, using reconstituted samples. An example is presented in Figure 4.4, which shows data of an oedometer test on reconstituted samples of Bothkennar clay, plotted in log p -e plane. A simple conversion to the ln p -e plane must be carried out to estimate the appropriate values of κ and λ i . mental data are not available, the value of M e can be estimated from Equation 4.25, assuming that ϕ e = ϕ c . This assumption gives more accurate results than considering M e = M c , which overestimates the strength in extension.

Parameter H 0 for the bounding surface

Calibration of H 0 is based on stress and stress-strain paths of monotonic loading on overconsolidated samples.

Constants controlling the variable α

The calibration method to identify the different constants controlling the variable α is detailed hereafter:

• The initial value of α corresponds to the equilibrium value of α under K 0 -consolidation. As reported before, the K 0 -consolidation is carried out in a triaxial cell by imposing constant stress ratio and lateral strain equal to zero. These conditions lead to Equation 4.31.

On the other hand, K 0 consolidation initially induces simultaneous rotation and expansion/contraction of the yield surface, in which the stress state ratio η and the variable α tend to an equilibrium value, while the preconsolidation pressure evolves according to Equation 3.12 along the normal consolidation line (NCL). Through this process, the elastic volumetric strain rate is defined in terms of the slope κ of the swelling line as:

εv = λ (λ -κ) εp v (4.33)
By neglecting the shear elastic strain, i.e. εe εe v δ , the expression of the dilatancy is found:

εp v εp d = 1 - κ λ (4.34)
According to the work dissipation equation the dilatancy is also expressed by Equation 3.10, reported hereafter, as function of the variable α and the stress ratio:

εp v εp d = M 2 -η 2 2 (η -α)
Combining the two previous equations at K 0 conditions, i.e., η = η K 0 and α = α K 0 , gives:

α K0 = η 2 K0 + 3 (1 -κ /λ) η K0 -M 2 3 (1 -κ /λ) (4.35)
Equation 4.35, demonstrated by Dafalias (1986a), does not depend on the RHR by virtue of its derivation. It is also independent on the yield surface expression. It can be therefore also applied in case the yield surface is different from the plastic potential, i.e. non-associated flow rule.

• The constant ω d is estimated using the results of K 0 -consolidtion. According to Wheeler et al. (2003b), zeroing the RHR ( α = 0), in the case of K 0 loading leads to :

ω d = 3 (4M 2 -4η 2 k0 -3η k0 ) 8 (η 2 k0 -M 2 + 2η k0 ) (4.36)
• The constant A is estimated by trial and error using experimental results on drained or undrained triaxial tests. Several calibration methods based on analytical and optimization criteria have been proposed, see for example [START_REF] Yin | Modeling time-dependent behavior of soft sensitive clay[END_REF]) and [START_REF] Yin | Optimization techniques for identifying soil parameters in geotechnical engineering: Comparative study and enhancement[END_REF].

• The constant µ controls the rotation of the yield surface for stress ratio |η| > |M |. Its calibration can be therefore made using a triaxial test under constant stress ratio |η| > |M |.

Constants controlling the destructuration law

Gens and Nova (1993) suggested a procedure for calibrating χ 0 , the initial value of χ, based on isotropic consolidation or oedometer experiments on reconstituted and intact samples until complete lose of structure. χ 0 is determined measuring λ and λ i , the slopes of the compression line of the intact and reconstituted samples respectively, since λ /λ i = 1 + χ 0 . The determination of the constants ξ and ξ d requires an optimization procedure, involving model simulations of laboratory tests on intact samples. The calibration is first performed by selecting the value of ξ based on a triaxial test with small values of η (where the shear strains can be neglected, as in compression consolidation). Then, triaxial tests involving higher stress ratio can be used to select the most appropriate value of ξ d .

Constants controlling the tensile strength p t

Calibration of the initial value of p t and the constant ξ t is performed by trial and error on simulation results and triaxial tests experiments.

The best method to estimate p t0 is to consider an isotropic consolidation test (or triaxial test with small value of η), so that the tensile strength can be assumed constant since the plastic shear strains are neglected. The method provides comparison of the critical state line in stress plane for intact and reconstituted samples. The CSL of the reconstituted and intact samples are parallel in the stress plane. The CSL of the former passes through the origin of stress, while the later intersects the hydrostatic axis in the initial value of p t .

The calibration of ξ t implies experiments with a high value of stress ratio, taking thus into account the effect of shear plastic strains.

Validation with laboratory tests

The performance of the SCLAY-1SB to predict the behaviour of natural clays is examined in this section. For this purpose, three different natural clays are chosen to examine the various aspects of the model, that are, the dilatancy of the OC clays, the anisotropy and the destructuration. First, the bounding surface ability to capture the behaviour of OC clay is examined by simulating drained triaxial tests on isotropic consolidated and reconstituted samples of London clay. Second, drained triaxial tests under constant stress ratio are employed to capture the evolution of soil anisotropy. Simulations are based on reconstituted and normally consolidated samples of Otaniemi clay [START_REF] Karstunen | Plastic anisotropy of soft reconstituted clays[END_REF]. Therefore, only the RHR is activated. Lastly, intact and structured samples of the sensitive Bothkennar clay [START_REF] Smith | The yielding of bothkennar clay[END_REF], and [START_REF] Allman | Mechanical properties of reconstituted bothkennar soil[END_REF]) are chosen to examine the destructuration mechanism.

London clay

Results of four drained triaxial tests for reconstituted samples of London clay with overconsolidation ratios 1, 1.5, 2, 4, 10 are presented. Samples are isotropicaly consolidated at mean effective stress p m0 = 234kP a, and afterwards unloaded up to the value p initial , that is reported in Table 4.3 for each test. Further details on the samples preparation and testing procedure are given in [START_REF] Gens | Stress-strain and strength characteristics of a low plasticity clay[END_REF]. Material parameters used for simulation are presented in Table 4.4. Only the constants related to the bounding surface (H 0 ) and anisotropy (α 0 , A, and ω d ) are calibrated for the London Clay, all the other parameters are reported by the cited authors. The initial conditions of the samples are given in Table 4.3 and the results of test simulation are given in Figure 4.5.

Table 4.3: Initial conditions and tests data on London clay

Sample OCR e 0 α 0 CDIC 1 10 0.392 0 CDIC 2 4 0.398 0 CDIC 3 2 0.404 0 CDIC 4 1.5 0.407 0 CDIC 5 1.0 0.41 0 In Figure 4.5, the presented results in q -ε a plane show a good matching with the experimental curves for all overconsolidation ratios. Furthemore, the general form of the dilatancy curve (Figure 4.5(b)) is well predicted. The biggest deviations can be observed at the overconsolidation ratio 4, for which about 20% of volumetric strain is overestimated. One manner to improve this point is to consider a non-associated flow rule. 4.4 summarizes the material parameters of the SCLAY-1SB model, table 4.5 the loading data of the chosen tests, while further details on the experimental campaign can be found in [START_REF] Karstunen | Plastic anisotropy of soft reconstituted clays[END_REF].

Table 4.5: Simulated tests on Otaniemi clay, from [START_REF] Karstunen | Plastic anisotropy of soft reconstituted clays[END_REF] Test No e 0 (-) q 0 (kPa) p 0 (kPa) The tests are carried out under constant stress-ratio as illustrated in Figure 4.6. The tests are designed to produce considerable yield rotation, so that the effect of the rotational hardening can be investigated. In addition, reconstituted and normally consolidated samples are employed to deactivate the effect of the structure, and dilatancy of OC clays, that is χ 0 = 0, p t = 0, b = 1, K p = K p and the constant H 0 are of no use.

The tests start with an initial consolidation stage, where all samples are loaded along the same stress path at the stress ratio η K 0 , until the stress level p 0max , and unloaded following the stress path η K 0 up to a mean effective stress p (between 6kP a and 12kP a). Then, two loading/unloading cycles are carried out with two different stress ratios η 1 and η 2 .

p' q p' 0 q 0 p' 0max p' 1max p' 2max ¡ K0 1 ¢ 2 Figure 4
.6: Triaxial tests procedure on Otaniemi clay as described in [START_REF] Karstunen | Plastic anisotropy of soft reconstituted clays[END_REF] Simulation of the aforementioned tests is limited to the last two cycles of loading under η 1 and η 2 . Furthermore, simulations are also carried with the MCC model, in order to highlight the effect of anisotorpy on the clay response.

Results are presented in Figure 4.7 in terms of the mean effective stress versus volumetric strain and shear stress versus deviatoric strain. The SCLAY-1SB model provides much better results compared to the MCC model. In particular, the MCC model overestimates the deviatoric shear strains for the four simulations compared to the SCLAY-1SB model.

For the tests where loading begins with compression, i.e. CAE 3516R and CAE 3519R, simulations with the SCLAY-1SB model provide reasonable results during the first loading stage, while the volumetric and deviatoric shear strains are slightly overestimated during the second loading/unloading phase. On the contrary, simulations of the CAE 3819R and CAE 3820R tests which start with extension are less accurate regarding volumetric strains (especially for CEA 3820R). It is clear that volumetric strains are underestimated. The reason could be that M e was estimated using the friction angle in compression in Equation 4.25.

ε d [%] ∆ q [kPa] CEA3820R 2 3 4 5 ln p [kPa] ε v [%] ¨5 0 5 ©60 30 0 30 ε d [%] ∆ q [kPa] CEA3516R ln p [kPa] ε v [%] ln p [kPa] ε v [%] ε d [%] ∆ q [kPa]

Bothkennar clay

An extensive experimental investigation have been performed by [START_REF] Smith | The yielding of bothkennar clay[END_REF] and [START_REF] Allman | Mechanical properties of reconstituted bothkennar soil[END_REF] in order to study the behaviour of the sensitive Bothkennar clay. The resulting data have been used for calibration and validation of numerous constitutive laws for natural and intact clays [START_REF] Rouainia | A kinematic hardening constitutive model for natural clays with loss of structure[END_REF], Taiebat et al. (2010b), Wheeler et al. (2003a), Karstunen and Koskinen (2004b) among others). In the following, we simulate some laboratory tests on intact and reconstituted samples in order to test the performance of the SCLAY-1SB model regarding to the sensitivity of intact clays. The tests include oedometer experiments as well as drained and undrained triaxial tests on reconstituted and intact samples. The calibration method is based on the procedure described in section 3 and the same set of parameters have been considered for all tests (see Table 4.4).

One-dimensional consolidation tests

The one-dimensional consolidation behaviour of Bothkennar clay has been studied by testing high quality intact and reconstituted specimens with two types of samples from Sheerbrooke and Laval Universities. According to [START_REF] Smith | The yielding of bothkennar clay[END_REF], Sherbrooke and Laval samples are two types of intact specimens with the former appeared to suffer marginally less disturbance (de-bonding).

The initial values of the internal variables used for the simulation with the SCLAY-1SB model are given in Table 4.6. The comparison of the model simulations with the experimental results are presented in Figure 4.8. The model response is fairly consistent with observations and a complete destructuration is correctly obtained at very large stresses. Sherbrooke and Laval tests have been simulated. It can be observed that the results do not show the accuracy that was obtained in one-dimensional consolidation, but they are still acceptable. The difference is more pronounced in extension than in compression tests, which may be explained by the complexity of the soil behaviour, that is at the same time anisotropic and sensitive. Plastic anisotropy is also reproduced relatively well, since the stress state at the critical state is close to the experimental data.

Conclusion

The bounding surface formulation is incorporated into the SCLAY-1S model to simulate plastic strains within the conventional yield surface during monotonic loading. The proposed model, named SCLAY-1SB, captures several important features of clays under monotonic loading as for example the response of overconsolidated clays, the initial and induced anisotropies and the sensitivity of intact and naturally deposited soils. Formulation of the SCLAY-1SB model is presented both for triaxial and multi-axial spaces. A method to calibrate the model parameters is also provided.

Validation of the model is given though comparison with experimental data. Reconstituted and intact as well as normally and overconsolidated samples have been simulated. Generally, reasonable results are obtained for the various tests. Simulations of the structured Bothkennar clay as well as the drained triaxial tests of London clay and Otaniemi clays have show good agreement with experimental curves.

The next chapter deals with the application of the SCLAY-1SB model to reproduce the non linear behaviour of offshore monopiles submitted to monotonic lateral loading.

Introduction

This chapter is about the hydro-mechanical modelling of offshore monopiles embedded in clayey soil subjected to a monotonic lateral loading. The SCLAY-1SB constitutive law is adopted for the simulations.

Monopiles are the most common foundations used to support offshore structures, due to their lower cost and simplicity. Their installation procedure is considered as low-risk and robust for most soil conditions. Comparing to relatively new foundation concepts, design methods for monopiles are simple and well detailed in design codes. Nevertheless, non-linear soil-pile interaction is still an object of study.

Most of offshore piles are designed to resit to lateral loading, whether monotonic or cyclic. When subjected to horizontal forces, pile foundations resist by bending, while still supporting the vertical load transmitted by the superstructure. Their mechanical behaviour is such that they resit a lateral load by mobilizing passive pressure in the surrounding soil. The distribution of soil's reaction depends on the stiffness of the pile and the soil, as well as on the fixity (boundary conditions) at the pile ends.

The most common design method for laterally loaded piles is the p-y spring technique, based on the Winkler approach, which models the load transfer between the pile and soil by a series of non-linear springs along the pile following uniaxial laws (p-y curves).

This chapter opens with a brief state of the art, presenting the main features of the offshore monopiles, such as the ultimate bearing capacity, the problem of soil-pile interaction, and some of the most popular p-y spring models. Afterwards, finite element analyses are presented. First, a two-dimensional finite element model is introduced, aiming to re-visit the p-y curves for natural clays and to cover several points such as the loading rate, soil and stress anisotropies and cohesion degradation of stiff clays. Second, a three-dimensional finite element analysis is carried out, where a centrifuge test of a laterally loaded pile is modelled. Lastly, the chapter ends with a numerical comparison of slender and short monopiles, by the mean of 3D Finite Element Modelling.

2 Monopiles in offshore wind turbine parks Nowadays, offshore wind industry is exponentially developing due to the significant contribution of the North Sea Wind Turbines sector to the energy production in Europe. Given that the EU's energy target is to increase the share of renewable energy by 2020, there is a great potential of offshore wind energy applications. In order to achieve this goal, the cost of the wind turbine parks should be significantly reduced. This objective can be pertly reached through optimisation of existing design methods and also by introducing new technologies. One of the areas where costs can be reduced is the support structure, knowing that their cost can attain 30% of the total cost of a wind turbine park project (EWEA ( 2017)).

The type and cost of offshore foundations are strongly related to the distance from the shore, which influences the width of the support structure and the geotechnical conditions of the seabed and the sea depth. The latter is one of the most important factors influencing the viability of an offshore wind farm, since the foundation cost significantly increases with the increase of the sea depth. Therefore, the prevailing support structure design for OWTs is by so far the monopile foundation. According to EWEA (2017), 90% of offshore wind farms adopted this type of foundation due to its multiple advantages; simple design and massive fabrication, an installation approach based on conventional impact driving which is relatively low-risk and robust for most soil conditions. 108CHAPTER 5. MONOPILES UNDER MONOTONIC LATERAL LOADING

Loads acting on offshore structures

Offshore wind turbines are often placed in adverse environmental conditions, with strong wind and wave loading. Figure 5.2 (a) shows the typical distribution of a wind turbine and wave pressures along the length of the tower. According to [START_REF] Bhattacharya | Similitude relationships for physical modelling of monopile-supported offshore wind turbines[END_REF], two types of lateral loads are distinguished as follows:

ii. cyclic loading is applied the top of the tower, as a result of wind shielding and the blade passing effect, due to the periodical obstruction caused by the spinning of the blade.

i. the bottom part of the tower, which is unobstructed by the spinning turbine blade, is subjected to monotonic force, coming from wind loads. The total environmental lateral load acting on an offshore wind turbine can be approximated as a horizontal load H of monotonic and cyclic components, acting at a distance d above the foundation level (see Figure 5.2 (b). H represents the resultant lateral load applied on the tower that must be resisted by the monopile foundation. In addition to the horizontal loads, the tower is subjected to a vertical force coming from the usual gravity load. However, the vertical load is usually negligible with respect to the lateral load since the OWTs have often a lightweight structure [START_REF] Bhattacharya | Similitude relationships for physical modelling of monopile-supported offshore wind turbines[END_REF].

Therefore, the design of OWT foundations is mainly controlled by the lateral load components. Efficient foundations can be designed by addressing uncertainties in the lateral loading design process (for example, the different types, magnitudes and localisations of loads).

Soil-pile interaction

When a monopile is subjected to horizontal forces, the equilibrium of the system (superstructure-pile-soil) is reached by the mobilization of the passive and eventually active pressures in the surrounding soil. The distribution of soil's reactions depends on the fixity of the pile ends and the relative stiffness of the soil-pile system, denoted K.

Theoretically, the notion of the relative stiffness allows distinction of slender and rigid piles. Under lateral loading, the rigid piles rotate without bending, whereas slender piles resist by bending.

The relative stiffness K is generally expressed in function of the transfer distance, l 0 [m] as follows:

l 0 = 4E p I p E s 0.25 (5.1)
where E p and E s are the elastic moduli of the pile and the soil, respectively, and I p is the pile inertia moment. For a circular pile K is hence expressed as:

K = 4 B l 0 4 (5.2)
where B is the pile diameter. Numerous recommendations of limits between rigid and slender piles are proposed in the past. For the theoretical case of an elastic and homogeneous soil, some of these recommendations are given in Table 5.1, function of the relative stiffness K. It is noteworthy however to mention that the choice is not trivial, because of the dispersion of the limit values and the unrealistic and simplified assumptions of the soil behaviour.

Ultimate lateral bearing capacity of piles

Two mechanisms contribute to the soil resistance around a laterally loaded pile, as illustrated in Figure 5.3. At shallow depth, soil fails in a conical wedge that extends to the surface. At a certain depth where soil failure is no longer influenced by the ground surface, soil failure corresponds to a localized flowaround mechanism, i.e., the soil flows around the pile in a non-gapping plane strain mechanism. At shallow depth where the conical wedge prevails, a gap may appear at the soil-pile interface in the active side (back side, see Figure 5.3) of the foundation. In this case, the resistance due to the shear strength and weight of the passive wedge in front of the pile is the only one mobilised. Otherwise, if a gap does not form, both active and passive wedges are mobilised. According to [START_REF] Jeanjean | A framework for monotonic py curves in clays[END_REF], these conditions can be referred also as no suction (when gapping) and suction (case of no gapping) conditions.

The lateral bearing capacity, p u , is usually expressed as function of the pile diameter B and the undrained shear strength of the soil, s u :

p u = N p .s u .B (5.3)
where N p is the factor of lateral bearing capacity. The undrained soil shear strength is related to the soil properties and the mode of shearing, whereas the value of N p results from the interaction of the soil-pile system and hence depends on the type of the failure mechanism. In the literature, the limit equilibrium and then plastic limit analyses were performed to provide values of N p . Experimental studies on real size foundations and laboratory tests were furthermore employed as well as finite element and finite difference methods. Several pioneer studies are presented hereinafter both for the flow-around and conical wedge mechanisms.

Flow-around mechanism

• Limit equilibrium analysis Reese (1958) performed a plane strain limit equilibrium analysis, where a square pile was assumed surrounded by eight soil elements with prescribed deformation patterns, as illustrated in Figure 5.5. As a result, the bearing capacity factor N p was derived, Reese found N p = 12 for fully rough conditions and N p = 10 for smooth contact surfaces.

Pile Load Figure 5.4: Limit equilibrium analysis of a square pile by [START_REF] Reese | Discussion on soil modulus for laterally loaded piles[END_REF], after [START_REF] Jeanjean | A framework for monotonic py curves in clays[END_REF] • Plastic limit analyses A plastic limit analysis of the laterally loaded pile problem was firstly proposed by [START_REF] Randolph | The limiting pressure on a circular pile loaded laterally in cohesive soil[END_REF]. They considered an isotropic, homogeneous and perfectly plastic soil, and developed a plane strain upper-bound and lower-bound solutions for the ultimate lateral bearing capacity. Under these assumptions, values of N p were derived for different interface roughness. The derived values are summarized in Table 5.2. Strictly speaking, the solution is exact for fully rough interfaces, i.e. r = 1.0, since it was demonstrated that the upper and lower bounds are in this case identical. The proposed solutions for r < 1 correspond to the lower-bound solution.

Refinement of the N p values, especially for the upper-bound solution have been discussed by [START_REF] Murff | P-ultimate for undrained analysis of laterally loaded piles[END_REF], [START_REF] Martin | Upper-bound analysis of lateral pile capacity in cohesive soil[END_REF] and [START_REF] Klar | Upper bound for cylinder movement using "elastic" fields and its possible application to pile deformation analysis[END_REF] among others. The upper-bound solution of [START_REF] Martin | Upper-bound analysis of lateral pile capacity in cohesive soil[END_REF] was found very close to the [START_REF] Randolph | The limiting pressure on a circular pile loaded laterally in cohesive soil[END_REF] lower-bound solutions. This finding led them to conclude that the solution of Randolph and Houlsby, i.e., values in Table 5.2 could turn out to be exact for every value interface roughness.

• Guidelines of API and standards ISO codes

The majority of the offshore geotechnical structures design codes follow the recommendations of [START_REF] Matlock | Correlations for design of laterally loaded piles in soft clay[END_REF], [START_REF] Broms | Lateral resistance of piles in cohesive soils[END_REF] and [START_REF] Reese | Discussion on soil modulus for laterally loaded piles[END_REF], for determining the value of N p for the flow around mechanism.

According to the last edition of API in 2003 API ( 2003), the ultimate bearing capacity factor of clays at depth, varies between 8 and 12. In particular, API proposed the following recommendations:

-Reese: N p = 10 for r = 0; N p = 12 for r = 1;

-Meyherhof: N p = 8.28 for r = 0; N p = 8.85 for r = 1;

-Tschebotarioff: N p = 11.4;

-Mackenzie: N p = 8.5

The value of N p = 9 was recommended by the API in the absence of specific criteria, according to the recommendations of [START_REF] Matlock | Correlations for design of laterally loaded piles in soft clay[END_REF].

The same tentative values were later adopted by [START_REF] Broms | Lateral resistance of piles in cohesive soils[END_REF].

Shallow wedge failure mechanism

Similarly, the soil resistance due to the wedge mechanism has been investigated using limit equilibrium and plastic limit analyses and empirical correlations. Several studies are summarised hereinafter.

• Limit equilibrium analysis

Based on a limit equilibrium analysis of the passive wedge, [START_REF] Reese | Discussion on soil modulus for laterally loaded piles[END_REF] derived the lateral bearing capacity factor N p for the theoretical case of Reese calculated all the forces acting along the inclined bottom of the wedge, for a specific angle of β = 45 • , and a pile square geometry when calculating the axial friction. In this specific case, the following equation of N p have been found in function of the depth z, density γ, pile diameter B, contact roughness r, and undrained shear strength s u :

N p = N p0 + N pw (5.4) where:

N p0 = 2 + 2.83 z B
for r = 0 (5.5)

N p0 = 3 + 2.83 z B for r = 1 (5.6) N pw = γz s u (5.7)
The first term in Equations 5.5 and 5.6 represents the resistance on the bottom of the wedge and the side of the pile, whereas the second term represents the soil resistance along the side of the wedge.

• Empirical correlations

Bearing capacity factors have also been derived using curve fitting of 1-g load tests results or laboratory tests combined with finite elements analyses (see [START_REF] Jeanjean | A framework for monotonic py curves in clays[END_REF] or the recommendations of the SOLCYP project based on centrifuge tests [START_REF] Puech | Design of Piles Under Cyclic Loading: SOLCYP Recommendations[END_REF]).

The recommendations of [START_REF] Matlock | Correlations for design of laterally loaded piles in soft clay[END_REF] are the basis of the API and ISO standards for soft clays. He used experimental results of a series of pile tests, conducted at Lake of Austin and Sabine River and he modified the Reese's framework for a fully rough pile. The N p factor is expressed as:

N p = 3 + J. z B + γz s u (5.8)
where J is model constant. for data of Lake Austin and Sabine River tests, the value of J = 0.22, and J = 0.5, for rough and smooth contact are recommended, respectively.

p-y curves for natural clays

In practice, analysis of laterally loaded pile is generally performed using the non-linear Winkler spring model, often known as p-y model, due to its simplicity, low computational cost and the ability to model layered soils. Following this approach, the pile is modelled as a beam and the soil-pile interaction is characterized by a set of independent, unidimensionnal and non-linear springs.

Each spring obeys to a p-y curve, where p is the soil reaction per unit of length that causes a deflection y of the corresponding pile section (see Figure 5.6). The p-y curves have been first applied in the gas and oil petroleum industry, for thin and slender piles (API-2003). The model calculates the foundation stiffness, the displacement at the pile head and the bending moments along the pile. More recently, new p-y curves for monopiles of large diameter (XL) have been proposed following the international PISA project Byrne et al. (2015b). They take into account the lower L /B ratio (where L being the pile embedded length), the larger diameter and the rigid body rotational movement of the XL monopiles.

The p-y curves depend on a number of factors, such as soil properties, pile stiffness, relative soil-pile stiffness and pile ends conditions. A brief review of some common p-y curves is presented hereafter.

• Matlock p-y curve for soft clay

The monotonic p-y curve for soft clay proposed by [START_REF] Matlock | Correlations for design of laterally loaded piles in soft clay[END_REF] originates from field testing, performed on lightly over-consolidated clays at Sabine River. The monotonic p-y curve relates the mobilised lateral resistance to the normalised lateral displacement with a power law relation:

p = 0.5p u y y c 0.33
; with y c = 2.5ε 50 .B (5.9)

where p u represents the ultimate lateral resistance per unit pile length, y the pile deflection, and B is the pile diameter. ε 50 is the axial strain at which 50% of maximum shear stress is mobilised in an undrained compression test. Although the API code does not specify the type of test to perform, unconsolidated undrained compression tests are often used in practice.

• Reese p-y curve for stiff clay

The p-y curve for stiff clays presents a peak shape due to the degradation of their mechanical properties. [START_REF] Reese | Field testing and analysis of laterally loaded piles om stiff clay[END_REF] explained the postpeak softening behaviour by the cohesion degradation and proposed a specific procedure to elaborate the p-y curve. As shown in Figure 5.7, Reese defined the p-y curve by dividing the curve into four parts. Initially, the curve is linear until a displacement of y = y c , followed by a non-linear increase of the soil reaction until the peak corresponding to a displacement of y = 6A s .y c . Then, the soil

The ultimate lateral bearing capacity per unit pile length p u , calculated as N p .s u .D is also related to the interface roughness factor by : N p = 9 + 3r (5.12)

Undrained shear strength from the SCLAY-1SB

As demonstrated herein-before, the computation of p-y curves requires a close estimation of the ultimate pile lateral bearing capacity, defined as p u = su B , where B is the pile diameter and s u the undrained shear strength of the clay. Therefore, the notion of undrained shear strength plays an important role for determining the undrained p-y curves for piles embedded in clayey soils.

Two numerical approaches are used to estimate the undrained shear strength s u : total stress-based and effective stress-based methods.

• The former group addresses the non-linear stress-strain relationship by empirically fitting experimental data. The total stress-based models employ generally a Von-Mises criterion, where the soil is defined with its Young modulus, Poisson's ratio (generally equal to ν = 0.49 to ensure constant volume conditions) and undrained shear strength. The basic disadvantage of these models is that the stress-strain response is decoupled from the pore-pressure generation, any hydro-mechanical coupling is therefore neglected.

• The effective stress-based models used for clay are mainly based on elastoplasticity and necessitate non-linear incremental solution schemes. Furthermore, the complexity of the constitutive laws and of the solidpore water interaction mechanism impose significant computational requirements. Finally, the undrained shear strength is not a direct input of these models; s u is indirectly determined and may depend on several parameters such as the coefficient of earth at rest K 0 , the overconsolidation ratio, the friction angle and dilatancy. For all these reasons, the effective stress-based models are not often used in engineering everyday practice, even though they allow a deep understanding of the complex soil-structure interaction problems, as well as the effect of some fundamental clay features such as anisotropy, dilatancy and consolidation [START_REF] Ahayan | Behaviour of laterally loaded pile[END_REF]).

The present research work deals with effective stress-based models. In particular, the SCLAY-1SB constitutive law is employed hereafter to reproduce the behaviour of laterally loaded piles. Therefore, the equation of the undrained shear strength is firstly derived before going further in the numerical analyses. .9 presents the idealized results of an undrained triaxial compression test reproduced with the SCLAY-1SB model for the case of an overconsolidated, anisotropic and intact clay sample. The specimen's state is defined with: the stress invariants (p ,q,θ), the state parameters (e,p mi ,χ, p t ) and the consolidation parameters κ and λ i . The undrained compression path on the specimen brings the specimen from Point J to failure at point L on the critical state line, at the same void ratio, i.e: e L = e J (5.13)

Points K and L both lie on the critical state line. It follows that:

e L -e K = λ ln p L p K (5.14)
On the other side, K and J lie on the swelling line, hence:

e J -e K = κ ln p J p K (5.15)
Combination of the previous three equations leads to the following relations:

-κ ln 5.16) which leads to:

p K p J = λ ln p L p K = λ ln p L p J + λ ln p J p K (
p L p J = p K p J (1-κ λ )
(5.17)

Knowing that K lies on the yield surface, the following equation is verified:

(M -α) 2 + M 2 -α 2 1 - p m p K = 0 (5.18) p m p K = 1 + (M -α) 2 M 2 -α 2 (5.19)
Therefore:

p m p K = 2M M + α (5.20)
Given that Point I refers to the normally consolidated state and lies on the virgin consolidated line, i.e η I = η K0 and:

(η K0 -α 0 ) 2 + M 2 -α 2 0 1 - p m p I = 0 (5.21) p m p I = 1 + (η K0 -α 0 ) 2 M 2 -α 2 0 (5.22)
where α 0 is the initial value of the state parameter α. This equation leads to:

p m = 1 + (η K0 -α 0 ) 2 M 2 -α 2 0 p I (5.23)
On the other hand:

p I p J = 1 /3 (1 + 2K 0 ) σ c 1 /3 (1 + 2K 0 ) σ v0 = σ c σ v0 = OCR (5.24)
where σ c is the preconsolisation pressure.

Finally, the combination of the two previous equations leads to the following expression:

p L = M + α 2M 1 + (η K0 -α 0 ) 2 M 2 -α 2 0 1-κ λ 1 + 2K 0 3 σ v0 (5.25)
The undrained shear strength under triaxial and plane strain conditions can be calculated as:

2D-FEA for a laterally loaded pile

The motivation of this section is to present a new framework for calculation of best-estimate p-y curves, based on finite element method. In particular, the influence of the hydro-mechanical behaviour of soil on the shape of p-y curves and on the lateral bearing capacity factor is investigated. The study is performed using a fully hydro-mechanical coupled approach within a 2D plane strain model. The plane strain simplification is justified, at least near the pile for depths where the flow-around failure mechanism prevails.

The effect of the loading rate is first discussed. Then, the influences of different clay properties (i.e, soil and stress anisotropies, clay structure and cohesion) on the N p factor and the shape of the soil reaction curves are extensively studied.

Finite element model

Geometry and boundary conditions

The 2D finite element model is shown in Figure 5.10. The model consists of a square section of soil with a rigid disc having a diameter B, representing the pile cross-section to its center. Only a half-pile segment is modelled due to the symmetry of the problem. The outside boundary of the model is fixed, while a horizontal load H is applied to the pile section, causing a rigid displacement of the cross-section along the axis of loading. The pile is loaded until reaching the ultimate state, when the load-displacement curve becomes flat. A coupled hydro-mechanical analysis is considered and a constant water pressure is imposed in the external boundary and inside the pile section. During loading, pore water pressures are thus able to dissipate (and/or generate) within the soil and soil-pile interface.

The outside boundary is extended to a distance of R = 7.5B from the pile section. Note that the location of the external boundaries is not trivial for 2D-laterally loaded pile problems, as it can affect the initial modulus of the load-displacement curve. This problem has been covered by many authors in the past [START_REF] Baguelin | Theoretical study of lateral reaction mechanism of piles[END_REF] and [START_REF] Klar | Upper-bound and load-displacement solutions for laterally loaded piles in clays based on energy minimisation[END_REF], among others), using either plastic limit methods, or finite element and difference finite methods. For instance, [START_REF] Baguelin | Theoretical study of lateral reaction mechanism of piles[END_REF] studied the case of an elastic soil and found that the initial soil reaction modulus decreases with the distance R from the outside boundary. To avoid this, the auhors proposed values of R depending on the pile diameter, the transfer distance of the soilpile system and the pile fixity conditions. The selected value of R = 7.5B corresponds to a flexible pile of diameter B = 0.9m with a free head, according to [START_REF] Baguelin | Theoretical study of lateral reaction mechanism of piles[END_REF].

Soil elements

For all the calculations presented hereafter, the soil is modelled with the SCLAY-1SB law, unless otherwise mentioned. Depending on the aim of the analysis different variants of the constitutive law are employed. Nevertheless, the same soil parameters referring to the Otaniemi clay and summarized in Table 5.3, are used.

Interface elements

Interface finite elements are introduced between the soil matrix and the pile section. The interface elements belong to the zero-thickness family (see Chapter 2, Section 2). The probable zone of contact respects the ideal contact constrain states:

g N ≥ 0, p N ≥ 0 and p N .g N = 0 (5.28)
where p N is the normal effective pressure and g N the gap between the two sides of the interface, as shown in Figure 5.11. This ideal constrain establishes that contact holds if the gap is equal to zero (closed), giving rise to a normal contact pressure. Otherwise, if the gap is positive (open), there is no contact pressure. The shear behaviour of the interface is described in a similar way. The maximum shear resistance along the interface is ruled by the Mohr-Coulomb criterion, according to:

τ max = σ n tan ϕ int + c int (5.29)
where φ int and c int are the steel-soil friction angle and cohesion, respectively. These parameters are selected depending on the nature of the roughness of the soil-pile surface. According to [START_REF] Kostkanova | Measurement of wall friction in direct shear tests on soft soil[END_REF], a fully rough soil-steel interface corresponds to ϕ int = 2 /3ϕ, where ϕ is the friction angle of the soil surrounding the pile. If the shear stress is lower than τ max , soil and pile are considered "stuck", and the relative tangential displacement g T is null. If the maximal shear stress is reached, both sides of interface encounter a relative displacement.

Fluid flow through and along the interface is captured due to the threenode discretisation of pressures (see Figure 5.11). By doing so, three different flows are distinguished:

• Transversal flow through each side of the interface. Between each side and the inside of the interface, f T w is related to the variation of pore water pressures as follows:

f t w = ρ w .T w .∆p w (5.30) where ρ w is the specific mass of water and T w is a transversal conductivity.

• The longitudinal flow along the interface, respecting the Darcy's law:

f l w = -ρ w k l µ w (∆p w + ρ w .g.∆z) (5.31)
where k l is the longitudinal permeability and µ w is the dynamic viscosity of water. k l = 10 -9 m 2 is selected for the following analyses.

• The storage flux, coming from the conservation mass equation, provided from the progressive opening (or closure) of the gap.

The first coupling between mechanical and flow quantities stems from the Terzaghi principal, by defining a total pressure p N acting on each side of the interface as p N = p N +u w , where u w is the pore water pressure in the interface. The second coupling term comes from the storage flux, that relay the variation of pore water pressure to the opening (or closure ) of the gap.

Validation

In order to validate the finite element mesh, the numerical value of N p is compared to the analytical solutions of [START_REF] Randolph | The limiting pressure on a circular pile loaded laterally in cohesive soil[END_REF]. A total stress approach is adopted and the theoretical case of an elastic, perfectly plastic soil is treated. Conditions of constant volume are imposed (ν = 0.499) and the criterion of Von-Mises, with constant undrained shear stress, is used both in soil and contact elements.

Table 5.4 compares the computed lateral bearing capacity factor N p with the theoretical solutions of [START_REF] Randolph | The limiting pressure on a circular pile loaded laterally in cohesive soil[END_REF] for different interface roughnesses. Under the adopted conditions, the contact roughness is equal to the ratio of the soil and soil-pile interface strengths. The maximum numerical error is of 8.07% for the smooth interface (r = 0) and the minimum 0.64% when the interface is fully rough (r = 1). The numerical error can be attributed to the mesh refinement but it is not corrected in the rest of the numerical analyses since the computed error is largely acceptable.

Effect of the loading rate

The loading rate may have significant effects on the behaviour of laterally loaded piles, especially for offshore foundations where the soil is always saturated. Full-scale tests on a laterally loaded pile in fine soil, presented by [START_REF] Tassios | Efforts repetitifs horizontaux sur pieux verticaux[END_REF], demonstrated that the loading rate affects the shape of the global response of the pile; the head force-displacement curve is stiffer when the load is applied faster.

In order to avoid introducing this additional parameter to the present study, undrained conditions are hereafter considered. This assumption is commonly considered in clayey soils, since clay permeability is generally very low compared to the conventional rate of loading. However, a clear limit between drained, partially drained and fully undrained conditions is not clear. Furthermore, the drainage conditions depend on the consolidation properties, such as the intrinsic permeability and the coefficient of vertical consolidation, which in its turn depends on the mean effective stress. [START_REF] House | Evaluating the coefficient of consolidation using penetration tests[END_REF] treated this problem based on t-bar penetration tests. In particular, the authors quantified the evolution of the normalized soil strength as function of the normalised rate of loading as follows:

V = v.B c v (5.32)
where B is the diameter of the t-bar and c v the vertical consolidation coefficient. They concluded that when the loading rate decreases, the strength increases since the pore water pressures evolve. Beyond a certain value of V , the strength remains unchanged. At this stage, the undrained conditions are reached.

Based on these results, they provided limits for drainage conditions. For kaolin clay the values are: V < 0.2 for drained conditions, V > 20 for undrained conditions and 0.2 < V < 20 for partially drained conditions. Nevetheless, these values cannot be applied directly on the current study, since the soil parameters correspond to the Otaniemi clay rather than Kaolin clay.

The rate of loading can be taken into consideration in two ways in a finite element analysis: (i) with rate-dependent constitutive laws, which is out of scope of the present research, and/ or (ii) with a coupled solid-pore water pressure approach in saturated or partially saturated medium. The present study stands with the latter case, knowing that different variants of the SCLAY-1SB model with the Terzaghi effective stress are used throughout this research. The undrained conditions must therefore be verified before going further.

To achieve this aim, simulations have been carried out using different values of the normalized loading rate V (V = 0.1;V = 1.0 and V = 30), in addition to drained conditions, by fixing the Pore Water Pressure (PWP) in all the mesh nodes. Results are presented in Figure 5.12 in terms of normalized p-y curves. The p-y curves are normalised with the same value of the undrained shear strength, corresponding to a normally consolidated, isotropic and reconstituted soil,i.e., α 0 = 0, K 0 = 1, A = 0, OCR = 1, χ 0 = and c = 0. The other soil parameters are summarized in Table 5.3. As expected, Figure 5.12 clearly shows the influence of the loading rate .12: p-y curves computed wth different loading rates on the p-y curves. The soil reaction modulus strongly depends on the loading rate; the response is stiffer when the loading rate is higher and beyond a certain limit (for V > 20) the p-y curve is independent of the loading rate. On the other hand, the lateral bearing capacity seems invariable, since the different curves tend toward the same ultimate value (the same plateau). Figure 5.13 (a) presents the variation of the effective normal pressures at the soil-pile interface, normalised with its initial value. The figure reveals the influence of the loading rate on the interface behaviour. In particular, the gap at the soil-pile interface, corresponding to a null value of the effective normal pressure, is strongly affected by the loading rate. For low values of the normalized velocity, i.e. V < 1.0, including drained conditions, zero-valued effective normal pressures are observed in the active side of the pile section. p n remains positive for V ≥ 1.0, which means that the pile displacement redistributes the contact pressure without any gap, leading to no volume change conditions. Figure 5.13 (b) depicts the distribution of PWPs at the soil-pile interface, normalized with their initial values, for the different velocities. The evolution of PWP around the pile varies with the loading rate. Generally, excess PWP is observed in the passive side of the pile, while a decrease of PWP is revealed in the active side. For V > 30 the PWP is the most mobilized, compared to PWP computed with smaller values of V .

It is also interesting to observe the variation of PWP in the soil near the pile section. Figure 5.14 depicts the PWP increments for undrained conditions (with V = 30), at different loading stages, i.e. at the beginning of loading y = 0.01B, y = 0.08B and at the ultimate state for y = 0.15B. The pile displacement redistributes the PWP in the soil such that a decrease of PWP is developed at the active side and an excess of PWP appears in the passive side of the pile. In addition, by being closer to the ultimate state the negative increments of PWP decrease until reaching a minimum value when the loaddisplacement curve reaches the plateau (ultimate state). This behaviour can be explained by the failure mechanism of the soil at the ultimate state. When in failure, the pile displacement induces solely the redistribution of pore water pressure without any increase of the pile capacity (when the flat plateau is reached).

Effect of clay features on the p-y curves

Undrained conditions are considered all along the following analysis. A lateral displacement is applied to the pile section with a constant velocity of v = 3cm/s, which corresponds to V = 30, knowing that the soil permeability is k = 10 -9 m/s. The soil-pile interface is considered smooth and the soil is modelled with different variants of the SCLAY-1SB law. The obtained p-y curve is normalised with the undrained shear stress, calculated from Equation 5.27 under plane strain conditions.

Effect of anisotropy

The results of the previous section were about an isotropic material. However, it is known that clay exhibits an anisotropic behaviour, coming from the anisotropy of stress due to the virgin consolidation and to the anisotropy induced by the pile loading. The question faced by design engineers is how to account for the initial and induced anisotropies, or in other words, are the values of N p and the shape of the p-y curves affected by soil anisotropy. Another question is what is the value of the shear strength to use to compute the normalized p-y curve.

The effects of initial and induced anisotropies are studied hereinafter. First, the effect of K 0 -consolidation only is analysed. Then, the effect of induced anisotropy is studied and finally a p-y curve is computed considering both stress and soil anisotropies. For each case, the p-y curve is normalized using the corresponding undrained shear strength.

Figure 5.15 (a) presents the normalized p-y curves, computed for different values of K 0 . The lateral bearing capacity factor seems to be invariant with the value of K 0 , since all the p-y curves tend toward the same plateau. However, the shape of p-y curves is strongly dependent on K 0 since the soil reaction modulus increases with increasing K 0 . 5.15 (b) shows the effect of soil anisotropy on the soil reaction curves. In this case, the initial stress conditions are considered isotropic while the rotational hardening rule (RHR) is activated to capture the anisotropy induced by the pile loading. The RHR is controlled by the parameter A that predicts the evolution of the internal variable α. The resulting normalized p-y curves have a closed shape and almost the same value of the bearing capacity factor. Furthermore, a smooth softening is observed, especially for the higher value of A. The softening behaviour may be attributed to the degradation of anisotropy when the soil reaches the failure conditions at the critical state.

Finally a more realistic case, where both soil and stress anisotropies are activated, is presented hereafter. The computed soil reaction curve is compared to the cases where initial and induced anisotropies are considered separately, see Figure 5.16. Again, the soil reaction curves tend toward the same plateau. It is concluded therefore that the bearing capacity factor is not affected by the stress or the soil anisotropy. However, the lateral ultimate capacity p u depends on the anisotropy, since the undrained shear strength is given in function of the anisotropy parameters.

Effect of cohesion and structure

In the following analysis the anisotropy is deactivated and the soil is considered normally consolidated. Only features related to the clay structure are activated, i.e, the state parameters p t , χ and their respective softening rules. Figure 5.17 (a) shows soil reaction curves computed with different values of initial cohesion and initial state of the structure. The results reveal the following points:

• When the soil softening is activated (for non-zero values of cohesion or χ), the p-y curves do not present a plateau, soil reaction decreases after reaching a peak value.

• Both the peak and the post-peak curve depend on the amount of cohesion and the state of structure. Furthermore, the peak decreases with increasing c and χ.

• The lateral bearing capacity factor, which is defined in this case with the initial values of c , is smaller than the N p value, computed when softening is deactivated, i.e. for c = 0 and χ = 0.

These results are convenient with the observations made by [START_REF] Reese | Field testing and analysis of laterally loaded piles om stiff clay[END_REF] for the p-y curve of stiff clays that relate its shape to the degradation of the cohesion. The 2D model provided interesting results regarding the response of laterally loaded piles. However, it should be kept in mind that the plane strain assumption is realistic for high depths, where the flow-around mechanism is activated. The response of a laterally loaded pile is actually a 3D problem, where several phenomena are interacting, such as the stress distribution surrounding the pile, the gap opening between the soil and the foundation near the surface, the soil deformation near the ground surface and the water flow in the soil and the soil-pile interface.

Accordingly, the following section presents a 3D modelling of a single pile embedded in soft clay subjected to monotonic lateral load. In particular, a centrifuge test is modelled using a fully hydro-mechanical approach and the SCLAY-1SB constitutive law. Details of the centrifuge test and the finite element model are presented hereinafter, followed by the calibration procedure and the results of the analysis.

Centrifuge model description

Within the framework of the PN -ANR SOLCYP project, [START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF] reported detailed centrifuge tests of laterally loaded piles embedded in Kaolin clay. The centrifuge tests were carried out on reduced models to 1 /50 th , in macro gravity of 50.g.

Geometry and instrumentation

A sketch of the centrifuge test is provided in Figure 5.18. A free-head pile with a hollow section is selected for the finite element analyses. Table 5.5 summarizes the geometries of the pile in the model and prototype scales (reduced and real scale, respectively) . The pile is instrumented with 21 levels of strain gauges to measure the bending moment along the pile. According to [START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF], the diameter of the pile is enlarged by 6% because of the coating on the strain gauges. Nevertheless, the pile stiffness remains unchanged, i.e., the bending stiffness E p .I p verifies the similarity rule between the model and the prototype ( 1 /50 4 ) .

Pile installation and loading

The pile was installed at 1.g by drilling. A slightly oversize hole was first created (2mm oversized), the pile was then placed into the pre-drilled hole.

Afterwards, the small soil-pile gap was closed by hand at the surface and supposed to be closed by the self-weight of the soil along the pile surface when ramped up on flight due to the soil deformation. The lateral load was applied at the pile head, at 2m of the mud-line. The monotonic tests were run in displacement-controlled mode at a constant velocity of 0.4mm/s (model scale), corresponding to 20mm/s at prototype scale.

Soil conditions

Normally consolidated Kaolin clay was reconstituted for the tests. Two types of soil conditions were prepared: saturated, slightly over-consolidated clay, and unsaturated and heavily overconsolidated clay. The numerical analysis of this manuscript deals with the former case.

The soil sample was prepared from Kaolin slurry mixed under vacuum to 90% of water content. To achieve the desired sample height, the samples were prepared in three layers, and then consolidated at 1.g and on flight at 50.g, before the pile installation.

Oedometer tests were performed in order to measure the soil consolidation parameter, and to determine the initial depth of the soil sample. Table 5.6 summarizes the consolidation parameters, [START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF]. Table 5.6: Consolidation parameters of Kaolin soil, from [START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF] ρ s (kg/m 3 ) C v (m 2 /s) n (-) w (%) κ (-) λ (-) 2650 3.10 -7 90 54 0.027 0.143

The consolidation pressure at 1.g was chosen according to the value of OCR targeted in the middle of each layer at flight. The OCR is defined as the ratio of the consolidation pressure applied at 1.g and the effective stress at 50.g:

OCR = σ c σ v 50g (5.33)
As shown in Figure 5.19, the soil sample is lightly overconsolidated, the OCR is close to 2 except in the shallow layer, where the OCR is around 7. The undrained shear strength was established from an in-flight cone penetration test using the following empirical correlation:

s u = q c 18.5
(5.34) where q c is the limit capacity of the CPT tests carried out under 1.g, and in flight, before the pile installation.

It is unclear which shear mode the correlated shear strength corresponds to. Hence, the calibration of the soil parameters is carried out directly on

Finite element model description

Geometry and boundary conditions

Figure 5.20 shows details of the geometry, meshing, and boundary conditions of the finite element model used to simulate the centrifuge test. Taking advantage of the symmetry, only half of the pile and the soil domain is discretised. Both horizontal and vertical displacements are blocked all along the lateral boundary of the finite element mesh, while the vertical displacements are fixed at the bottom. The dimensions of soil domain and the pile are identical to the prototype.

A fully hydro-mechanical analysis is performed to simulate the accumulation of excess PWP both in soi and interface elements. During the entire finite element analysis, the water table is maintained at the ground surface by imposing a null-value of PWP at 0.4m from the mud-line. PWP are also fixed in the bottom of the mesh in order to simulate the thin layer of sand that exists in the centrifuge model. In addition, PWP are fixed throughout the pile, supposing that the foundation is totally impermeable.

The clay and the pile are modelled using hydro-mechanical volumetric 8nodes finite elements. The soil-pile interaction is modelled by hydro-mechanical zero-thickness interface elements at the soil-pile contact surface. Relative shear displacement as well as gapping at the interface are therefore allowed. In addi- Because of the relatively large displacements, large deformations are employed, adopting the updated Lagrangian configuration, i.e. the deformed geometry of each element is updated after each loading increment. By doing so, the significant changes in the position and shape of soil elements are captured (see Chapter 2).

Pile loading and drainage conditions

Folllowing the centrifuge test loading conditions, a horizontal displacement is applied to the pile head at a distance of 2m from the mud-line, with a constant velocity of 20mm/s. According to [START_REF] Oh | A review of foundations of offshore wind energy convertors: Current status and future perspectives[END_REF], the selected loading rate represents the rate of the wind field measured in the North Sea, knowing that the total duration of a storm is between 2-10 seconds, and an extreme event can mobilize up to 20% of the pile diameter.

On the other side, according to the soil consolidation parameters, summarized in Table 5.7, the selected loading rate ensures undrained conditions, since the normalized velocity V = 6.10 4 is bigger than the threshold limiting the undrained condition, i.e. V > 20. The pile response is hence independent from the loading rate.

136CHAPTER 5. MONOPILES UNDER MONOTONIC LATERAL LOADING

Soil parameters

The soil matrix is modelled using the SCLAY-1SB law. Only four out of the 13 constants of the soil model are provided directly from the experiments, namely κ, λ, e and OCR. Moreover, constants related to the cohesion and the state of structure are deactivated since the Kaolin clay had been reconstituted. Anisotropy is also deactivated in order to reduce the number of parameters. Consequently, two more constants have to be selected, namely the friction angle and the constant H 0 controlling the bounding surface. The calibration procedure of the soil parameters is detailed hereinafter.

• OCR profile

The overconsolidation ratio is provided by the experimental data, as shown in Figure 5.19 (b). In the finite element model, the OCR is considered constant for the second and the third layer and equal to the experimental OCR at the middle of each layer, i.e OCR = 2.2, and OCR = 2.1 in the second and third layer, respectively. In the upper layer however, the preconsolidation pressure is considered constant to the provided experimental value instead of the OCR (σ c = 55kP a, as presented in Fire 5.19 (a)). By doing so, a best-fitting of the experimental profile of OCR is guaranteed.

• Friction angle

The friction angle is selected from the results of the CPT test provided by [START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF] using the simplified analogy between the CPT and an axially loaded pile. The limit capacity of the CPT q c is interpreted as the limit value of the axial bearing capacity of a pile having the same diameter as the CPT cone. Accordingly, the value of q c is expressed as follows: 5.35) where N q , N g and N c are the axial bearing capacity factors, depending only on the friction angle. B c is the diameter of the CPT, γ the effective specific density of the soil and c the cohesion. Knowing that the soil is reconstituted and the diameter of the CPT is negligible, the second and the third terms of the Equation 5.35 are hence zeroed. The value of N q is expressed using the Meyerhof model for deep foundations. Accordingly, the limit resistance q c is rewritten:

q c = N q .σ v0 + N g .γ B c 2 + N c .c (
q c = exp (7 tan ϕ) q 0 (z) (5.36)
where ϕ is the friction angle of soil at depth z corresponding to plane strain mode of shearing, i.e sin 3θ 0. The friction angle of soil is finally retrieved:

ϕ = tan -1 1 7 ln q c (z) σ v0 (z) (5.37)
Using Equation 5.37, the friction angle is determined all along the pile diameter.

In order to verify this procedure, the undrained shear strength is calculated under plane strain conditions as function of the friction angle, the OCR, and the consolidation parameters (κ and λ) using equation 5.27. The resulting profile is then compared to the measured profile from CPT tests provided by [START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF] in figure 5.21. A good agreement is observed between the measured and the calculated profiles. The maximum difference reaches 7% at 15m depth, which is largely accepted. Hence, an average value of φ = 10 • is considered for the finite element analysis. Neither the profile of q c or the value of s u can be employed for selecting a well-suited value of the constant H 0 . Hence, calibration of this constant was made directly on the pile global response. The selected value is given in Table 5.7.

Table 5.7 summarizes the soil parameters used for the finite element analysis for the different soil layers. The same soil parameters are selected for the three layers, given that they had been prepared with the same Kaolin clay. Differences exist only for the upper layer, where a cohesion of 5kPa is considered to avoid numerical problems, in addition to the elastic parameter κ which has been increased to better reproduce the force-displacement curve at the pile head (Figure 5.22).

Interface parameters

The soil-pile contact is considered fully rough, i.e the soil-pile friction angle is taken equal to two third of the soil friction angle, and the cohesion is considered null. The permeability of the contact surface is assumed equal to the soil permeability (k = 10 -9 m 2 ).

Results analysis

In the following, results of the finite element modelling are compared to the experimental results of the centrifuge tests. All results refer to the prototype model, unless otherwise mentioned.

Global response of the pile

The global force-displacement curve at the pile head (2m from the mudline) has been used during the calibration procedure for selecting H 0 . Figure 5.22 shows the computed and experimental curves. The computed forcedisplacement curve shows a very good agreement with the experimental results, the numerical curve fits almost perfectly the experimental one. At the end of loading, both numerical and experimental curves show an ultimate lateral capacity of 330kN , which is considered as the ultimate lateral capacity of the pile, since the corresponding displacement reaches 20% of the pile diameter [START_REF] Broms | Lateral resistance of piles in cohesive soils[END_REF].

The calibrated finite element model is able to reproduce the global response of the pile. This response is the resultant of the redistribution of the applied force in the soil surrounding the pile. It is interesting therefore to observe the evolution of the soil reaction along the pile. In particular, p-y curves as well as the profile of soil reaction are analysed hereinafter.

p-y curves

The procedure to determine the soil reaction curves from 3D finite element results is implemented in post-processing routines as follows. The computed total normal and vertical shear stresses in the interface elements are integrated around the pile circumference to determine the local lateral load p. These data are related to the local pile lateral displacement y. Consequently, the p-y curves are elaborated at each integration point of the interface element along the pile. Figures 5.23 (a) and (b) compare the numerical to the experimental p-y curves. Although the comparison is qualitative, several conclusions can be made. First of all, the experimental p-y curves show a smooth peak for a depth higher than 2m, while the computed curves have an increasing shape. In addition, the soil capacity at the ground surface seems to be underestimated by the numerical model. At the upper depth (higher than 1m), the computed soil reaction tends toward a flat curve of about 15kP a, whereas the experimental curves are continually increasing. This behaviour is very pronounced for the p-y curve at the ground surface, for which the soil reaction reaches 150kP a.

The analysis of these results is complicated since the p-y curves are not plotted for the same displacement of the pile head. It is more interesting hence to plot, for a given displacement of the pile head, the distribution of the lateral displacement and the soil reaction along the pile. Accordingly, Figure 5.24 (a) compares the computed and experimental profiles of lateral displacements for a pile head displacement of 0.15%B. It can be observed that the numerical lateral displacements fit perfectly the experimental results. This is particularly true for depths higher than 12m, while at the top of the foundation slight differences are observed; a null value of lateral displacement is provided experimentally, while the numerical results predict a displacement of -0.42%B. These differences can be attributed to the procedure of plotting the experimental p-y curves. According to [START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF], the lateral displacements along the pile were computed by double integration of the measured moments, while the distribution of soil reactions was estimated by double derivation of the measured moments. This experimental procedure requires two boundary values at the ends of the foundations. Therefore, the pile bottom was considered fixed to extract the p and y profiles.

The numerical soil reaction profile is compared to the experimental results in Figure 5.24 (b). Qualitatively, the numerical profile of the soil reaction presents a good agreement with the experimental curve. The maximum and minimum values of the soil capacity are well predicted by the finite element model. Furthermore, the evolution of the soil capacity with depth is also well captured, except for the upper layer (until 2m from the mud-line). As previously observed from the p-y curves, the numerical model underestimates the soil capacity at the ground surface, compared to the experimental results. This could be explained by the undrained shear strength of soil at the surface. Even if the soil is highly overconsolidated at the surface, s u remains small, because the vertical effective stress is almost null at the surface. This gives rise to very small values of s u and then p u at the ground surface, especially when the cohesion is equal to zero.

Moments along the pile

The evolution of bending moments are one of the important ingredient for the design methods. In particular, the value of the maximum moment and its corresponding depth play an important role in engineering practice, since piles, and particularly slender piles, resist to lateral loads by bending. Therefore, it is interesting to analyse the evolution of the moments with depth for different values of lateral reaction at the pile head. The procedure to calculate the bending moments along the pile is explained hereinafter. The bending moment is related to the vertical stresses in the pile element by the mean of the classical equation of the Euler-Bernoulli beam, given for a hollow circular section as follows:

M (z) = σ max (z) .I p B /2 (5.38)
where I p is the inertia moment of the section, which is equal to:

I p = π 64 B 4 -(B -2.e) 4
for a hollow and circular section. σ max is the maximal value of the vertical total stress in the pile. The value of σ max is computed considering a linear interpolation of the vertical stress in the integration points of the pile finite 142CHAPTER 5. MONOPILES UNDER MONOTONIC LATERAL LOADING elements, given that the pile behaves linearly and the vertical stresses are null at the pile center. Figure 5.25 compares computed and experimental pile bending moment profiles for different values of horizontal load. The agreement between the experimental and numerical curves is very satisfactory. The model predicts well the evolution of the moment profile with the lateral reaction. It also gives a good estimation of the maximal moments and their corresponding depths. However, the numerical calculations slightly overestimate the moment for deep depths, particularly for depths higher than 8m. It is unclear why the overestimation is limited to high depths. One of the reasons could be related to the mesh refinement, given that the finite element mesh is less refined in the bottom of the pile compared to the upper depths. It can be concluded that the finite element model is capable of reproducing main results of the centrifuge tests and this with a very satisfactory precision. The numerical global pile response and the moment, displacement and soil reaction profiles show a good agreement with the experimental curves.

However, the study can not be concluded before analysing the deformability aspects of the soil surrounding the pile, as well as the behaviour of the soil-pile interface in relation with the soil-pore pressure interaction, which is one of the main advantages of the adopted finite element model.

Soil deformations

As previously mentioned, two failure mechanisms prevail in the soil surrounding a laterally loaded pile. At high depths, the soil fails in a flow-around

3B 2B B 2B
Figure 5.26: Velocity field in the soil at the end of loading mechanism, i.e the soil flows in a non-gapping plane strain mechanism. Near the ground surface, a conical wedge extending to the surface can be observed. Figure 5.26 depicts the velocity field in the soil at the end of loading, i.e for a head pile displacement y = 20%B. As expected, both conical wedge and flowaround mechanisms can be seen. The conical wedge is almost symmetrical around the pile and extends horizontally to a distance of about 1.2B and vertically until three times of the pile diameter, which corresponds to almost 17% of the pile length. Obviously, the wedge dimensions depends of the amount of the pile head displacement. However, the wedge does not extend to high depths in the slender piles, since they are often dominated by the flow-around mechanism.

Figure 5.27 shows the distribution of the volumetric p v and deviatoric p d plastic strains in the soil surrounding the pile at shallow depths, plotted in a deformed configuration of the soil. At the ground surface, a vertical displacement of the soil can be deduced in the positive side of the pile, while the soil settles in the negative side to follow the displacement of the foundation due to the lateral load. This behaviour means that both the passive and active wedge are mobilised around the pile. Hence, soil in the negative and positive sides resist to the lateral pile movement. Expectedly, the shear plastic strains prevail Figure 5.27: Volumetric p v and deviatoric p d plastic strains in the soil surrounding the pile compared to the volumetric plastic strains. The former are concentrated at shallow depths, smaller than 6m, with maximum values around the pile section at the ground surface, particularly in two areas: in the active wedge, where the soil is subjected to the maximum extension, and in the passive wedge, at an angle of about 45 • from the axis of symmetry. The plastic volumetric strains are also concentrated at the ground surface in the soil surrounding the pile. They are almost symmetrical with respect to the pile center, so the condition of constant volume is verified.

Soil-Pile interaction

Figure 5.28 shows the distribution of the effective normal p n and shear τ stresses along the soil-pile interface. The figure depicts the most important mechanical aspects of soil-pile interaction, namely, the opening of the gap between the soil and the pile and the normal friction along the pile. As previously explained in this chapter, these elements are of the uttermost important when dealing with the problem of laterally loaded piles, and in particular for the estimation of the bearing capacity factor.

Obviously, the normal and shear stresses have a similar distribution. They are almost proportional, with a ratio of about 0.1, corresponding to the imposed friction angle coefficient tan ϕ int = 0.1. Furthermore, it can be observed that the distribution of the stresses around the pile varies along the depth. At high depth, values of p n and τ are maximum in the positive side of the pile. At shallow depths, higher values are attained in the negative side. This obser- vation is related to the deformed shape of the pile. Indeed, when a horizontal load is applied at the pile head, at shallow depth, the pile moves horizontally in the same direction as the applied load, whereas a negative displacement can be detected at the bottom of the pile. Therefore, the distribution of the stresses is unversed at high depth compared to those near the ground surface. Finally, it is important to observe that null values of p n are limited to a very reduced area in the positive side at the ground surface, which corresponds to a gap opening at the soil-pile interface. Similar observations have been mentioned by [START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF] while observing the soil surface at the end of the experiment. However, it should be kept in mind that the gapping behaviour depends on the amount of the applied loading. For instance, [START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF] observed the forming of a gap when the pile was loaded up to 1D. In this case, the gap extended horizontally to a distance of about 2 times of the pile diameter. The gapping behaviour is mainly related to the solid-pore water pressure interaction at the soil-pile interface. Figure 5.29 (a) presents the distribution of the PWP around the pile for different sections. In this figure, PWP are normalised with their corresponding initial values. Negative pore pressures are observed in depth smaller than 6m, while slight variations of the PWP are observed for higher depths. Suction prevails in the positive side of the pile and its variation increases when getting closer to the ground surface (PWP are ten 146CHAPTER 5. MONOPILES UNDER MONOTONIC LATERAL LOADING times higher than their initial values at the ground surface). A rigid monopile with a small value of L/B ratio is treated hereafter. More specifically, the following section presents a 3D FEA of a rigid pile of 4m diameter and of 16m length. In order to compare the behaviour with the previous studied rigid slender monopile, all the other parameters are kept the same.

Finite element model

A monopile of B = 4m and L = 16m, embedded in the Kaolin clay, whose parameters are presented in the Table 5.7, is considered. Figure 5.30 illustrates the 3D mesh used for the simulations. Both meshing and boundary conditions are the same to those considered in the previous analysis. As for the slender pile, a horizontal displacement is applied to the pile head 148CHAPTER 5. MONOPILES UNDER MONOTONIC LATERAL LOADING at a distance of 2m from the ground surface, with a constant rate of 20mm/s, which respects the undrained conditions.

5.2 Response of the short pile to a monotonic lateral loading • Given that the rotation center is above the pile base, a 'kick-back' behaviour takes place, accompanied by fairly large displacements mobilised at the pile base.

Theses observations are commonly accepted for short rigid piles, see for example Byrne et al. (2015b) who carried out full-scale tests on laterally loaded short monopiles. .33 shows the bending moment distributions along the pile depth. Compared to Figure 5.25, the maximum moments in the short pile are significantly larger that those in the slender pile. Moreover, the depth corresponding to the maximum moment is deeper that of the slender pile (0.45L for the former and 0.35L for the latter). Therefore, the pile diameter affects not only the maximum moment, but also its corresponding depth.

Conclusions

This chapter aims at describing the salient features of laterally loaded single piles of the OWT sector. A fully hydro-mechanical transient analysis of this type of foundation was carried out.

The chapter opens with a brief state of the art summarizing the main criteria affecting the lateral capacity of piles. In particular, the main failure mechanisms in the soil and the various methods for determining the lateral bearing capacity have been presented. Furthermore, attention has been paid to the p-y curve spring model, by presenting the most popular p-y curves for soft as well as stiff clays.

The second section presents a 2D finite element analysis of a laterally loaded pile embedded in a natural clay. The objective of this section is to present a best-estimate p-y curve based on a fully coupled hydromechanical analysis, where the plasticity of soil is modelled by the SCLAY-1SB model that captures anisotropy and the destructuration of natural clays. The 2D analysis has led to several findings:

• The effect of the loading rate was studied to define the velocity that respects the undrained conditions. It was found that the loading rate affects the initial modulus of the p-y curves; the curves are stiffer for higher loading rates. Beyond a certain limit, the p-y curves remain unchanged, which corresponds to the undrained domain.

• The drainage condition is related to the variation of the PWP at the interface, as well as the opening of the gap in the interface. No-gapping behaviour is observed when the undrained conditions are reached. Vari-ation of the PWP reaches maximal values with negative variation in the positive side of the pile and positive variation in the negative side.

• By computing the undrained normalised p-y curve, the effect of anisotropy and structure of the clay have been studied. It was concluded that the anisotropy of soil or of the stress do not affect the value of N p . However, the shape of the p-y curve is strongly affected since these features influence the value of the undrained shear strength. When the degradation of cohesion and structure is considered, the computed p-y curve presents a peak followed by a softened curve, representing well the stiff clay reaction curves under undrained conditions.

The third section presents a 3D finite element analysis of a centrifuge test carried out in IFSTTAR laboratory, within the framework of the PN-ANR SOLCYP project. The geometry of the test consists of a 0.9m diameter and 16 embedded depth pile with a hollow section of 15mm, loaded laterally at a distance of 2m from the mud-line.The analysis deals merely with monotonic loading and takes into account the plasticity of soil, modelled by the SCLAY-1SB law, and the plasticity of the steel-soil interface, defined by a Mohr-Coulomb criterion. A fully hydro-mechanical approach has been considered in both soil and soil-pile interface, where both flows inside and through the interface are taken into account.

The analysis begins with the calibration of the model parameters. The soil parameters provided by [START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF] were well respected, except in the upper layer of soil, where the elastic parameter κ has been enlarged. The numerical results show generally a very good agreement with the experimental curves. In particular, the global displacement-force curve presents almost a perfect fitting with the experimental curve.

The calibration was followed with an analysis o the soil and soil pile interaction, In particular, the observation of the soil strains around the pile reveals two failure mechanisms, namely the conical wedge and the flow-around mechanism. Moreover, the analysis of the distribution of the effective stresses in the interface demonstrated that both passive and active wedges were mobilised while the gap opening is limited to a narrow surface at the mud-line. Finally, the pore water pressure in the interface have been studied and the presence of suction at shallow depth in the positive side of the foundation has been related to the non-gapping behaviour.

The chapter ends with a numerical analysis of a rigid short pile laterally loaded. The analysis shows the differences of the slender and rigid piles with respect to the pile capacity, strains and the bending moment and soil reaction profiles.

Introduction

Soil is generally subjected to complex loading conditions that cannot be always simplified to monotonic forces. Considering for example offshore foundations, they are subjected to different types of cyclic loading due to harsh environmental conditions, such as wind load, sea waves and earthquake excitations. These significant cyclic loads can lead to permanent deformations of soil surrounding the pile and may cause significant accumulated displacements. Neglecting the permanent displacements and the plasticity/damage may lead to catastrophic failure of the pile-supported structure. For instance, [START_REF] Lin | Damage of piles caused by lateral spreading-back study of three cases[END_REF] performed a back study of three well documented historical cases caused by ground lateral spreading and concluded that the piles failure was mainly due to a significant bending moment and large head displacement in addition to pile buckling, that were underestimated in the foundation design.

Hence, the cyclic behaviour of soil occupies an important place in the field of constitutive modelling. Various concepts have been proposed to replicate important features of the soil behaviour under cyclic loading. Remarkable examples include multi-surface plasticity with kinematic hardening [START_REF] Prevost | Plasticity theory for soil stress-strain behavior[END_REF], hypoplastic models [START_REF] Mašín | A hypoplastic constitutive model for clays[END_REF], subloading surface plasticity [START_REF] Hashiguchi | Elastoplastic constitutive equation of soils with the subloading surface and the rotational hardening[END_REF] and bounding surface plasticity Dafalias (1986b). This manuscript will focus on the last framework, due to its mathematical simplicity and its significant success in representing the cyclic behaviour of soils.

The SCLAY-1SB model, presented in Chapter 4, incorporates the bounding surface plasticity concept in the original SCLAY-1S model Wheeler et al. (2003b) to enhance the prediction of the behaviour of heavily overconsolidated clays. The SCLAY-1SB model has proven capable of modelling different features of clays, such as initial and induced anisotropies, destructuration and softening of sensitive soils and the smooth peak behaviour of heavily overconsolidated clays. However, it is limited to monotonic loading.

This chapter presents an extension of the SCLAY-1SB model for cyclic loading. Three key aspects are added. First, plasticity is modelled at early stages of unloading/reloading cycles by introducing a moving projection center with an appropriate kinematic law. Second, the prediction of both ratcheting (increase of cyclic strain) and shakedown (decrease of cyclic strain, known also as soil densification) during cycling is done via a unique function in the plastic modulus. Last, a non-associated flow rule is adopted to predict the 'butterfly' shape that characterizes the undrained cyclic stress path of clayey soils.

To valide the performances of this extended final version of SCLAY-1SB, the model is calibrated against laboratory tests for Georgia kaolin clay [START_REF] Sheu | Modelling of stress-strain-strength behavior of a clay under cyclic loading (soil dynamlics)[END_REF], tested under different stress histories. Then, the model is validated by simulating in-situ self-boring pressuremeter tests using a 2D-axisymmetric finite element model. The tests data refer to an experimental investigation carried out by [START_REF] Reiffsteck | Utilisation des essais d'expansion cyclique pour définir des modules élastiques en petites déformations[END_REF] in the framework of the SOLCYP project [START_REF] Puech | Design of Piles Under Cyclic Loading: SOLCYP Recommendations[END_REF].

Model formulation under axisymmetric conditions

For the sake of simplicity, the constitutive relations are formulated hereafter with reference to triaxial conditions. Accordingly, the model equations are written in terms of mean effective stress p = (σ a + 2σ r ) /3, and deviatoric stress q = σ a -σ r . Their work-conjugate counterparts are the volumetric strain ε v = ε a + 2ε r and deviatoric strain, ε d = 2 (ε a -ε r ) /3, respectively. Subscripts a and r denote the axial and radial components of stress and strain, while v and q denote volumetric and deviatoric terms, respectively. All the stress variables are considered as effective and compression follows the positive convention both for stress and strain measures.

Intrinsic formulation for reconstituted soil

Since the intrinsic clay behaviour is the focus of this work, the destructuration mechanism of the SCLAY-1SB model is not activated in this section. Only the intrinsic model is formulated in the following. Details regarding the damage of the clay structure are presented and discussed later on in the chapter.

Base-line model: SCLAY-1SB

A detailed formulation of the SCLAY-1SB model for monotonic loading is presented in Chapter 4. The key equations related to the cyclic behaviour are hereafter the only recalled. Figure 6.1 presents the bounding and loading surfaces in triaxial plane and multi-stress space.

The size of the bounding surface is determined by p m , which grows or shrinks in proportion to the plastic volume change, in the same way as in the MCC model. The internal variable α, governs the rotation and distortion of the bounding surface, and represents the degree of plastic anisotropy. The parameter M denotes the stress ratio at the critical state, depending on the shearing mode (i.e., M = M c for triaxial compression and M = M e for triaxial extension).

The image stresses, denoted (p, q) lie always on the bounding surface and are determined by the radial mapping rule proposed by Dafalias (1986b). As shown in Figure 6.1 (a), a projection center (p c , q c ) is used to map radially the current stress (p, q), to (p, q) on the bounding surface. Mathematically, this

q = q c + b (q -q c ) (6.2)
where b varies from 1 to ∞. These two end-values are attained when the current stress coincides with the image stress (b = 1) or the projection center (b = ∞).

Projection center

The version of the model suitable for monotonic loading presented in Chapter 4 adopts a fixed projection center at the origin of stress. Accordingly, the plasticity is allowed at an early stage of virgin loading. However, reversible strains are overestimated during unloading/reloading cycles compared to the commonly observed behaviour of clays, as illustrated in Figure 6.2 (a).

To overcome this deficiency, a moving projection center is adopted in the following. This choice is commonly employed by the bounding surface family of models, as explained in [START_REF] Khalili | A bounding surface plasticity model for cyclic loading of granular soils[END_REF]. It implies that the projection center position has rest to the current stress at any stress reversal, thus allowing plastic strains to occur even at early stage of the unloading/reloading process (see Figure 6.2 (b)). According to [START_REF] Khalili | A bounding surface plasticity model for cyclic loading of granular soils[END_REF], the criterion used to detect a stress reversal is expressed as: L. σtrial > 0 (6.3) where σtrial is the trial stress rate obtained assuming that the strain rate ε is elastic, and L is the loading direction. Following Dafalias (1986b), the loading direction at the current stress is assumed to be the gradient of the bounding surface at the image stress. It reads:

L = ∂F ∂p + ∂F ∂q (6.4)
In parallel, the size of the loading surface has to be reduced to an infinite small value at the stress reversal state, by equalling the distance ratio b to an infinite value, so the radial mapping rule (Equation 6.2) is still verified.

However, a moving projection center has to be considered with a certain accuracy, knowing that the mapping rule should guarantee the uniqueness of the image stress [START_REF] Dafalias | Saniclay: simple anisotropic clay plasticity model[END_REF]. Consequently, an evolution rule has to be defined to control the movement of the projection center between two stress reversals. Various evolution rules have been proposed in the past, one of the simplest rules, proposed by Dafalias (1986b), consists in considering that the increment of the mean component p c of the projection center is proportional to the preconsolidation pressure p m , while the deviatoric component q c is assumed constant between two stress reversals. It reads: ṗc = ṗm p m p c (6.5) qc = 0 (6.6)

The previous evolution rule is only efficient in the case of isotropic soils, when the rotational hardening related to α is deactivated. In the general case, the relative position of the projection center should also be maintained upon rotation and distortion of the bounding surface. To consider this point, [START_REF] Seidalinov | Bounding surface saniclay plasticity model for cyclic clay behavior[END_REF] proposed an evolution rule by assuming that the relative location of the projection center with respect to the bounding surface is constant. It is expressed as:

ṗc = p c p m ṗm (6.7) qc = q c p m ṗm + p c -X p c (p m -p c ) α [(M 2 -α 2 ) .p c . (p m -p c )] 1/2 α (6.8)
where X is a ratio defined in function of the distance of the projection center from the rotation axis that remain constant during loading. Note that the first and second term in the rate equation of q c are separately derived by assuming that the only changed internal variable is p m or α, respectively. Therefore, the additive superposition in the equation 6.8 does not include the effect of the simultaneous change of the aforementioned internal variables.

A novel evolution law for the projection center is proposed in the following to get over this limitation. The idea is to maintain the projection center within the bounding surface by enforcing the incremental nullity of the bounding function at the projection center σ c , which implies:

Ḟ (p c , q c ) = ∂F ∂q c qc + 1 3 ∂F ∂p c ṗc + ∂F ∂κ κ = 0 (6.9)
where κ represents the set of the internal variables of the model, i.e. κ = [p m , α]. Trivial manipulations of the Equation 6.9 lead to the following equation:

M 2 -α 2 p m . ṗc + ṗm .p c p 2 c + 2 (η c -α) . ηc -2 (η c -α) α + 2α (η c -α) 2
M 2 -α 2 α = 0 (6.10) where η c = qc pc is the ratio of the projection center state and ηc its increment ηc = qc qc -ṗc pc η c .

The first term in the left of the Equation 6.10 is related to the mean component of the projection center, while the second term is related to the deviatoric component. By zeroing the two terms, the following system of equations is obtained:

ṗc = p c ṗm p m (6.11) qc = p c p m ṗm -p c 1 - η c -α M 2 -α 2 : α . α (6.12)
Since equations 6.11 and 6.12 are derived by enforcing that the projection center remains always on or inside the bounding surface, the condition of the uniqueness of the image stress is guaranteed regardless of the changes in size, inclination and distortion of the bounding surface.

Plastic modulus and strain accumulation effect

The plastic modulus proposed by Yin et al. ( 2011) is adopted. It reads:

K p = K p + h.M 4 1 + e λ -κ b.p 3 m 1 - 1 b (6.13)
where K p refers to the plastic modulus at the image stress, and results from the consistency condition applied on the image stress. The second term of equation 6.13, named the distance function, includes three key aspects:

-the term 1 + e /λ -κ is in accordance with the same term appearing in the modulus K p and results from the hardening of F = 0 as a function of the plastic void raio change, as suggested by Dafalias (1986b).

-the dependency on the value of M introduces the effect of the mode of shearing in similar way to the variable M . This feature is necessary for reproducing the path-dependence on the plastic compliance of the material, so different values of the plastic modulus are predicted in function of the Lode angle.

h was considered as a constant on the first version of the model. In order to introduce the effect of strain accumulation on the plastic modulus, h is function of a new variable d, as follows:

h = H 0 1 + d (6.14)
where H 0 is a model constant. The new variable d is introduced to take into account the deformation features of cyclically loaded clays, such as the ratcheting and the shakedown effects. The increment of such a new variable takes the following expression:

ḋ = (a d -a v (1 + d)) εp a (6.15)
where

ε p a = ε p d 2 + ε p v 2 (6.16)
Accumulation of both the volumetric and deviatoric components of the plastic strain is considered, via the variable ε p a , defined as the norm of the plastic strain tensor.

-The contraction behaviour of the clay sample is well captured by the constitutive model, -The mean effective stress at the critical state is largely underestimated, compared to the experiment.

-The shape of the stress path at the critical state is quite different from the experiment. A fixed value of p is predicted at the critical state, while the experiment shows that the mean effective stress varies within a butterfly shape when the current stress state approaches the critical state.

In order to overcome this deficiency, the current work adopts a non-associated flow rule. More particularly, the hybrid flow rule proposed firstly by [START_REF] Whittle | Formulation of mit-e3 constitutive model for overconsolidated clays[END_REF], and introduced to the bounding surface models by [START_REF] Shi | A hybrid plastic flow rule for cyclically loaded clay[END_REF] is employed. Accordingly, the increment of the plastic strain tensor is expressed as:

εp v = λ.R v (6.17) εp q = λ.R q (6.18)
where λ is the plastic multiplier. R v and R q are the volumetric and deviatoric component of the flow rule direction, denoted R. As in the conventional associated flow rule of bounding surface models, the deviatoric component R q is expressed solely in function of the image stress. It reads:

R q = ∂F ∂q = 2p (η -α) (6.19)
The volumetric component R v is defined in function of the image and the current stresses. More particularly, [START_REF] Shi | A hybrid plastic flow rule for cyclically loaded clay[END_REF] assumed that R v is a combination of the plastic flow at the image stress, denoted R v and at the current stress, R c v . Such a combination is expressed as:

R v = R i v . 1 b r + R c v 1 - 1 b r (6.20)
where r is a novel constant that controls the relative contribution of the flow directions, and its allows the image stress flow rule to be recovered as a special case when r = 0. R i v and R c v are derived from the bounding surface equation as follows:

R i v = ∂F ∂p = p (M -η) (6.21) R c v = p (M -s l .η) (6.22)
where the value of s l depends on the loading direction and the current stress ratio:

s l = 1 if (η -α) η ≥ 0 (6.23) s l = -1 if (η -α) η < 0 (6.24)
According to [START_REF] Shi | A hybrid plastic flow rule for cyclically loaded clay[END_REF], the parameter s l is introduced to reproduce the plastic contraction following stress reversals. Figures 6. 4 (a) and (b) compare the plastic flow computed through the conventional flow rule and the hybrid flow rule for the case of an undrained cyclic lading. Figure 6. 4 (a) shows clearly the limitation of the conventional flow rule. In fact, the mean effective stress is almost constant when the current stress is being close to the critical state. The results obtained from the same simulation using the hybrid flow rule are presented in Figure 6. 4 (b). The typical butterfly shape of the stress path is correctly captured and the results of the change from contraction to dilation when approaching to the critical state.

General formulation for structured clays

When dealing with a structured material, the destructuration mechanism of the baseline model should be activated. In other words, the internal variables χ and p t , related respectively to the state of structure and the tensile strength should be considered. Figure 6.5 shows the bounding surface for structured clays where p t is the tensile strength of the material, and p mi the preconsolidation pressure at the intrinsic state, when the structure of the clay is totally damaged. p mi is related to the preconsolidation pressure of the material by the mean of the state parameter χ. Chapter 4 describes in details the hardening rules controlling the state parameters p t and χ. In the following, only changes related to the cyclic features are explained. The equations related to the plastic modulus and the flow rule are unchanged. However, modifications have to be considered regarding the evolution rule of the projection center, as explained hereafter.

Let us rewrite equation 6.9 by substituting the set of states parameters κ = [p m , α] by κ = [p mi , α, p T , χ]:

Ḟ (p c , q c ) = ∂F ∂q c qc + 1 3 ∂F ∂p c ṗc + ∂F ∂p mi ṗmi + ∂F ∂p t ṗt + ∂F ∂χ χ = 0 (6.25)
Development of each term in equation 6.25 leads to the following equation:

M 2 -α 2 p m . ( ṗc + ṗt ) + ṗm . (p c + p t ) (p c + p t ) 2 +2 (η * c -α) . η * c -2 (η * c -α) α+2α (η * c -α) 2 M 2 -α 2 α = 0 (6.26)
where η * c is the stress ratio, written as:

η * c = q c p c + p t (6.27)
Therefore, the following system of equations is retrieved:

ṗc = (p c + p t ) ṗm p m -ṗt (6.28) qc = p c + p t p m ṗm -(p c + p t ) 1 - η * c -α M 2 -α 2 : α . α (6.29)
3 Calibration procedure

For the sake of brevity, only the calibration of parameters related to cyclic features is detailed hereafter, with reference to data on slightly sensitive Georgia Kaolin clay. Table 7.3 summarizes the clay parameters used for simulation.

ν κ λ i M c M e ω ω d µ ξ v ξ d ξ H 0 a d a v r 0
.2 0.036 0.121 0.88 0.87 5 0.5 -0 0 0 1.2 10 0 40 Figures 6.6 and 6.7 present the calibration results of the soil parameters related to the monotonic behaviour. The reader is referred to Chapter 4 for a detailed description of the calibration procedure of the parameters related to monotonic features.

Since the Georgia clay is considered as a slight sensitive material, the state variables related to the structure and tensile strength are taken equal to zero. Figure 6.6 illustrates the calibration of the consolidation parameters κ λ and ν, while parameters related to the critical state and the rotational hardening rules have been calibrated on triaxial tests, as depicted in Figure 6.7, on normally consolidated samples. The constant H 0 controlling the plastic modulus has been calibrated on cyclic triaxial tests because of the absence of monotonic triaxial test on overconoslidated samples (see Figure 6.8).

Parameters related to strain accumulation

The parameters a d and a v can be calibrated through observations of the strain developed during cyclic loading. Figure 6.9 (a) presents the accumulation of axial strain for the first three cycles of an undrained triaxial test at q cyc = 165.5kP a, carried out at p 0 = 345kP a. It can be observed that the rate of strain accumulation is increasing with the number of cycles, thus implying cyclic ratcheting, hence, a v = 0. Among the three simulations, the analysis based on a d = 10 provides the best performance. It is noteworthy that a transition behaviour from ratcheting to shakedown is also conceivable, but all the available experiments with number of cycles smaller than 50 reveal a

Parameter controlling the hybrid flow rule

The parameter r influences the stress path during the undrained cyclic loading. The stress path of a cyclic test on the Georgia clay is therefore used to calibrate r. More particularly, the test at q cyc = 121.4kPa has been used for the calibration. A value of r = 40 is selected, as presented in The simulation results show a satisfactory agreement with the experiments. The mean effective stress at the critical state is well captured, especially for the test at q cyc = 121.4kPa. The butterfly shape of the stress path is also well captured as soon as the stress path migrates close to the critical state. The final value of the mean effective stress is satisfactory in tests q cyc = 136kPa and q cyc = 140.7kPa. Nevertheless, the predicted stress paths for the two tests seem to migrate more rapidly than in the experiment. For both tests, the critical state is reached almost after one cycle, compared to five cycles in the measured stress paths. The measured and computed stress-strain curves are presented in Figure 6.11 for the same three tests on Georgia kaolin. The strain accumulation is well predicted in both extension and compression. Especially, for the test computed with q cyc = 121.4kPa, where data of 22 cycles are available, both the computed and measured stress-strain curves present a larger accumulation in extension. For q cyc = 136kPa, where 15 cycles have been measured, the accumulation of ε a is slightly larger in extension than in compression in both experiments and simulations. It is generally accepted that the in-situ tests offer one of the best means of determining the soil properties. In particular, the pressuremeter test in one of the most promising procedures among in-situ tests that is characterized by welldefined boundary conditions, and provides a great ability for measurement of the initial in-situ horizontal stresses. Furthermore, several design methods of foundations have been based on pressuremeter tests, either directly, by employing soil parameters such as the limit pressure and the deformation modulus, or, indirectly by estimating soil intrinsic parameters.

The simulation of pressuremeter tests have been selected to validate the proposed model. More specifically, the cyclic pressuremeter tests carried out by [START_REF] Reiffsteck | Utilisation des essais d'expansion cyclique pour définir des modules élastiques en petites déformations[END_REF] within the framework of the SOLCYP project [START_REF] Puech | Design of Piles Under Cyclic Loading: SOLCYP Recommendations[END_REF]) are modelled and analysed hereinafter.

Pressuremeter apparatus

The pressuremeter is an expandable tube which is either pushed into the soil (self-boring pressuremeter), or inserted into a pre-bored hole (pre-boring pressuremeter) and inflated under controlled conditions, as shown in Figure 6.12 [START_REF] Robertson | In situ testing and its application to foundation engineering[END_REF]. .12: Simplified scheme presenting the apparatus of the self-boring pressuremeter, from [START_REF] Robertson | In situ testing and its application to foundation engineering[END_REF] Plots of pressure as function of volume increase, referring to as the pressure expansion curve, are obtained from the tests, from which values for soil parameters can be determined. Different approaches have been proposed to analyse the expansion of a cylindrical cavity. [START_REF] Hughes | Pressuremeter tests in sands[END_REF] developed a closed form solution for an elastic perfectly plastic undrained clay. Numerically, the problem is often simulated using plane strain or axisymmetric models [START_REF] Fahey | A finite element study of the pressuremeter test in sand using a nonlinear elastic plastic model[END_REF].

Finite element model desription

Mesh and boundary conditions

In the following analysis, simulations of self-boring pressuremeter (SBP) tests are performed using a 2D axisymmetric finite element model, given the symmetry of loading and the transverse anisotropy of the soil. Figure 6.13 shows the geometry and boundary conditions of the model. The vertical displacements are blocked at the top and bottom of the mesh. The lateral boundaries are stress-defined; the stresses at the external lateral boundary are fixed at the initial in-situ horizontal stress value, and are considered isotropic according to the experimental data, while a cyclic loading is applied on the pressuremeter hole, at the lateral boundary on the left side of the mesh, as illustrated in Figure 6.13.

Drainage conditions

The finite element model takes into account the hydro-mechanical coupling in order to capture the variation of pore water pressures (PWPs) in the soil. The water table is fixed at 2m from the ground surface, as described in the experiments, and pore water pressures are considering free to dissipate within the soil mesh following the Darcy's law. The soil permeability is assumed isotropic and equal to k = 10 -7 m/s, and the loading is applied with a constant rate equal to 10 -2 ml/s.

Cyclic loading of the pressuremeter hole

Four SBP tests are selected for simulation. The loading procedure is carried out in four stages, as illustrated in Figure 6.14. First, the probe is expanded until 0.5% of volumetric strain. Then, three cyclic pressures are applied in five stages with different amplitudes, each stage contains 50 unloading/reloading cycles (see Figure 6.14). The selection of the parameters related to the strain accumulation, namely a d and a v , is made using the results of a pressuremeter test. More particularly, the test carried out at depth of 6m, namely SBP 1 , is selected for the calibration. Accordingly, the first 20 cycles of the test have been simulated for three different values of a d and a v , as shown in Figure 6.16. The first simulation (the red curve in Figure 6.16) has been performed with a set of parameters producing the ratcheting effect (increase of plastic strain with cycles), the second set of parameters (the green curve) simulates the shakedown behaviour, while the third simulation (the black curve) considers a transition effect from ratcheting to shakedown behaviour, as the two values of a d and a v are taken positive. The case where shakedown is produced presents a good agreement with experiments. Consequently a d = 0 and a v = 10 are selected for the following simulations.

ν κ λ i M c M e ω ω d µ ξ v ξ d ξ H 0 a d a v r 0.

Cyclic pressuremeter curves

Simulations of the four cyclic pressuremeter tests are performed using the set of parameters presented in Table 6. 3 and Table 6.4. Figures 6.18 (a) and (b) present the measured and computed hole pressure against injected volume and hole volume in function of the loading time respectively. Simulations are in good agreement with the experiments, except that the virgin loading is slightly underestimated. This can be related to the calibration procedure; in fact all models parameters, except those for the strain accumulation, have been calibrated on laboratory tests. However, stress conditions, dimension of samples, and so on, are different from in-situ tests conditions.

Conclusion

This chapter presents the modifications introduced to the SCLAY-1SB model to accurately model the cyclic features of clay. In particular, three modifications have been proposed. First, a moving projection center with a novel evolution rule has been developed in order to predict the plasticity in an early stage of unloading/reloading cycles, while respecting the uniqueness condition of the image stress. Second, a new variable related to the strain accumulation has been introduced in the plastic modulus. The results show that both degradation of the plastic modulus and densification are well reproduced as well as the transition behaviour from degradation to densification. Last, the hybrid flow rule of [START_REF] Shi | A hybrid plastic flow rule for cyclically loaded clay[END_REF] has been adopted to accurately capture the mean effective stress at the critical state. The formulation of the model has been presented in the triaxial plane for both reconstituted and structured materials. The calibration procedure of the new parameters has been described and applied to the reconstituted Georgia kaolin clay. The calibration results show good agreement with the experiments. However, the model has shown its limitation to reproduce the stress path before reaching the critical state.

Finally, the model validation has been carried out on a set of cyclic selfboring pressuremeter tests, carried out on the in-situ and intact Merville clay. The results revealed the densification behaviour of the clay under cyclic loading, and the model has shown its capacity to predict the strain accumulation behaviour of in-situ soils.

Introduction

Offshore piles are subjected to cyclic lateral loads due to harsh environmental conditions, such as wind loads, sea waves and earthquakes. These tough conditions lead to permanent deformations in the soil surrounding the pile, causing significant accumulation of pile head displacements and rotations, which have to be limited to avoid exceeding the serviceability limit of the wind turbine. According to [START_REF] Achmus | Behavior of monopile foundations under cyclic lateral load[END_REF], the permanent rotation of a monopile at the ground surface should not exceed 0.5 • . However, the estimation of the accumulated displacement is considered one of the most challenging tasks for engineers.

Therefore, in addition to the design of the ultimate lateral capacity of monopiles, fatigue design is a substantial aspect for offshore structures [START_REF] Wichtmann | Prediction of longterm deformations for monopile foundations of offshore wind power plants[END_REF], [START_REF] Zhang | Centrifuge modeling of the cyclic lateral response of a rigid pile in soft clay[END_REF], [START_REF] Yang | Numerical study on accumulated deformation of laterally loaded monopiles used by offshore wind turbine[END_REF],ect). More particular, the effect of cyclic loading has to be considered. Neglecting the intrinsic cyclic features of the soil can lead to catastrophic failure at the foundation level, causing overturning or collapse of the major pile-supported structure of the wind turbine [START_REF] Lin | Damage of piles caused by lateral spreading-back study of three cases[END_REF].

Even though several methods have been proposed to estimate the lateral displacement of monopiles under cyclic lateral loading, none has been generally accepted. For instance, the most common design code, API 2010, recommends the soil resistance degradation method to account for cyclic loading. This method reduces the soil reaction value determined from the monotonic curve. However, neither the number of cycles nor the deformability of soil are considered.

Accordingly, understanding the cyclic behaviour of offshore monopiles is of the great interest. Different experimental programs have been carried out in the last decades, such as [START_REF] Hong | Cyclic lateral response and failure mechanisms of semi-rigid pile in soft clay: centrifuge tests and numerical modelling[END_REF], [START_REF] Puech | Design of Piles Under Cyclic Loading: SOLCYP Recommendations[END_REF][START_REF] Puech | Design of Piles Under Cyclic Loading: SOLCYP Recommendations[END_REF][START_REF] Liao | Response of flexible monopile in marine clay under cyclic lateral load[END_REF], among others. In particular, the SOLCYP project [START_REF] Puech | Design of Piles Under Cyclic Loading: SOLCYP Recommendations[END_REF] presented new recommendations based on a series of centrifuge and full-scaled tests on instrumented piles. However, this research project focused on flexible and semi-flexible piles, with no specific attention to OWT monopiles. On the numerical level, two approaches are often used to take into account cyclic modelling, semi-empirical approaches to describe the lateral deformation response of a pile with respect to the number of load cycles and the classic finite element modelling, where the soil is usually presented with a sophisticated material law, and the cyclic loads are simulated in an incremental way.

In this chapter, all the previously presented developments are integrated into an application of a monopile foundation under cyclic lateral loading. The cyclic behaviour of the OWT monopiles is modelled using the implicit finite element method, where the cyclic version of the SCLAY-1SB model is employed to simulate the soil behaviour, and a fully coupled hydromechanical approach is considered both in soil and interface elements.

The chapter is divided into two main parts. First, a 2D-FEA of single pile is modelled and analysed under cyclic lateral loading. Second, a 3D-FEA of a single monopile of 4m of diameter and 16m of length is presented and analysed.

2D-FEA of a pile

A 2D finite element modelling of an offshore pile embedded in clayey soil and submitted to lateral and cyclic loading is presented. The aim is to understand the effect of strain accumulation on the response and to compare the soil deformability around a section of the laterally loaded pile with the soil deformability around a cylindrical expansion cavity. In particular, the cyclic pressuremeter test (SDB 1 performed at 12m of depth) simulated in the previous chapter is compared to the cyclic response of the pile under lateral loading. For this reason, the soil parameters, the geometry and depth of the pile section, as well as the pile loading are considered similar to the aforementioned test.

Finite element model description

2D FE model and soil parameters

As for the pressuremeter test SBD 1 , a pile section at 12m of depth is considered hereinafter. The finite element model is similar to the model presented in chapter 5 (see Figure 5.10). The mesh, the pile geometry and the boundary conditions are identical. The Merville clay, whose parameters are summarized in Table 7.1 and Table 7.2 is adopted. The marine Merville clay is known as a stiff and intact clay. Moreover, as demonstrated in Chapter 6, the Merville clay has a densified behaviour, i.e, the soil becomes stiffer when cyclically loaded, which is the reason why the soil parameters a d = 0 and a v = 10 are selected.

Pile loading

The cyclic loading is a force-controlled signal, applied at the pile section in two phases. First, a monotonic loading is applied until reaching a mean load, denoted H m , that corresponds to lateral displacement of about 3% of the pile diameter. Afterwards, a one-way signal of 50 cycles is applied to the pile section, with a constant amplitude equal to 30% of the mean load, i.e. H c = 0.3H m .

The amount and amplitude of the cyclic signal are deduced from the loading of the cylindrical cavity carried out during the pressuremeter test. In particular, only the first stage of cyclic signal has been considered.

Drainage conditions

Undrained conditions are assumed all along the analysis, by applying a very high loading rate (V = v /B.cv > 20, B = 1m the pile diameter and c v = 10 -7 s -1 the coefficient of vertical consolidation of the soil), which is obviously more representative for offshore piles embedded in clay.

Effect of strain accumulation

Figure 7.1 presents the lateral force-displacement curve of the pile section under cyclic lateral load. The loading curve in the first cycles exhibits the evident non-linearity of soil behaviour. In the following cycles, the unloading and reloading curves are relatively linear and stiffer than those in the first cycles. As shown in Figure 7.1, the secant subgrade modulus of the soil, K s,i , increases with cycles and tends towards the elastic stiffness. This behaviour corresponds to the 'plastic shakedown' of the pile. According to [START_REF] Swane | Shakedown analysis of a laterally loaded pile tested in stiff clay[END_REF], this phenomenon occurs in every pile embedded in an elastoplastic soil, whose properties remain unaltered after a certain number of cycles. When the accumulated permanent deformations 'stabilize', the shakedown occurs.

In general, the shakedown phenomenon is a result of two mechanisms: (a) the cyclic behaviour of the soil due to deterioration or densification of the soil properties with cycles, and/or (b) the soil-pile interaction due to the gaping behaviour at the soil-pile interface and the vertical load-spreading along the pile. The former is well captured by the numerical model as the constitutive model for soil captures the strain accumulation effect; the latter cannot be simulated in the present simulations since the numerical model is bi-dimensional.

The densification/degradation of soil under cyclic loading is controlled by the model parameters a d and a v , as explained in Chapter 6. In order to assess the influence of these parameters, a simulation is carried out with a set of parameters a d = 20 and a v = 0 to generate degradation of the soil plastic modulus. Result are presented in Figure 7.2 in terms of lateral load- In the reference case (a d = 0; a v = 20), the increment of pile displacement decreases and tends to stabilize with cycles, since the soil becomes stiffer due to strain accumulation. On the latter case (a d = 20 and a v = 0), the increment of pile displacement increases rapidly with loading until failure at the 10 th cycle. The pile response becomes softer and no stabilization is observed.

Comparison of the pile section and the pressuremeter test

One obvious application of the pressuremeter test is in the design of laterally loaded piles because of the loading analogy. The similarities between the two problems have been explained by several authors [START_REF] Briaud | Using the pressuremeter curve to design laterally loaded piles[END_REF], [START_REF] Byrne | Prediction of P-Y curves from Pressuremeter Tests and Finite Element Analyses[END_REF] and [START_REF] Robertson | Design of laterally loaded piles using the pressuremeter[END_REF]. According to [START_REF] Briaud | Using the pressuremeter curve to design laterally loaded piles[END_REF], the analogy however is not complete and the pressuremeter curve is not identical to the p-y curve. As shown in Figure 7.3, only the component of pile capacity related to the passive side of the section is similar to the soil behaviour around the pressuremeter cavity.

Undrained conditions Pressuremeter

Risisting soil mass Lateral load

Initial cavity

Expanded cavity

Active side Passive side

Figure 7.3: Analogy between a laterally loaded pile section and the soil cavity expansion, from [START_REF] Briaud | Using the pressuremeter curve to design laterally loaded piles[END_REF] Therefore, most of the pressuremeter-based methods propose a similarity ratio to convert the pressuremeter curve to the p-y curve. For instance, [START_REF] Robertson | Design of laterally loaded piles using the pressuremeter[END_REF] suggested that the pressuremeter curve should be increased by 2.1 to 2.5, depending on the soil properties and the pile depth, to get the correct p-y curve. This was confirmed with finite element analyses, as demonstrated in [START_REF] She | Prediction of P-Y curves from finite element analyses[END_REF] by means of 2D-FEA.

The pressuremeter test may also be useful to evaluate the behaviour of piles under cyclic lateral loading, although this potential has not yet been investigated; only monotonic behaviour has been studied in the past [START_REF] Cosentino | Standardizing the pressuremeter test for determining py curves for laterally loaded piles[END_REF]. To meet this objective, results of the 2D-FEA of the laterally loaded pile under undrained conditions are compared to the simulation of the cyclic self-boring pressuremeter, performed with the same soil parameters and at the same depth (the simulations presented in the section 4). The lateral displacements of the pressuremeter tests and the 2D pile simulations have almost a similar evolution. The 2D-FEA pile displacement is slightly bigger than the radial displacement of the pressuremeter hole, but differences does not exceed 5%. This result is encouraging for investigating a possible application of the cyclic pressuremeter test for evaluation of the accumulated displacement of laterally loaded piles. Based on monotonic loading, the lateral displacement y 1 could be estimated. Then, the evolution of the lateral displacement could be estimated from cyclic pressuremeter tests, performed along the pile depth. Validation needs further research, incorporating cyclic pressuremeter tests in different clays, and also several numerical studies, especially three-dimensional analysis, to take into account the different mechanisms affecting the behaviour of laterally loaded piles.

The 2D-FEA of monopiles under plane strain conditions and cyclic lateral loading is a first step analysis, which is very efficient in terms of computational demands. It is important however to keep in mind that a 2D analysis cannot explain several phenomena in laterally loaded monopiles. For instance, at shallow depth, soil fails in a conical wedge that extends to the surface. Whereas, at a deeper depth where soil failure is no longer influenced by ground surface, soil failure corresponds to a localised flow around mechanism. Thus, a 2D analysis is only valid at deeper depths where the flow around mechanism prevails. Furthermore, for short piles (L/B<5) that represents the majority of the existing OWT monopiles, conical wedge failure dominates, and flow around mechanism may completely disappear. Finally, in the cyclic 2D analysis presented, a stabilization of the monopile behaviour with cycles is observed, which is known as shakedown phenomenon. This phenomenon is related to both the strain accumulation of the soil and the stress-spreading with depth due to the cyclic loads. If the strain accumulation is well captured by the soil model, the stressspreading cannot be captured because of the 2D plane-strain assumption. A 3D modelling is therefore necessary to simulate more accurately the behaviour of OWT monopiles under cyclic lateral loading.

3D-FEA of OWT monopiles

In this section, a tri-dimensional modelling of an OWT monopile (B = 4m and L = 16m) embedded in clay soil is presented. The case of cyclic loading is studied within a fully hydromechanical approach.

Figure 7.5 illustrates the finite element model employed for the analysis. The model is identical to the one presented in Chapter 7 for the analysis of the monotonic behaviour of short monopiles, in terms of mesh, pile geometry, mechanical and hydraulic boundary conditions, as well as drainage conditions.

The pile is loaded at 2m from the ground surface. The hight selected represents the moment eccentricity, M /H where H is the lateral load and M the moment in the pile at the ground surface, that is representative of wind and wave loading on a full scale of wind turbine structure. In practice, values of M /H between 1 and 5 represent wave loading, while higher values, i.e. 5 < M /H ≤ 15, represent wind loading Byrne et al. (2015a). Georgia Kaolin clay is selected as soil material all along the analysis. Soil is modelled using the SCLAY-1SB law, with the set of parameters calibrated on cyclic undrained traixial test (Chapter 7). Table 7.3 summarizes the soil parameters adopted for the simulation. Georgia Kaolin is selected for the analysis due to its cyclic features; as mentionned in Chapter 7, this clay presents a ratcheting behaviour when laterally loaded. Hence, it seems interesting to observe its effect on the global response of the pile when subjected to cyclic loading.

Table 7.3 The shear behaviour at the soil-pile interface is simulated with the Mohr-Coulomb criterion, where the soil-pile contact is considered fully rough, i.e. the soil-pile friction angle is equal to 2/3 of the soil friction angle and cohesion is null.

Monotonic loading response

In order to select and apply a cyclic loading on the monopile, a displacementcontrolled monotonic load is applied on the pile head. Figure 7.6 (a) and (b) show the resulting load-pile head displacement and pile rotation versus the pile moment at the ground level, respectively. The maximum displacement reached is equal to 8.6% of the pile diameter, which is equal to 34.4cm (the code did not converge after this level). The corresponding ultimate lateral pile capacity is thus considered equal to P u = 5M N .

Cyclic loading response

Cyclic loading

The cyclic loading is a force-controlled signal applied in two phases. First, a monotonic loading is applied until reaching the mean load H m , followed by a one-way cyclic signal with a constant amplitude, noted H c . The mean load is taken equal to H m = 0.3P u = 1.5M N , while the cyclic amplitude is equal to H c = 0.4P u = 2M N with a total number of cycles N = 25.

Ideally, simulation should be carried out for thousands cycles, given that the number of cycles due to wind and wave could exceed 10 8 over the lifetime of the structure [START_REF] Achmus | Behavior of monopile foundations under cyclic lateral load[END_REF]). However, this procedure is not feasible herein, since the considered method, namely the implicit finite element method, simulates the behaviour cycle by cycle while considering the complexity of the soil behaviour, which hugely increases the computation time. Furthermore, the accumulation of the numerical error could be significant compared to the increment of pile displacement which generally become smaller with increasing cycles.

Lateral displacement and rotation at the pile head Figure 7.7 presents the lateral force-displacement curve at the pile head and the rotation-moment at the ground surface. The lateral displacement curve exhibits the evident non-linearity of soil behaviour, especially for the first cycles. In the following cycles, the unloading and reloading curves are slightly stiffer than those in the first cycles, and the increment of lateral displacement decreases with increasing cycles. The evolution of the pile rotation with the applied moment shows a similar behaviour. The increment of pile rotation decreases with increasing cycles, which leads to a stabilized behaviour. of the pile, normalized with their value at the first cycle, in function of the number of cycles. Both lateral displacement and rotation are continuously increasing with cycles, tending towards a plateau (however, more cycles are required to reach the plateau). Moreover, the evolution of rotation and lateral displacement is quite similar, but the rate of rotation is slightly faster. It may be concluded that pile rotation prevails that the displacement, when the pile is cyclically loaded . Obviously, the H-Y curve is similar to that retrieved from the 2D analysis shown in Figure 7.1. However, the behaviour is quite different; in the 2D analysis, the shakedown is merely related to the soil behaviour, which has been modelled to be stiffer with increasing number of cycles i.e. densification. In the 3D analysis, the soil stiffness deteriorates with increasing cycles. Accordingly, in this case, the shakedown results from the vertical load distribution along the pile. This point is explained in details hereafter.

Displacement and moment profiles

Figure 7.9 (a) depicts the displacement profiles for different cycles. As previously observed, the increment of pile displacement decreases with cycles. After the first cycle, the lateral displacement at the ground level reaches 2% of the pile diameter, and moves only with almost 0.75% of the pile diameter after 25 cycles.

As a consequence of the resulting rigid motion, the pile rotates around a rotation center positioned at 3m from the pile base, that is 0.75 of the pile diameter. Furthermore, the position of the rotation center remains unchanged for the different cycles. Consequently, a 'kick-back' behaviour takes place, with The moment profiles are presented in Figure 7.9 (b) for different cycles. A slight increase of the maximum moment is observed. Precisely, the maximum moment increases of about 5% between the 1 st and the 25 th cycle. Furthermore, the corresponding depth becomes slightly deeper, varying from 8.8m in the first cycle to 10.2m at the end of loading.

Cyclic p-y curves

As already explained in Chapter 4, the conventional p-y approach can not predict the complex behaviour of short monopiles but it can be helpful to explain some aspects of the cyclic behaviour of short and rigid piles, especially the vertical redistribution of loads along the pile. Accordingly, the p-y curves are extracted for different sections along the pile, as illustrated in Figure 7.10. Figure 7.11 depicts six soil reaction curves at different depths both for monotonic and cyclic loading. The p-y curves are normalized by the corresponding undrained shear strength, which is calculated from soil parameters using Equation 5.27.

In Figure 7.10, it can be observed that the reaction curves at smaller depth than 7.0m present a softening behaviour, while at deeper depths they evolve following an increasing path. Figure 7.11 reveals three different behaviours, depending on the soil depth:

• At shallow depth, less than 3m in this case study, the mobilized reactions decrease under the effect of cycles. Moreover, the envelope curve tends towards a flat plateau. At the end of loading, the lateral capacity ratio is reduced from 3.1 for the monotonic curve to 1.6 for the envelope curve of the cyclic soil reaction.

• For relatively bigger depths, between 3m to 5.5m, the soil reactions are continuously decreasing, and no plateau is observed at the end of loading.

• At deeper depth, higher than 6.1m, the soil reactions increase with cycles, while following the monotonic p-y curves.

This behaviour is common for cyclically loaded monopiles, either slender or short and rigid piles. It is due to the transfer of soil reactions under the effect of cycles from the shallow to deeper depths. With this type of analysis, the reduction coefficient that should be applied to the lateral capacity ratio, usually used in design codes, and the depth limit from which the soil reaction stop decreasing can be deduced. It is noteworthy that these two parameters are related to the properties of cyclic loading, namely, the mean load and the amplitude of cycles, as well as the number of cycles. can be concluded that the rotational flow mechanism replaces the flow around mechanism that prevails in deeper depths for slender piles. These observations are broadly consistent with the experimental observations presented by [START_REF] Hong | Cyclic lateral response and failure mechanisms of semi-rigid pile in soft clay: centrifuge tests and numerical modelling[END_REF] from centrifuge tests on laterally loaded short piles.

Conclusion

This chapter describes the behavior of offshore piles under cyclic lateral loading. A fully hydro-mechanical coupled numerical analysis was carried out.

In the first section, a bi-dimensional finite element model of a laterally pile has been performed. In particular, the effect of strain accumulation due to cyclic loading has been studied. Afterwards, an attempt to relate the accumulated displacement of the pile to the radial displacement of the pressuremeter cavity has been discussed. The analysis show encouraging results. A possible application of the cyclic pressuremeter test for evaluating the accumulated displacement of laterally loaded piles is feasible, but needs further experimental and numerical research.

The second section of the chapter presents a 3D analysis of a short monopile under cyclic lateral loading. The numerical analysis reveals a shakedown behaviour of the pile, both the increments of pile displacement and rotation decrease with increasing cycles. This behaviour is related to the distribution of the load along the pile depth. Indeed, the comparison of the monotonic and cyclic soil reaction curves offer a detailed explanation. In particular, a decrease of soil reaction is observed in shallow depths, while, for deep depths, the soil reactions increase while remaining on the monotonic p-y curves.

Finally, the failure mechanism in the soil near the foundation has been analysed. It results that a wedge mechanism prevails at shallow depths extending to the ground surface, whereas a rotational failure mechanism appears close to the pile base.

This chapter summarises and concludes the main findings of this research work. First of all, this research is fundamentally based on numerical modelling of the offshore foundations using finite element code of LAGAMINE. Accordingly, the basics of the finite elements employed during this work are elaborated.

Research assumptions

Finite element approach

The general approach involves representing the soil (clay) as a saturated porous material. The quasi-static analysis is considered throughout this research work, either for monotonic and cyclic loading condition. In this context, the mechanical equilibrium equation is expressed in the usual differential local form, and the weak form is written in the current configuration, which corresponds to the updated Lagrangian configuration. The hydraulic equilibrium is based on the equation of mass conservation, which is written in the local and global forms in function of Eulerian variables.

The flow of the fluid mass is described by the Darcy's law. The flow is hence considered laminar, and the soil permeability is assumed isotropic. Whereas, the stress acting on the solid skeleton are assumed equal to the effective stress of Terzaghi.

The soil-pile interface has been modelled using the hydro-mechanical contact element that belongs to the zero-thickness family, which takes into account the following:

• Sliding, captured by the Mohr Coulomb criterion;

• Fluid flow parallel to the interface, which is described by the generalized Darcy law;

• Fluid flow from the solids through the interface;

• Storage flow, which ensures the fluid mass conservation, and gives a hydro-mechanical coupling term.

Constitutive model for clay under monotonic and cyclic loading

Along this research work, the constitutive model SCLAY-1SB has been used to predict the soil behaviour. The SCLAY-1SB model, which has been developed within the framework of the PhD thesis, is an elastoplastic constitutive law and belongs to the bounding surface family of models, where the common model for clay SCLAY-1S is taken as the baseline model. As a result, the developed law is able to predict the following features of natural clay:

• Plastic anisotropy,

• Destructuration,

• Smooth transition from elastic to elastoplastic behaviour

• Effect of the third invariant of stress or the Lode angle,

• Strain accumulation effects (Ratcheting, densification, and a transition, behaviour from ratchting to densification,

• the butterfly shape of the undrained cyclic stress path.

The model has been validated on different clay and various experiments, in order to test its performance regarding the different clay behaviour, and different stress paths. In particular, oedometer and triaxial tests have been simulated and compared to experimental data, in addition to undrained cyclic traixial tests and cyclic pressuremeter in-situ experiments. Generally, the model has shown a good agreement with the experimental results.

Contributions and outcomes

This thesis contributes to the body knowlegde by introducing a cyclic constitutive law for clay in the finite element code of LAGAMINE, which has been achieved through sequential steps. First, a deep understanding of the clays behaviour and their intrinsic features has been intiated. Afterwards, an analyses of the existing constitutive laws and theories has been thoroughly analysed, in order to unravel their most important elements and detect their limitation. Finally, a constitutive law has been developed in both the triaxial plane and the multiaxial space, and then validated on various laboratory and in-situ tests carried out on different soils.

The implementation of the constitutive law in the finite element code of LAGAMINE opened up new horizons that are concerned with the capability of reproducing the behaviour of monopiles embedded in clayey soils.

Indeed, studying of the offshore monopiles can be pushed froward. The numerical analyses performed all along this research can be distibuted under three sections:

• The 2D and 3D analysis, where plane strain and tri-dimensional finite element models has been employed;

• Tackling monotonic and cyclic loading of the foundation;

• Analysing and comparing the behaviour of slender and rigid piles.

2D analyses

Loading rate and drainage conditions It has been found that the loading rate affects the initial modulus of the py curve, where the curve is stiffer. When the rate is higher and beyond a certain limits, the p-y curve is unchanged, which corresponds to the undrained domain.

Effect of soil behaviour

By computing the undained normalised p-y curve, the effect of anisotropy and structure of the clay have been studied. It was concluded that neither the anisotropy of soil nor the stress affects the value of Np . However, the shape of the p-y curve is strongly affected since these features influence the value of the undrained shear strength. When the degradation of cohesion and structure are considered, the computed p-y curve presentes a peak, followed by a softened curve, which is a well represention of the stiff clay reaction curves under undrained conditions.

Cyclic response of pile section

The aim of this analysis was to confront the lateral displacement of the pile section and radial displacement of the pressuremeter cavity. The results were promising and call for further research to study the application of the cyclic pressuremeter test on the design of offshore monopiles.

3D analyses

Modelling of centrifuge test

The main improvement was the confrontation of the numerical modelling with laboratory tests . In particular, a centrifuge test of laterally loaded pile has been simulated and compared to the experimental results. This application has been limited to the monotonic loading and slender piles, even though it was helpful to validate the numerical model.

Behaviours of rigid and slender monopiles under monotonic and lateral loading

Second, the behaviours of rigid and slender monopiles have been compared by the mean of 3D modelling. As a result, differences on pile deflection, pile capacity and distribution of load and moments along the pile depth have been observed.

Rigid monopile under cyclic lateral loading

The final numerical appreciation consists of modelling a rigid pile under cyclic loading. The simulation has revealed several observations. First of all, the analysis of the numerical p-y curves has explained the shakedown behaviour of the pile, since pile reactions decrease with cycles for the shallow depths and increase near the pile base. Second, the maximal moment of the pile is slightly increasing with the cyclic loading and its corresponding depth is varying with cycles. Finally, the observation of soil flow around the pile has revealed the existance of two different mechanism: the wedge mechanism that extends to the ground surface, and the rotational mechanism in deeper depth.

Strength and limitations

The specificity of the finite element analyses carried out in this research work is the consideration of a fully hydro-mechanical approach in both soil matrix and interface element, where both flow inside and through the interface sides have been taken into account . Through this approach, it has been possible to achieve the best approximation of the real behaviour of pore-water pressure around the piles. Furthermore, the development of the soil model SCLAY-1S has allowed the consideration of the complexity of soil behaviour. However, the sophistication of soil model has several drawbacks. In fact, the considered finite element approach, namely the implicit method, is restricted to few number of cycles since the soil behaviour is simulated cycle by cycle. Thus, the pile displacement become smaller with the increasing number of cycles, and therefore the numerical error could exceed the pile displacement and may lead to some failures during the analysis process.

Future research

The soil model SCLAY-1SB, as presented in this research, is a rate-independent constitutive law. Therefore, an extension for the constitutive, which counts the rate-dependency, would be of the great interest for future research. Furthermore, the enhancements of the rule that control the strain accumulation is highly recommended. For example, the introduction of the memory surface in the model, or a pseudo-creep concept, could largely enhance the prediction of strain accumualtion due to cyclic loading. Last but not least, numerous questions related to the behaviour of rigid monopiles remain unresolved. The analysis presented in this research represents only a first step of a deep and rich investigation, where the effect of pile geometry and soil behaviour are prone to further analyses. Nowadays, offshore wind energy industry is developing exponentially, due to the significant contribution of the North Sea wind turbines energy production to the total consumed energy in Europe. Given that the EU's energy target is to increase the share of renewable energy by 2020, there is a great potential of the offshore wind energy applications towards this direction. In order to achieve this, the support of scientific research is crucial. Monopiles have been by far the most support structure for offshore turbines, nowadays becoming applicable also for complex site conditions. The main objective of this PhD thesis is to study the soil-foundation interaction problem for offshore wind turbines monopiles embedded in clays. We focus therefore on the numerical modelling of natural clay behavior. We aim to develop a constitutive model for clay soils, which allow developing new p-y curves that could be widely applied in offshore wind turbines monopiles

31

  constitutive modelling . . . . . . . 5 1.4 Finite element modelling of offshore Monopiles . . . 6 2 Research aim and objectives . . . . . . . . . . . . . 7 Research outline . . . . . . . . . . . . . . . . . . . . 8 Motivation and context of the research work 1.1 European energy targets and wind turbines With regards to the European Union's (EU) aim to achieve 20% of Green House Gas reduction, increase of the energy production from renewable sources, and improve energy efficiency, the European Commission draws the roadmap, represented by Directive 2009/28/EC, for research and action plans across Europe.

Figure 1 . 1 :

 11 Figure 1.1: Projected costs from a study done in the early 90's in Vendeby, from EWEA (2017)

Figure 1 . 2 :

 12 Figure 1.2: Distribution of bottom sediments of the North Sea, from Bjerrum(1973a)

Figure 1 . 3 :

 13 Figure 1.3: Behaviour at the soil-pile interface

Figure 2 . 1 :

 21 Figure 2.1: Description of a standard porous medium (REV -Representative Element Volume)

Figure 2 . 3 :

 23 Figure 2.3: Updated Lagrangian configuration

  17

Figure 2 . 4 :

 24 Figure 2.4: Position of the stress state in the principal stress space (a) and in the deviatoric plane (b). σ * 1 ,σ * 2 ,σ * 3 are the projection of the principal stresses on the deviatoric plane.

Figure 2 . 6 :

 26 Figure 2.6: Hydraulic characterization of the interface between two deformable bodies, from Cerfontaine (2014)

Figure 2 . 7 :

 27 Figure 2.7: Types of hardening rules

Figure 2 .

 2 Figure 2.9: Drained triaxial test on London clay[START_REF] Gens | Stress-strain and strength characteristics of a low plasticity clay[END_REF]

Figure 2 .

 2 Figure 2.10: Critical state and critical state line in p -q plane and e -p plane

Figure 2 .

 2 Figure 2.11: Yield surface of Cam-Clay models

Figure 3 . 2 :

 32 Figure 3.2: Basic crystral structural units of phyllosilicates (Mitchell and Soga (2005))

Figure 3 . 4 :

 34 Figure 3.4: Classification of fabric, from[START_REF] Sides | The microstructure of dispersed and flocculated samples of kaolinite, illite, and montmorillonite[END_REF]

Figure 3 . 5 :

 35 Figure 3.5: Volume change of dense and loose soil samples, after Reylonds (1885)

Figure 3 . 6 :

 36 Figure 3.6: Contractancy and dilatancy chart, from[START_REF] Shimizu | Effect of overconsolidation on dilatancy of a cohesive soil[END_REF]

Figure 3 SimulationsFigure 3 . 8 :

 338 Figure 3.7: Stiffness strain curve for reconstitutied clay at different OCR, from Jardine et al.(1984)

Figure 3 . 9 :

 39 Figure 3.9: Failure envelope in the Hvorslev-Modified Cam-Clay model[START_REF] Zienkiewicz | Finite element studies of soils and porous media[END_REF]

Figure 3 .

 3 Figure 3.10: Schematic illustration of the Bounding Surface with the radial mapping rule in the stress space, from[START_REF] Kaliakin | Theoretical aspects of the elastoplasticviscoplastic bounding surface model for cohesive soils[END_REF]

Figure 3 .

 3 Figure 3.11: Variation of undrained shear strength with angle of stress reorientation for different specimens undr compression, from O'Rourke et al.(1976)

  Figure 3.12: Yield surface determined experimentally on St. Alban and St.Louis clays, from[START_REF] Tavenas | Effects of stresses and time on yielding of clays[END_REF]

  Figure 3.13: Triaxial tests under constant stress-ratio[START_REF] Gao | Compression and extension yield of an anisotropically consolidated soil[END_REF]

Figure 3 .

 3 Figure 3.14: Inclined yield surface proposed byDafalias (1986a)

Figure 3 .

 3 Figure 3.15: Undrained stress paths in compression and extension after K 0 consolidation, capturing the uniqueness or lack of CSL, depending on the RHR, from Dafalias and Taiebat (2013) It can be concluded that rate equations of α given as function of εp v (or ∂g ∂p ) do not predict a unique critical state in ln p -e plane. This condition is however not sufficient, since the value of the mean effective stress at the critical state p c is not yet examined. This value should be unique, not dependent on the shearing mode.The p c value can be determined from the yield surface equation (by replacing q c = M p c ), which leads to:

Figure 3 .

 3 Figure 3.16: Sensitivity chart after[START_REF] Leroueil | Discussion on: Correlations between index tests and the properties of remoulded clays[END_REF]

  (a)-(c), which depict several results of triaxial and oedometer tests on intact and reconstituted specimens of different origins.

Figure 3 .

 3 Figure 3.17: Experimetal tests on intact and reconstituted samples,[START_REF] Leroueil | The general and congruent effects of structure in natural soils and weak rocks[END_REF]: (a) one-dimensional compression curves for intact and remoulded samples of Mexican clay; yield curves of intact and remoulded samples of (b) Backebol clay and (c) Atchafalaya clay.

Figure 4 . 1 :

 41 Figure 4.1: Schematic illustration of the intrinsic and bounding surfaces

  surface is illustrated in figure 4.2 in the deviatoric plane. The Van Eekelen criterion is built from the Drucker-Prager criterion, for which m = M c = M e , by introducing a dependence on the Lode angle, with the purpose of matching more closely the Mohr-Coulomb criterion as illustrated in Figure 4.2 (a). Consequently, lower stresses are permitted in extension compared to compression (M e < M c for ϕ e = ϕ c ).

  Figure 4.2 (b) shows the difference between the Mohr-Coulomb (where M is considered constant), Von-Mises and Van-Eerkelen surfaces for a friction angle of 20 • .

Figure 4 . 2 :

 42 Figure 4.2: Friction angle models in the deviatoric plane, after[START_REF] Barnichon | Finite element modelling in structural and petroleum geology[END_REF]

Figure 4

 4

Figure 4 . 3 :

 43 Figure 4.3: Calibration of ν based on K 0 -consolidation test (data of Bothkennar clay from Smith et al. (1992))

3. 3 ¡Figure 4 . 4 :

 344 Figure 4.4: Calibration of constant ν based on K 0 -conosolidation test (data of Bothkennar clay from Smith et al. (1992))

Figure 4 . 5 :

 45 Figure 4.5: Simulation and experiments of drained triaxial tests on London clay

100CHAPTER 4

 4

Figure 4 . 7 :

 47 Figure 4.7: Simulation of drained triaxial tests on Otaniemi clay

Figure 4 .

 4 Figure 4.9 and Table4.7 summarize the data of undrained and drained triaxial tests carried out on intact samples of Bothkennar clay. The samples were first consolidated back to the in-situ stresses (the point A in Figure4.9) along the K 0 -line, which retraces their lightly overconsolidated stress histories. From the

Figure 4 .

 4 Figure 4.8: Simulation of oedometer test on reconstituted and intact samples on Bothkennar clay

Figure 4 . 9 :

 49 Figure 4.9: Stress paths of drained triaxial tests on Bothkennar clay

  (a)-(b), in terms of q -ε d , p-ε v , and ε v -ε d curves, respectively. Numerical results are in good agreement with the experimental curves. The constitutive model captures relatively well the volumetric and the deviatoric shear strains under different stress-ratios.

Figure 4 .

 4 11 compares model simulations with the experimental results of undrained trixial compression and extension on Bothkennar samples. Both

Figure 4 .Figure 4 .

 44 Figure 4.10: Simulation and experiments of drained triaxial test of Bothkennar clay

Figure 5 .

 5 Figure 5.1 shows the current trend of foundations used in the commercial offshore wind sector with respect to the seabed depth. According to[START_REF] Oh | A review of foundations of offshore wind energy convertors: Current status and future perspectives[END_REF], gravity and monopile types of foundation are mainly used in shallow waters, due to their reliability. In transitional water, monopile, caisson and tripod are commonly employed. Whereas for water depth over 30 meters, multi-pod (either multi-caisson or jacket structures) are usually selected due to their lower costs. More recently, a new generation of Monopiles of 4 -10m of diameter, named XL Monopiles, are studied to be used in deep water(Byrne et al. (2015b)).

Figure 5 . 1 :

 51 Figure 5.1: Types of foundations of OWTs, from EWEA (2017)

Figure 5 . 2 :

 52 Figure 5.2: Schematical representation of loads acting on OWT supported on a monopile, after Bhattacharya et al. (2011)

110CHAPTER 5 .Figure 5 . 3 :

 553 Figure 5.3: Failure mechanisms of laterally loaded pile

Figure 5 . 5 :

 55 Figure 5.5: Limit equilibrium analysis of a passive soil wedge by Reese, after Jeanjean et al. (2017)

Figure 5 . 6 :

 56 Figure 5.6: Winkler approach (p-y model) for a laterally loaded pile, after Jeanjean et al. (2017)

Figure 5 . 7 :

 57 Figure 5.7: Undrained p-y curve proposed by Reese et al. (1975)] for stiff clays

118CHAPTER 5 .

 5 Figure 5.9 presents the idealized results of an undrained triaxial compression test reproduced with the SCLAY-1SB model for the case of an overconsolidated, anisotropic and intact clay sample. The specimen's state is defined with: the stress invariants (p ,q,θ), the state parameters (e,p mi ,χ, p t ) and the consolidation parameters κ and λ i .

Figure 5 . 9 :

 59 Figure 5.9: Undrained triaxial compression test reproduced with the SCLAY-1SB model for an overconsolidated, anisotropic and intact clay sample; yield envelope and critical state line

Figure 5 .

 5 Figure 5.10: 2D finite element model

Figure 5 .

 5 Figure 5.11: Interface finite element

  Figure 5.12: p-y curves computed wth different loading rates

Figure 5 .

 5 Figure 5.13: Normal pressure and PWPs around the pile section at the end of loading (y = 20%B) for different loading rates

128CHAPTER 5 Figure 5 .

 55 Figure 5.14: PWP variation in the soil when (a) y = 0.01B, (b) y = 0.05B, and (c) y = 0.15B

Figure 5 .

 5 Figure 5.15: Effect of (a) stress anisotropy, (b) soil anisotropy on the p-y curve

Figure

  Figure 5.15 (b) shows the effect of soil anisotropy on the soil reaction curves. In this case, the initial stress conditions are considered isotropic while the rotational hardening rule (RHR) is activated to capture the anisotropy induced by the pile loading. The RHR is controlled by the parameter A that predicts the evolution of the internal variable α. The resulting normalized p-y curves have a closed shape and almost the same value of the bearing capacity factor. Furthermore, a smooth softening is observed, especially for the higher value of A. The softening behaviour may be attributed to the degradation of anisotropy when the soil reaches the failure conditions at the critical state.Finally a more realistic case, where both soil and stress anisotropies are activated, is presented hereafter. The computed soil reaction curve is compared to the cases where initial and induced anisotropies are considered separately, see Figure5.16. Again, the soil reaction curves tend toward the same plateau. It is concluded therefore that the bearing capacity factor is not affected by the stress or the soil anisotropy. However, the lateral ultimate capacity p u depends on the anisotropy, since the undrained shear strength is given in function of the anisotropy parameters.

Figure 5 .

 5 Figure 5.16: Effect of stress and soil anisotropies on the p-y curve

Figure 5 .

 5 Figure 5.17: Effect of (a) cohesion and (b) structure on the p-y curve

Figure 5 .

 5 Figure 5.18: Schema of the centrifuge test

134CHAPTER 5 Figure 5 .

 55 Figure 5.19: Profiles of (a) preconsolidation pressure and (b) OCR in the soil before pile installaton, after[START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF]

Figure 5 .

 5 Figure 5.20: Mesh and boundary conditions of the 3D-FE model

Figure 5 .

 5 Figure 5.21: Measured and calculated profiles of undrained shear strength

Figure 5 .

 5 Figure 5.22: Force-displacement curve at pile head, measured curve from[START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF]

Figure 5 .

 5 Figure 5.23: (a) Computed and (b) experimental p-y curves at different depths, from measured curve from[START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF]

Figure 5 .

 5 Figure 5.24: Profile of (a) lateral displacement (b) soil reaction for a pile head displacement of 0.15%B, experiment from[START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF]

Figure 5 .

 5 Figure 5.25: Pile moment profiles for different head forces

Figure 5 .

 5 Figure 5.28: Effective normal p n and shear τ stresses along the soil-pile interface

Figure 5 .

 5 Figure 5.29: Distribution of pore water pressures (PWP) at the soil-pile interface

Figure 5 .

 5 Figure 5.30: Mesh and geometry of the short monopile

Figure 5 .

 5 Figure 5.31 (a) shows the horizontal load versus lateral displacement applied at the pile head. Comparing the results with Figure5.22 for the slender pile, it can be observed that the short pile capacity and displacements are larger. According toByrne et al. (2015b), the ultimate short pile capacity is defined at 10% of the pile diameter instead of 20%B for slender piles (since the displacement is more important for short piles). Therefore, the ultimate capacity in the considered short pile is 1.8M N .When laterally loaded, a displacement and a rotation take place in short piles. The two curves are depicted in Figure5.31. It can be observed that the pile does not rotate until M = 0.5M N/m. Afterwards, the pile rotation increases until 1 • , which exceeds the maximum allowable rotation angle for the Ultimate Limite State (ULS) (equal to θ max = 0.5 • according to[START_REF] Achmus | Behavior of monopile foundations under cyclic lateral load[END_REF]).

Figure 5 .

 5 Figure 5.31: Global response of the short pile in terms of (a) force-displacement curve, and (b) moment-rotation at the pile head

Figure 5 .

 5 Figure 5.32: Short monopile: Profiles of lateral displacement

Figure 5 .

 5 Figure 5.33: Short monopile: bending moment profiles

Figure 6 . 1 :

 61 Figure 6.1: SCLAY-1SB model: intrinsic and bounding surfaces in (a) triaxial plane and (b) multi-stress space

Figure 6 . 2 :

 62 Figure 6.2: Undrained stress path during virgin and unloading process predicted with (a) fixed projection center and (b) moving projection center

Figure 6 . 4 :

 64 Figure 6.4: Computed plastic flow direction in an undrained cyclic loading test with different plastic flow rules: (a) image stress flow rule; (b) hybrid flow rule, from Shi et al. (2018)

Figure 6 . 5 :

 65 Figure 6.5: Bounding and intrinsic yield surfaces of structured clays

Figure 6 Figure 6

 66 Figure 6.6: Isotropic consolidation of reconstituted sample of Georgia clay: experiment from Sheu (1985) and simulation

Figure 6

 6 Figure 6.8: Calibration of the constant H 0

Figure 6

 6 Figure 6.10: Cyclic undrained triaxial test on samples of Georgia Kaolin clay: experiment and simulation

Figure 6

 6 Figure 6.11: Cyclic undrained triaxial test on samples of Georgia Kaolin clay: experiment and simulation

  Figure6.12: Simplified scheme presenting the apparatus of the self-boring pressuremeter, from[START_REF] Robertson | In situ testing and its application to foundation engineering[END_REF]

Figure 6

 6 Figure 6.13: Finite element model of the SBP tests: geometry and boundary conditions

  Figure 6.14: SBP tests: Loading phases of the cylindrical cavity

  lack of undrained cyclic triaxial tests on Merville clay, determination of the constant r controlling the flow rule is unfortunately infeasible. Accordingly, an associate flow-rule is considered in this analysis, i.e r = 0.

Figure 6 Figure 6

 66 Figure 6.15: Undrained triaxial tests on reconstituted samples of Merville clay Experiments from Han (2014) and simulations

Figure 6

 6 Figure 6.17: Results of simulation of the test SBD 1 in terms of (a) pressuremeter curve (b) evolution of the injected volume with loading time

Figure 7 . 1 :

 71 Figure 7.1: Response of the pile section in terms of lateral load-displacement curve

Figure 7 . 2 :

 72 Figure 7.2: Influence of the strain accumulation on the lateral pile displacement

  Figure 7.4 presents the results where the maximum lateral displacement in each cycle is plotted as function of the number of cycles. The curves are normalized with the maximum displacement at the first loading cycle.

Figure 7 . 4 :

 74 Figure 7.4: Comparaison of the lateral displacement of the pile section and the pressuremeter cavity

Figure 7 . 5 :

 75 Figure 7.5: Boundary conditions and mesh of the 3D model

  : SCLAY-1SB parameters for the Georgia clayν κ λ i M c M e ω ω d µ ξ v ξ d ξ h 0 a d a v r 0

Figure 7

 7 Figure 7.6: (a) Load-displacement curve at the pile head and (b) Moment at the ground surface-pile rotation, for monotonic loading

Figure 7

 7 Figure 7.7: Load-displacement curves at the pile head for cyclic loading

Figure 7 . 8 :

 78 Figure 7.8: Evolution of pile displacement and rotation with the number of cycles

Figure 7

 7 Figure 7.9: Pile head displacement profiles for different cycles

Figure 7 .

 7 Figure 7.10: Reaction curves along the pile under 25 cycles of lateral loading

Figure 7 .

 7 Figure 7.12: Velocity field in the soil at the 25 th cycle of loading

Titre:

  Loi de comportement pour les argiles naturelles: de la modélisation des essais de laboratoire au comportement des monopieux offshore Mots clés : Loi de comportement, argile, monopieu, charge latérale, chargement cyclique La contribution des parcs éoliens offshore en termes d'énergie renouvelable ne cesse de croître. L'électricité produite notamment par les éoliennes en mer du Nord représente de ce fait une part importante de l'énergie consommée en Europe. Pourtant, les objectifs de la Commission Européenne en termes de transition énergétique à l'horizon de 2020 sont loin d'être réalisés sans l'optimisation des techniques existantes et le développement de technologies de plus en plus innovantes. Les monopieux sont les fondations les plus populaires dans l'industrie des éoliennes offshores. Ce sont des pieux larges supportant des éoliennes de fortes puissances et situées dans des eaux plus profondes. L'objectif général de cette thèse de doctorat est d'étudier les interactions solstructure des monopieux sous sollicitation latérale, fondés dans les sols argileux. Une attention particulière est donnée aux comportement du sol dans le but de modéliser finement le comportement des sols argileux. Title A constitutive Model for natural Clays: From Laboratory Testing to Modelling of Offshore Monopiles Keywords : constitutive lw, clay, monopile, lateral load, cyclic loading

Table 3 .

 3 3: Summary of advantages and limitations of the chosen RH rules Equation N • Critical state Excessive rotation

	3.14
	3.15
	3.16

  • It does not avoid the excessive rotation of the yield surface for high values of the stress ratio. The RHR of Wheeler et al. (2003a) is henceforth not applicable for stress ratio exceeding the critical state (i.e, |η| > |M |).

Table 4 .

 4 1: Constants and calibration procedure of SCLAY-1SB model

	Constant Meaning	Calibration
	ν	Poisson's ratio	K 0 -consolidation
	κ	Slope of the rebound line	1D-consolidation
	λ i	Slope of the swelling line	1D-consolidation on reconstituted sample (RS)
	M c	Slope of the CSL compression Compression shear triaxial test on RS
	M e	Slope of the CSL extension	Extension shear triaxial test on RS
	H 0	shape function's constant	Triaxial test on OC and RS
	A	RHR's constant	Stress-controlled drained triaxial test on RS
	ω d	RHR's constant	Equation 4.36 on RS
	µ	RHR's constant	Shear triaxial test with |η| > |M | on RS
	ξ	Destructuration's constant	1D-consolidation on structured samples (SS)
	ξ d	Destructuration's constant	Triaxial test on SS
	ξ t	Tensile strength's constant	Triaxial test on SS

Table 4 .

 4 2: Internal variables of the SCLAY-1SB model

	Internal variable Meaning
	e	Void ratio
	p mi	preconsolidation pressure of reconstituted sample
	b	State parameter of the mapping rule
	σ c	Projection center
	α	Anisotropy tensor
	χ	State of structure
	p t	Tensile strength

Table 4 . 4 :

 44 Constants of the SCLAY-1SB model for the different clays

		Elasticity Standard (MCC)	RH		Bonding	OC
	Clays	ν	κ	λ i	M c	M e A ω d	µ	ξ	ξ d	ξ t	H 0
	Otaniemi	0.25 0.04 0.26	1	0.49 10 1.5 228 -	-	-	-
	Bothkennar 0.20 0.04 0.25 1.542 1.02 10 0.93 100 10 0.4 0.5	-
	London	0.25 0.008 0.065 1.2	0.87 0	0	-	-	-	-0.02

Table 4 .

 4 6: Initial values of the internal variables used for the simulation of the oedometer test on Bothkennar clay

	Sample	e 0	p mi 0	α 0 χ 0	p t 0
	Laval	2.01 55kPa 0.6 3.5 8kPa
	Sherbrook 1.86 65kPa 0.6 4 8kPa
	Drained and undrained triaxial tests

Table 4 .

 4 7: Details of tests data on Bothkennar clay[START_REF] Smith | The yielding of bothkennar clay[END_REF]

	Test	Drainage	θ( • )
	LCU1 Undrained test in compression	-
	LCU3	Undrained test in extension	-
	LCD30	Drained test	30
	LCD55	Drained test	55
	LCD70	Drained test	70
	LCD110	Drained test	110
	LCD315	Drained test	315

Table 5 . 2 :

 52 Values of the lateral bearing capacity factor for different interface roughness, after[START_REF] Randolph | The limiting pressure on a circular pile loaded laterally in cohesive soil[END_REF]

	Roughness r (-) N p (-)
	0.00	9.886
	0.04	10.869
	0.6	11.295
	0.8	11.673
	1.0	12.017

Table 5 .

 5 3: SCLAY-1SB parameters for the 2D-FEA

	Elasticity Standard (MCC)	RH		Bonding	OC
	ν	κ	λ i M c	M e	A ω	µ	ξ	ξ d	ξ t	H 0
	Otaniemi 0.25 0.04 0.26 1.2	0.88	10 1.5 228 10 0.1 0.1 1.0

Table 5 . 4 :

 54 Computed and theoretical values of N p

	Roughness r Computed N p Theoretical N p Error
	(-)	(-)	(-)	(%)
	0.00	9.886	9.142	8.07
	0.04	10.869	10.531	3.21
	0.6	11.295	11.088	1.87
	0.8	11.673	11.563	0.95
	1.0	12.017	11.940	0.64

Table 5 .

 5 5: Summary of pile geometry and properties in prototype and model scales

	Parameter Prototype	Model
	Material	Steel	Aluminium
	L	18m	360mm
	D	16m	320mm
	B	0,9m	18mm
	e	15mm	1mm
	E p	2.10 5 MPa 0,74.10 5 MPa
	E p .I p	895MN.m 2 143MN.mm 2

Table 5 .

 5 7: Soil parameters used for 3D-FEA

	Soil layer	ν	κ	λ	ϕ c	n 0 α χ	c	H 0
	Layer 1	0.25 0.060 0.143 10 • 0.54 0 0	0	7
	Layer 2	0.25 0.027 0.143 10 • 0.54 0 0	0	7
	Layer 3	0.25 0.027 0.143 10 • 0.54 0 0 5kPa 7

Table 6 .

 6

1: SCLAY-1SB parameters of Georgia clay

Table 7

 7

.3. Figures 6.10

present the measured and computed stress path and stress-

Table 6 . 2 :

 62 Loading conditions of the SBP testsThe soil is modelled with the SCLAY-1SB model. The model parameters related to the monotonic features have been calibrated using five undrained triaxial tests on intact and overconsolidated samples of Merville clay. Results of the calibration are presented in Figures 6.15 in terms of stress paths and stress-strain curves. The model parameters are summarized in Table6.3 and Table 6.4.

	Test	Depth	p 0	p 1	∆p 1	∆p 2	∆p 3
		(m)	(kPa) (kPa) (kPa) (kPa) (kPa)
	SBD 1	12	318	667	135	232	335
	SBD 2	10	283	290	350	350	1052
	SBD 3	8	210	240	200	200	300
	SBD 4	6	133	182	200	200	300
	4.3 Soil parameters
	Table 6.3: SCLAY-1SB parameters of Merville clay

Table 6 .

 6 4: Initial values of the SCLAY-1SB's state parameters for Merville clay

	e	OCR α χ	p t
	1.44	5	0 5 100kPa

Table 7 .

 7 1: SCLAY-1SB parameters for the Merville Intact clay

	ν	κ	λ i	M c M e ω	ω d	µ	ξ v	ξ d	ξ h 0 a d a v r
	0.2 0.04 0.28 1.00 0.90 15 0.15 100 10 0.20 5 60 0 10 0
	Table 7.2: Initial values of the SCLAY-1SB's parameters for the Merville clay
				e	OCR α χ		p t
				1.44	5	0 5 100kPa

Chapter 4

The SCLAY-1SB Model for Monotonic Loading

Monopiles Under Monotonic Lateral Loading

Contents reaction decreases linearly until reaching a plateau at a displacement of y = 18A S .y c . The parameter A s is given in function of the depth and y c is defined in the same manner as in Matlock curve. More details are provided in [START_REF] Reese | Field testing and analysis of laterally loaded piles om stiff clay[END_REF] and [START_REF] Khemakhem | Etude expérimentale de la réponse aux charges latérales monotones et cycliques d'un pieu foré dans l'argile[END_REF].

• JeanJean p-y curve

Based on centrifuge tests detailed in JeanJean [START_REF] Jeanjean | A framework for monotonic py curves in clays[END_REF]), the following p-y curve is proposed : 5.10) where G max is the initial shear modulus and s u the undrained soil shear strength.

• Zhang p-y curve

Based on extensive laboratory studies accompanied with 2D finite element analysis, [START_REF] Zhang | Scaling of lateral pile p-y response in clay from laboratory stress-strain curves[END_REF] proposed a p-y model that allows construction of p-y curves directly from Direct Simple Shear (DSS) tests, by upscaling the soil stress-strain response found in laboratory. The upscaling of the stress-strain curve to the p-y curve is illustrated in Figure 5.8. A point on the stress-strain curve corresponds to a point on the p-y curve with the same level of mobilisation in lateral bearing pressure, i.e. p /pu = τ /su. The upscaling to the corresponding displacement is done from the shear strain using two scaling coefficients ξ 1 and ξ 2 . Zhang employs a constant coefficient ξ 1 and assumes ξ 2 depending on the interface roughness r:

When a triaxial path is considered, the undrained shear strength is equal to the half of the deviatoric shear strength at the critical state, i.e, s u = 1 2 q L . It follows that:

where M c is the value of the slope of the critical state in p -q plane, either under compression or extension, depending on the triaxial path. p L in equation 5.26 is the mean effective stress at the critical state, given by equation 5.25.

• Plane strain conditions

When the plane strain conditions are verified, the undrained shear strength is equal to the value of shear stress τ xy at the critical state. It reads:

M in this case is determined from equation 4.25, by considering a null value of the Lode angle, i.e sin 3θ = 0. p L in equation 5.27 is the mean effective stress at the critical state, given by equation 5.25.

Chapter 6

The SCLAY-1SB Model for Cyclic Loading CHAPTER 6. THE SCLAY-1SB MODEL FOR CYCLIC LOADING

The new model parameters a d and a v take only positive or null values. When a v is selected equal to zero, the state parameter d is increasing with the plastic strain, leading to a decrease of the plastic modulus, which produce the ratcheting effect. When selecting a positive value of a v and a d = 0, the shakedown effect is predicted due to the accumulation of the volumetric strain. When a v and a d have non-zero values, a transition behaviour from ratcheting to shakedown is predicted. Accordingly, the proposed equation can predict ratcheting and shakedown effects as well as the transition behaviour from shakedown to ratcheting.

Flow-rule

In bounding surface plasticity, the image stress is often chosen as the reference point used to define the plastic flow direction, as is the case of the first version of the proposed model. Within this assumption and according to the critical state theory, when the image stress is located inside the contraction zone (dry-side), plastic shearing under undrained condition induces a reduction of the effective mean stress (equivalent to positive pore pressure variation), whereas increasing of mean effective stress is generated when the image stress is located inside the dilation zone. When the critical state is reached, i.e, η = M , the mean effective stress does not evolve any more, since the rate of volumetric plastic strain is null. This behaviour is illustrated in Figure 6.3, which represents the prediction of an undrained cyclic triaxial tests on a reconstituted sample of Georgia clay, compared to the experimental stress path of the same test. The comparison of the experimental and the predicted behaviour reveals the following points: .12 depicts the flow mechanism in the soil surrounding the pile in the vertical symmetric plane. Two distinct soil flow mechanisms can be identified, namely the wedge mechanism near the ground surface, extending to a depth of 2.5B, and a rotational soil mechanism near the base of the foundation, which is approximately equal to 1.5 times the pile diameter. Therefore, it