Keywords: Homogénéisation périodique, non linéaire, calcul à la rupture, domaine de résistance macroscopique, programmation conique Periodic homogenization, nonlinear, limit analysis, macroscopic strength domain, Second-Order Cone Programming. 2.2 . . . . . . . . . . . . . . . . . . . . . . 64 3.8 Quadratic element: 10-node quadratic tetrahedron. . . . . . . . .

Les capacités de résistance de l'argilite Callovo-Oxfordian (COx), qui est une roche hôte potentielle pour le dépôt souterrain profond de déchets radioactifs de haute activité en France, sont étudiées. À une échelle microscopique, des micros pores peuvent être observés dans la matrice. Une première étape d'homogénéisation a été réalisée afin d'évaluer le critère de résistance de la matrice. L'analyse microstructurale de ce matériau à quelques centaines d'échelle, référencée échelle échelle mésoscopique, montre une matrice argileuse et une distribution aléatoire d'inclusions minérales (quartz et calcite).

Dans le but de déterminer le domaine de résistance à l'argilite COx, un premier outil numérique a été développé dans le contexte du comportement élastoplastique de la matrice. Plusieurs modèles morphologiques du volume élémentaire représentatif ont été considérés, et soumis à un chargement incrémental dans des conditions périodiques jusqu'à la charge limite. A la suite de ce calcul élastoplastique, un point de la frontière du domaine de résistance est obtenu. Ce dernier est alors obtenu par des calculs élastoplastiques successifs.

Une alternative aux simulations élastoplastique directes, des approches cinématiques et statiques du calcul à la rupture sont réalisées. Une méthode du type éléments finis basée sur la construction d'un champ de contrainte (dans l'approche statique) et d'un champ de vitesse (dans l'approche cinématique) est développé dans un outil numérique permettant de calculer une limite inférieure et une limite supérieure de domaine de résistance.

I

t is well known that large classes of natural materials e.g. soils, rocks, chalk, hard/ soft tissues, etc. but also artificial ones such as concrete or plaster, bioengineered tissues, etc., are heterogeneous and multi-component. They often have a rich, complex internal structure. They can be, among other things, composite solids, porous matrices, mixtures, multi-component or contain randomly dispersed inclusions. Some geomaterials are composed of a clay matrix and a random distribution of mineral particles (e.g. quartz, calcite). These are referred as general composites which result from the combination of two or more constituents. They can be continuous reinforced composites (fibre) or discontinuous composites (pores, inclusions) depending on how they are arranged.

Despite the increased demand for capacity studies of natural geomaterials and the composition optimization of artificial ones for performance and cost, the determination of the mechanical properties of these heterogeneous materials is always one of the major challenges in fields related to the engineering sciences: it must be ensured that they can serve as geological barriers for the deep storage of wastes, both solid (radioactive waste) and gaseous (storage of CO2) or applied in the field of civil engineering to studying and emerging techniques in the soils to improve (or/ and reinforce) their resistance capacities for other construction purposes where the internal structure, the microstructure, plays a key role in studying and understanding macroscopic behavior of such materials.

In general, while the use of these heterogeneous materials has become widespread and diversified, the theory of their imitation is greatly hampered by their complex heterogeneity. This makes it difficult to determine their behaviors by conventional mechanical techniques. Consequently calculation and simulation methods for the mechanical behaviors of such materials still require numerous developments, both theoretically (techniques of homogenization) and numerically.

Many studies have mainly concentrated on the use of micromechanics to derive suitable models of the ductile failure of porous materials. Indeed, by using realistic and typical objects and conditions, the so-obtained models and theories are studied, and explicit theoretical achievements are extracted either by analytical or numerical tools then a comparison with experimental sampling is executed.

On the one hand, sedimentary clay-rich formations are considered by many countries as potential hosts for radioactive waste disposal facilities: Callovo-Oxfordian (France), Opalinus Clay (Mont Terri, Switzerland), Boom Clay (Bure and Mol, Belgium) ( [START_REF] Robinet | Effects of mineral distribution at mesoscopic scale on solute diffusion in a clay-rich rock: Example of the callovo-oxfordian mudstone (bure, france)[END_REF], [START_REF]Argile -synthesis -evaluation of the feasibility of a geological repository in an argillaceous formation[END_REF]). More concretely, by [START_REF]Argile -synthesis -evaluation of the feasibility of a geological repository in an argillaceous formation[END_REF], the Callovo-Oxfordian argillites is identified as a potential host rock for the repository which is planned to be situated at the surface in the middle of the Callovo-Oxfordian clay. On the other hand, recent studies have been carried out on this subject by the Multi-Echelle team, Navier laboratory, École National des Ponts et Chaussée for years; it has been possible to evaluate the macroscopic criterion of a porous medium with a purely coherent matrix by limit analysis approaches referred to the conceptual model of Gurson ([43]) who firstly presented a strength criterion of ductile porous media based on a kinematic approach of limit analysis on a hollow sphere model in which the solid is rigid-plastic and conforming to a von Mises criterion ( [START_REF] Mises | Mechanik der festen körper im plastisch deformablen zustand[END_REF]), which is exactly suitable for spherical stress states. The scientific context of this thesis is the determination of the macroscopic properties of composite materials by numerical studies, viz nonlinear homogenization calculations integrated into a code of finite element computation in three-dimensional space with porous matrix conforming to an elliptic criterion then especially applied to a Callovo-Oxfordian claystone at so-called meso scale (hundreds of µm).

The introductory chapter begins with a quick overview of porous materials, concretely Callovo-Oxfordian claystone which is a suitable host rock for the deep underground repository of high-level radioactive waste in France.

The second part is then devoted to the presentation of the calculation methods commonly used by engineers and being the subject of current research. We can distinguish two main families, i.e. analytical methods based on homogenization techniques and finite element numerical approaches. The description of these different methods and the difficulty in implementing them are also briefly discussed. Finally, this bibliographic study concludes with an exposition of the objective of this thesis and its general structure is then described.

An argillite porous geomaterial, Callovo-Oxfordian claystone

An argillite porous geomaterial, Callovo-Oxfordian claystone

Callovo-Oxfordian argillite is the main claystone from the eastern region of France and its geological formation has been selected as the host formation for installation of a future underground radioactive waste disposal facility in France. This sedimentary clayey formation, about 150m thick from bottom to top, is an argillaceous rock which lies at ca.-500m depth ( [START_REF]Argile -synthesis -evaluation of the feasibility of a geological repository in an argillaceous formation[END_REF]). A typical geological schema (from [START_REF]Argile -synthesis -evaluation of the feasibility of a geological repository in an argillaceous formation[END_REF] and [START_REF] Enssle | Fullscale 3d modelling of a nuclear waste repository in the callovo-oxfordian clay. part 1: thermo-hydraulic two-phase transport of water and hydrogen[END_REF]) of the eastern region of France is illustrated as figure 1.1.

As the studies of [START_REF] Chiarelli | Modeling of elastoplastic damage behavior of a claystone[END_REF] and [START_REF] Robinet | Effects of mineral distribution at mesoscopic scale on solute diffusion in a clay-rich rock: Example of the callovo-oxfordian mudstone (bure, france)[END_REF] have confirmed, the basic characteristics of [START_REF]Argile -synthesis -evaluation of the feasibility of a geological repository in an argillaceous formation[END_REF] and [START_REF] Enssle | Fullscale 3d modelling of a nuclear waste repository in the callovo-oxfordian clay. part 1: thermo-hydraulic two-phase transport of water and hydrogen[END_REF]). this claystone depending on stratigraphic sequence are 15% -30% tectosilicates (mainly quartz and feldspars), 20% -25% carbonates (mainly calcite and dolomite), 35% -60% clays rate (mainly illite and interstratified illite/smectite, kaolinite, mica, and chlorite) and 0% -3% pyrite.

From the material composition point of view, Callovo-Oxfordian argillite can Chapter 1. General introduction be conceptualized at two scales ( [START_REF] Cariou | An original constitutive law for callovooxfordian argillite, a two-scale double-porosity material[END_REF], [START_REF] Jougnot | Spectral induced polarization of partially saturated clay-rocks: a mechanistic approach[END_REF]): At meso scale (millimeter to micrometer), is generally described by a connected fine-grained shale matrix surrounding coarser non-clay inclusions (tectosilicates, carbonates and heavy minerals with grain size in the range of 1µm to 50µm). Furthermore, as pointed out by [START_REF] Robinet | Upscaling the porosity of the callovo-oxfordian mudstone from the pore scale to the formation scale; insights from the 3h-pmma autoradiography technique and sem bse imaging[END_REF], the spatial distribution of porosity is mainly controlled by the spatial distribution of the clay matrix (the porosity of the clay matrix is 40% -45% on average) and low porosities (0% -4%) of quartz and carbonate mineral grains. At this scale Callovo-Oxfordian argillite is a matrix-reinforced rock. At the microscopic scale, the clayed matrix itself can be represented as a composite of a solid phase of clay matrix (formed mainly by clay mineral particles and calcite microcrystals) and pores (see figure 1.2 for the microstructure of Callovo-Oxfordian argillite described at different scales). sample at the laboratory [START_REF] Jacquier | The influence of mineral variability of callovo-oxfordian clay rocks on radionuclide transfer properties[END_REF], mineral map and considered as composite material of porous clay matrix and inclusions at meso scale and the clay matrix itself presented as solid phase matrix and pores at micro scale (pointed out in [START_REF] Robinet | Upscaling the porosity of the callovo-oxfordian mudstone from the pore scale to the formation scale; insights from the 3h-pmma autoradiography technique and sem bse imaging[END_REF], [START_REF] Robinet | Effects of mineral distribution at mesoscopic scale on solute diffusion in a clay-rich rock: Example of the callovo-oxfordian mudstone (bure, france)[END_REF], [START_REF] Cariou | An original constitutive law for callovooxfordian argillite, a two-scale double-porosity material[END_REF] and compiled from figures which are the property of Andra [START_REF]Argile -synthesis -evaluation of the feasibility of a geological repository in an argillaceous formation[END_REF]).

1.3. Brief introduction of the micromechanical approach of the nonlinear behavior of geomaterials 1.3 Brief introduction of the micromechanical approach of the nonlinear behavior of geomaterials 1.3.1 Analytical estimation using homogenization approach Historically, the determination of mechanical behavior of heterogeneous materials based on the micro mechanical properties at the microscopic scale (matrix and inclusions) is the first solution. The method of homogenization appears as an interesting alternative one for heterogeneous material properties determination that allows to address the problems in a very strict way. This method is developed in order to substitute a heterogeneous medium by an equivalent homogeneous fictitious medium. This homogeneous medium then behaves on average as the heterogeneous medium providing the macro mechanical properties in a macroscopic scale. Several analytical approaches based on homogenization techniques, notably on the Eshelby model, have been developed in [33], [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF], [START_REF] Halpin | Effects of environmental factors on composite materials[END_REF], [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF], [START_REF] Hashin | Analysis of composite materials[END_REF], [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF], [START_REF] Tsai | Theory of Composites Design[END_REF], etc. Nonlinear homogenization approaches to study porous media (as reviews in [START_REF] Bilger | Effect of a nonuniform distribution of voids on the plastic response of voided materials: a computational and statistical analysis[END_REF], [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF]) has also been introduced by [23] based on the variational principles and by [START_REF] Suquet | Overall Properties of Nonlinear Composites[END_REF] using the estimation techniques.

The theory of limit analysis calculation as formulated by [START_REF] Salençon | Calcul à la rupture et analyse limite[END_REF] constitutes a suitable framework for conducting a strength homogenization, since it requires only the knowledge of the local resistance criterion of the solid phase. However, the direct (analytical or numerical) resolution of the problem on the composite structure is difficult due to the strong heterogeneities of the material constituent. This difficulty can be overcome when it can be assumed that, on the one hand, the size of heterogeneities are small compared to the dimensions of the studied structure, and on the other hand, these heterogeneities are regularly distributed.

In the context of periodic composite materials, the effective properties of those can be determined using the periodic homogenization method. The theoretical basis of this method can be found in [START_REF] Papanicolaou | Asymptotic Analysis for Periodic Structures[END_REF], [START_REF] Sanchez-Palencia | Non-homogeneous media and vibration theory[END_REF] and [56]. It is then possible to determine its resistance capacities on a macroscopic scale by the limit analysis or yield design method applied to the periodic media. The basics of this method have been introduced by [START_REF] Suquet | Analyse limite et homogénéisation[END_REF] in a general framework or by [START_REF] De | Approche fondamentale du calcul à la rupture des ouvrages en sols renforcés[END_REF] in the context of reinforced soils mechanics from the solution to a limit analysis boundary value problem relative to the unit periodic cell. Indeed, this method makes it possible to evaluate the macroscopic resistance criterion of the composite, based on the knowledge of material constituent.

Chapter 1. General introduction

In particular, it is an effective method able to obtain the anisotropy of resistance of such materials.

The first study by [START_REF] Dormieux | An extension of gurson model incorporating interface stresses effects[END_REF] is possible to determine a resistance domain coincides with the exact result, numerically determined by [START_REF] Turgeman | Approches numériques des charges limites pour un matériau orthotrope de révolution en déformation plane[END_REF], in the case where the base cell is subjected to isotropic compression, isotropic tensile or pure shear, but remains very coarse with respect to any loading. The object of the thesis is to determine numerical approaches to the macroscopic criterion by implementing the elastoplastic and the limit analysis using homogenization approach in the context of the periodic media with a porous matrix.

Numerical implementation

As studied and pointed out by [START_REF] Bignonnet | Fourier-based strength homogenization of porous media[END_REF], most of the strength homogenization simulations where discretized using the finite element method and solver either by using incremental elastoplastic procedure (e.g. by [START_REF] Hassen | A homogenization approach for assessing the yield strength properties of stone column reinforced soils[END_REF], [91] according to the periodic homogenization method implemented in the context of an elastoplastic behavior (e.g. [START_REF] Suquet | Elements of homogenization theory for inelastic solid mechanics[END_REF], [START_REF]Boundary conditions[END_REF]) or by directly second-order programming solvers (e.g. [START_REF] Gueguin | A homogenization approach for assessing the behavior of soil structures reinforced by columns or trenches[END_REF], [START_REF] Gueguin | Numerical assessment of the macroscopic strength criterion of reinforced soils using semidefinite programming[END_REF]). Exceptions to these in [START_REF] Bilger | Effect of a nonuniform distribution of voids on the plastic response of voided materials: a computational and statistical analysis[END_REF], [START_REF] Bilger | Bounds and estimates for the effective yield surface of porous media with a uniform or a nonuniform distribution of voids[END_REF], [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF], [72] are carried out by FFT-based method based on the efficient numerical method. This is a direct mesh free method initially introduced in linear elasticity by [75] and [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] with it complex microstructures can be simply solved.

Most of these numerical studies are restricted to the exploration of loading modes comprising a hydrostatic part combined either to a pure shear mode or to an axi-symmetric shear mode, thus involving only two or three values of the Lode angle. The lack of studies combining both complex microstructures and complex loading modes is mainly due to the high computational cost of the involved non linear problems. To over come this limitation, [START_REF] Bignonnet | Fourier-based strength homogenization of porous media[END_REF] was successfully proposed a slight adaptation of the mesh free, FFT-based method (by [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF] and [73]) using the Augmented Lagrangian and Uzawa scheme.

We recall here the basic steps of finite element method (FEM) introduced in [START_REF] Felippa | Introduction to finite element methods[END_REF], applied for macroscopic behavior estimations. These steps are shown schematically in figure 1.3.

First step, mathematical model, or idealization, is a very important process to 1.4. Contents transfer from a physical system to a mathematical model of the system, where the term model is understood. Indeed, the representative elementary volume (section 3. meshed by discretizing it into two (or three)-dimensional finite elements depending on the expression of mathematical model chosen. The boundary conditions and loadings are also introduced in this step. These conditions are generally associated with macroscopic loading modes or in agreement with periodicity conditions in the case of periodic composite materials. Finally, by means of an averaging of the microscopic fields, the macroscopic strain and stress tensors are evaluated. Mechanical properties of the representative elementary volume could then be estimated.

In this thesis, both elastoplastic (chapter 2) and direct second-order cone programming (chapter 3) numerical codes are developed to estimate the macroscopic behaviors of the heterogeneous material studied.

Contents

In order to study the nonlinear behavior of the heterogeneous geomaterials by micromechanical approach, analytical estimation and numerical deployment, the thesis is structured into three main parts. The first deals with the estimation of Chapter 1. General introduction the macroscopic resistance of the composite materials as a result of an elastoplastic procedure. The second one involves macroscopic criterion analytical studies of heterogeneous geomaterials based on the theory of limit analysis using the homogenization method. Finally, the third part, based on adapting numerical aspects, completes the study of the macroscopic strength domain of Callovo-Oxfordian claystone consisting of a matrix with elliptic criterion reinforced by rigid inclusions.

In detail, the thesis is organized as follows:

1.4.1 A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation

The developed finite element method tool is presented in chapter 2 in the context of the compressible behavior of the matrix. Several morphological patterns of the representative elementary volume are considered and subjected to incremental loading in periodic conditions until collapse occurs. As a result of such elastoplastic calculation, one point of the boundary of the strength domain is obtained. Other points can then be reached by successive elastoplastic calculations.

Limit analysis of heterogeneous materials with an elliptic

resistance criterion matrix based on the homogenization approach Chapter 3 presents a numerical development of static and kinematic (using continuous and discontinuous velocity fields) of limit analysis formulation for heterogeneous material with an elliptic resistance criterion matrix using the homogenization approach.

After pointing out the limitations of such an analytical approach, we proceed to the second section to develop a numerical method for the determination of the macroscopic resistance criterion, based on both static and kinematic definition. We show how, by restricting to finite-dimensional subspaces of the space of stress and velocity fields of a base cell, the resolution of these problems are reduced to problems of minimization of a convex function of a finite number of scalar variables.

Contents

Numerical estimates of the Callovo-Oxfordian claystone strength domains

The strength capacities of Callovo-Oxfordian claystone which is a potential host rock for the deep underground repository of high level radioactive waste in France is considered as meso mechanical model by the consideration porosity [START_REF] He | Cohesive mohr-coulomb interface effects on the strength criterion of materials with granular-based microstructure[END_REF], is the main investigation of this chapter. The clay matrix of this material is considered as a porous plastic homogenized material which is formed by a solid phase and inclusions.

Micro-pores can be observed in the matrix at the so-called micro-scale. The microstructure analysis of this material at some hundreds of µm shows a clay matrix and a random distribution of mineral inclusions.

A first strength homogenization step has been performed (in [START_REF] He | Strength properties of a drucker-prager porous medium reinforced by rigid particles[END_REF]) in order to evaluate the porous matrix strength criteria. Aiming the determination of the Callovo-Oxfordian claystone macroscopic strength domain, the second strength homogenization step is performed. The stress-based (static approach) and the velocity-based (kinematic approach) finite element method are used to perform the numerical calculation for a lower bound and upper bound of strength domain approximation of Callovo-Oxfordian claystone which concerns porous materials and a clay matrix reinforced by rigid inclusions, respectively.

After showing how the properties of symmetry and periodicity make it possible to simplify the boundary conditions and the loading imposed on the unit cell (the periodicity conditions replaced by conventional boundary conditions), the macroscopic strength domains are then derived on the basis of numerical limit analysis calculations performed on the representative elementary volume for the different morphological patterns (the rigid inclusions are placed into the matrix following a simple periodic Primitive Cubic and Face-Centered Cubic arrangement) in three-dimensional space including comparison with the results of elastoplastic estimations. The significant influence of the inclusion distribution on the macroscopic strength properties and on the local mechanical fields is studied as well.

Chapter 2

A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation

O

ver at least thirty years, many studies have been developed with the object of predicting the overall strength properties of heterogeneous material from the understanding of resistance characteristics of their individual components (matrix and inclusions) related to such key parameters as reinforcement volume fraction, especially for inclusion-reinforced composites.

In order to study the macroscopic strength domain properties of the composite medium consisted of a matrix and inclusions, both the porous one and the matrix reinforced with rigid cores are investigated.

To begin, the review of the homogenization method is briefly recalled (detailed presentation can be founded in [START_REF] Suquet | Analyse limite et homogénéisation[END_REF], [START_REF] De | Approche fondamentale du calcul à la rupture des ouvrages en sols renforcés[END_REF], [START_REF] Salençon | Yield strength of reinforced soils as anisotropic media[END_REF], [START_REF] Salençon | Stability analysis of reinforced soil structures through a homogenization method (in french)[END_REF], [START_REF]Boundary conditions[END_REF] and [START_REF] De | A homogenization approach to the yield strength of composite materials[END_REF]) and then, according to the periodic homogenization method implemented in the context of an elastoplastic behavior (introduced by [START_REF] Suquet | Elements of homogenization theory for inelastic solid mechanics[END_REF] and [START_REF] De | Approche fondamentale du calcul à la rupture des ouvrages en sols renforcés[END_REF]), we are interested in a numerical formulation of the composite model in the framework of the elastoplastic behavior of the different components.

Next, a finite element method (FEM) based program designed to simulate macroscopic strength domain properties of composite structures with an elliptic matrix criterion, as a result of an elastoplastic procedure, is developed under periodic boundary conditions. Lastly, as an application, the strength domain properties of several morphologies made up of both pore and rigid inclusions surrounded by studied matrix is concerned and the effects of the third invariant of its macroscopic stress are quantitatively discussed.

Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation 2.2 A homogenization approach for the macroscopic strength domain as a result of the elastoplastic procedure 2.2.1 Theory review of the strength homogenization

Representative elementary volume (REV )

To examine macroscopic structure behaviors, we only need to study the material properties at the macroscopic scale if it is statically homogeneous. The value of a macroscopic variable at a point is obtained by averaging the microscopic values of this variable over a certain volume of heterogeneous material. One of the main interesting aims in micromechanics is to obtain the material properties, at the macro structure, from the results of the studies in microscopic structures.

In classical micromechanics, there are two hierarchical level elements, macroelement and micro-element, which belong to two-level hierarchical mechanical structure: macroscopic structure and microscopic structure respectively (see figure 2.1). At the macroscopic scale, a continuum is created of many material points related to the microscopic scale where we can use continuum equations and copy the characteristics of the material as a whole. At the microscopic scale, a REV should be compounded of a huge number of micro-elements (inclusions, matrices, pores, etc.) such that it can be a representative element of the local continuum properties.

As mentioned in [31], if d is presented as the characteristic scale length of the local continuum, typically the inclusion size in heterogeneous materials, then the 2.2. A homogenization approach for the macroscopic strength domain as a result of the elastoplastic procedure condition:

d ≪ l (2.1)
permits us to consider that the elementary volume is representative, where l is the characteristic length of the elementary volume and meets the condition:

l ≪ L (2.2)
Combination of two conditions (2.1) and (2.2) above, then the two-condition of the size of the REV must be:

d ≪ l ≪ L (2.3)

REV periodic boundary condition

The REV boundary, denoted ∂Ω, is composed of two parts: a positive part, ∂Ω + , and a negative part, ∂Ω -, with the associated outward normal vector n + = -n -at a corresponding point x + ∈ ∂Ω + and x -∈ ∂Ω -, respectively (see figure 2.2), that must be fulfill the condition:

{ ∂Ω + ∪ ∂Ω -= ∂Ω ∂Ω + ∩ ∂Ω -= ∅ (2.4)
Mathematically, the periodic boundary condition expresses the periodicity of the velocity vector and anti-periodicity of the stress field as follows:

{ σ (x + ) n + = -σ (x -) n - U per (x + ) = U per (x -) , ∀x + ∈ ∂Ω + and matching x -∈ ∂Ω - (2.5)

Averaging operation

The averaging operation of the field a(x) on the REV is defined as follows:

⟨a⟩ = 1 |Ω| ∫ Ω a (x) dΩ (2.6)
Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation This integration is over the volume Ω of the REV.

n + W ¶W + ¶W - s( ) x n + + s( ) x n - - U per + ( ) x x + x - U per - ( ) x n -

Homogenous boundary conditions

The homogenization is an association of a micro-level tensor field with a macro-level tensor through an averaging operation (more details can be founded in [START_REF] Wang | Introduction Micromechanics and Nanomechanics[END_REF]). In order to establish this connection between the two behaviors, microscopic and macroscopic, we consider the problematic of the macroscopic loading definition applied to the REV.

The first problem is so-called "homogeneous stress" condition (see figure 2.3.(a)), imposed at the point x on the border ∂Ω of the REV, and defined by:

T (x) = Σ.n. (x) , ∀x ∈ ∂Ω (2.7)
Where n (x) is the normal vector at the considered point and Σ is the macroscopic strain.

If the tress tensor σ (x) is statically admissible field with the macroscopic strain Σ and, by definition, is subject to the above boundary conditions and verifies the equilibrium condition: divσ = 0. It can be shown that the imposed macroscopic strain Σ is equal to the average of the stress field σ (x) over the domain Ω:

2.2. A homogenization approach for the macroscopic strength domain as a result of the elastoplastic procedure

Σ = ⟨ σ ⟩ = 1 |Ω| ∫ Ω σ (x) dΩ (2.
ξ (x) = ∈.x, ∀x ∈ ∂Ω (2.9)
Where ∈ is interpreted as the macroscopic strain and ξ (x) denotes the local displacement at point x. We note ε (x) = ∇ s ξ (x) is a kinematically admissible strain field with macroscopic strain ∈ prescribed to the border ∂Ω of the REV. The following average relation is verified:

∈ = ⟨ ε ⟩ = 1 |Ω| ∫ Ω ε (x) dΩ (2.10)

Hill's lemma

Homogeneous boundary conditions (in strain or stress) make it possible to establish the Hill's lemma ([50] and [START_REF] Mandel | Plasticité classique et viscoplasticité[END_REF]). The form of an energy coherence rule concerning the quantity of "strain work" as follows:

⟨ σ * : ε ′ ⟩ = ⟨ σ * ⟩ : ⟨ ε ′ ⟩ (2.11)
Where σ * is a static admissible stress field which satisfying the condition divσ * = 0 and ε ′ is a kinematic admissible strain field which satisfying the condition of geometric Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation compatibility. The statement (2.11) occurs if either the stress fields σ * or the strain fields ε ′ satisfies a uniform type boundary condition or if the periodic boundary conditions are satisfied.

In the case of homogeneous stress imposed on the boundary of the REV or of the periodic boundary conditions, this lemma makes it possible to demonstrate that the macroscopic strain ∈ is equal to the average of the field of the microscopic strains ε follow equation (2.8):

∈ = ⟨ ε ⟩ .
Similarly, for the case of homogeneous strain imposed on the boundary of the REV or of the periodic boundary conditions, the macroscopic stress must be defined as the average of the field of the microscopic stresses σ follow equation (2.10):

Σ = ⟨ σ ⟩ .
Finally, in order to determine the effective rigidity of a medium, a macroscopic stress Σ is imposed on a REV and the macroscopic strain ∈ is calculated by taking the mean of the microscopic strain fields in the REV, or a macroscopic strain ∈ and the macroscopic stress Σ is then calculated by taking the mean of the microscopic stress fields in the REV. Let us denote C hom the effective stiffness tensor of the studied material, the effective law of this behavior is expressed in the following form:

Σ = C hom : ∈ (2.12)

Yield strength properties of the matrix

A composite material as shown in figure 2.4 may be modeled as a 'geo-composite' material constituted of a regular array of inclusions embedded into the matrix. The strength criterion of the matrix is isotropic and can be expressed as follows:

f ( σ ) ≤ 0 (2.

2.2.

A homogenization approach for the macroscopic strength domain as a result of the elastoplastic procedure

Macroscopic strength condition

The macroscopic strength criteria F and the corresponding homogenized domain G hom can be derived from solving a yield design problem attached to the REV and called a auxiliary problem (more expression of the auxiliary problem may be found in [START_REF] Gueguin | A homogenization approach for evaluating the longitudinal shear stiffness of reinforced soils: column versus cross trench configuration[END_REF]). More precisely, the macroscopic strength condition is defined as follows:

Σ ∈ G hom ⇔ F (Σ) ≤ 0 ⇔ { ∃σ statically admissible with Σ ∀x ∈ Ω, f ( σ(x) ) ≤ 0 (2.14)
where Ω is the REV sub-domain occupied by the matrix.

A stress field σ is statically admissible with Σ if it complies with the following conditions:

⋄ In absence of body forces, the equilibrium state (of the stress tensor σ(x) in Ω) is governed by the equilibrium equation:

{ divσ(x) = 0 ∀x ∈ Ω (2.15)
⋄ The stress vector remain continuous across any possible discontinuity surfaces of the stress field:

[ σ ] .n = 0 (2.16)
where

[ σ ]
is the stress jump of the stress tensor σ across such a surface following its unit normal n (see figure 2.5).

n( )

x ⋄ σ.n is anti-periodic ∀x on ∂Ω ⋄ The macroscopic stress Σ is equal to volume average of the microscopic stress field σ over the REV :

[ ] = s s s ( ) ( )-( ) x x x 1 2 2 1
Σ = 1 |Ω| ∫ Ω σdΩ = ⟨ σ ⟩ (2.17)
Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation Note that, the third condition, anti-periodic means σ.n takes opposite values at any couple of points located on the opposite sides of the REV and the forth condition which is related to the macroscopic stress Σ are both expressed in the context of the periodic homogenization.

Macroscopic strength domain as a result of an elastoplastic procedure

The determination of the macroscopic strength domain is carried out by means of a numerical elastoplastic procedure performed on the REV, leading to the evaluation of limit loads along prescribed loading paths. More precisely, the REV is subject to a strain controlled loading path. This means that a macroscopic strain of the form:

∈(t) = λ(t)∆ (2.18)
is prescribed to the REV, where λ(t) is a scalar multiplier increased from zero to its maximum value corresponding to the limit load and ∆ is the macroscopic strain tensor.

The solution of the elastoplastic auxiliary problem consists of finding at each time of the loading path a velocity field u(x) kinematically admissible with the loading path ⟨ ε( u) ⟩ = λ∆ along with a statically and plastically admissible stress field σ associated in each point to the velocity field through the elastoplastic constitutive behavior of the matrix.

t ∈ [0, T ] : ∈(t) = λ(t)∆ → { u(t), σ(t) } → Σ(t) = ⟨ σ(t) ⟩ (2.19)
The corresponding loading path in the space of stresses is pictured in figure 2.6 with the limit load Σ * identified as the intersecting point with the macroscopic strength surface. Such limit loads are characterized by the occurrence of an un-contained plastic flow mechanism on the REV, which means that the stress field in equilibrium with Σ * remains constant while the load multiplier λ can be arbitrarily increased. The associated plastic flow rule being assumed for the elastoplastic constituent materials at the microscopic scale, the macroscopic strain rate and then ∆ is outward normal to the macroscopic strength domain at point Σ * .

Elastoplastic formulation for the composite medium composed by a matrix and inclusions

S kl

S ij G hom
Î l D = .

stress path limit load S * F( 0 S)= . .

Elastoplastic formulation for the composite medium composed by a matrix and inclusions

A volume Ω occupied by a composite medium consisting of a matrix obeying to a perfectly plastic elastic behavior and a rigid inclusion is considered. As pointed out in [START_REF] Ghazi | Modélization multiphasique pour le calcul des ouvgages renforcés par inclusions rigides[END_REF], let's recall here the equations governing the evolution of such a system:

⋄ Equilibrium equation div σ (x, t) + ρF (x, t) = 0 (2.20)
Where σ is the stress tensor of the matrix and ρF is the density of external forces.

⋄ Elastoplastic behavior

The strain tensor ε is expressed as a function of the displacement field ξ as follows:

ε = 1 2 ( grad ξ T + grad ξ ) (2.21)
The strain ε decomposes into the sum of an elastic part ε e (reversible) and a plastic part ε p (irreversible):

ε = ε e + ε p (2.22)
Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation

The stress tensor of the matrix is related to the elastic strains as follows:

σ = C : ( ε -ε p ) (2.23)
where C is the elastic tensor of the matrix.

⋄ Plastic criterion

Assuming that the matrix is elastic perfectly plastic, the elastic domain C of this matrix is conventionally defined as the set of stress fields, within which all evolution is reversible (elastic). This domain is characterized by the function f , so-called load function, as follows:

f ( σ )      < 0 if σ is inside C = 0 if σ is on the border of C > 0 if σ is outside C (2.24)
The plastic criterion of the matrix is then written:

f ( σ (t) ) ≤ 0 ⇔ σ (t) ∈ C (t) (2.25)
⋄ Flow rule The rate of plastic strains is given by the flow rule which is written:

εp = η ∂f ∂σ ( σ ) with η { > 0 if f ( σ ) = ḟ ( σ ) = 0 = 0 else (2.26)
Where η is a plastic multiplier and in this case, the flow rule is said to be associated, reflecting the normality of the rate of plastic strain with respect to the boundary of the elastic domain.

⋄ Boundary conditions σ (x, t) .v (x, t) = T pd (v) on ∂Ω T (2.27) ξ (x, t) = ξ pd (x,

Temporal discretization of the loading process

Where T dp (v) and ξ dp are the prescribed effort and displacement on the complementary parts (∂Ω T and ∂Ω ξ ) of ∂Ω.

Temporal discretization of the loading process

Given the incremental nature of the equations describing the elastoplastic behavior of the composite medium, the loading process is conventionally subdivided into loading steps.

Note that {Q} is a vector of the total loading parameters and {∆Q} (i) is the loading increment applied to the composite medium at time t i , i = 1, ..., n,defined by:

{∆Q} (i) = {Q (i + 1)} -{Q (i)} (2.31)
Assuming that the solution of the problem of evolution at time t i , under the loading {Q i }, in terms of the field of displacement ξ (t i ), stress fields σ (t i ) and plastic strain field σ p (t i ), the solution of the problem of evolution at time t i+1 by the data of the increments of the displacement fields and the stress fields ∆ξ (t i ) such that the total fields ∆σ (t i ) satisfied the equations (2.20) to (2.30). These fields are given by:

ξ (t i+1 ) = ξ (t i ) + ∆ξ (t i ) (2.32) σ (t i+1 ) = σ (t i ) + ∆σ (t i ) (2.33) ε p (t i+1 ) = ε (t i+1 ) - ( C ) -1 : σ (t i+1 ) (2.34)
Assuming the increments of plastic strain are defined by:

∆ε p (t i+1 ) = ε p (t i+1 ) -ε p (t i ) (2.35)
In the following, the dependence of time is omitted. The resolution of the problem of incremental elastoplasticity is equivalent to that of the following problem of elasticity:

{ ∆σ ∆ε p } = Elas [ {∆Q} ; { ∆ε p }] (2.36)
That means a problem of elasticity in which the plastic strain field increments is given.

Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation

The increment of the stress field in the matrix, from (2.23), is written:

∆σ = C : ( ∆ε -∆ε p ) (2.37)
And the increments of the plastic strain ∆ε p are given by the plastic flow rule (2.26), written in incremental form as follows:

∆ε p = ∆η ∂f ∂σ ( σ + ∆σ ) with ∆η { > 0 if f ( σ + ∆σ ) = 0 = 0 if not (2.38)
The two follow cases of evolution are considered, elastic or elastoplastic, or in an analogous manner depending on the position of the stress field σ + ∆σ with respect to the convex C:

• Stress field inside the convex is:

f ( σ + ∆σ ) < 0 (2.39)
The evolution undergone by the matrix is elastic. The increment of the plastic strain ∆ε p is zero.

• Stress field on the border of the convex is:

f ( σ + ∆σ ) = 0 (2.40)
The matrix underwent elastoplastic evolution. The stress field is then formulated, by combining (2.37) and (2.38):

   σ + ∆σ = σ + C : ∆ε = σ + C : ∆ε -C : ∆η ∂f ∂σ ( σ + ∆σ ) and ∆η > 0 (2.41)
Furthermore, the plastically admissible stress field defined as follows:

σ p.a. = proj C . { σ + C : ∆ε } (2.42)
Where notation proj C represents the orthogonal projection on the convex C in the sense of the scalar product defined from the quadratic form of the elastic energy of the matrix, namely:

⟨ σ, σ ′ ⟩ = 1 2 σ : C : σ ′ (2.43)

Temporal discretization of the loading process

It follows that σ p.a. is the solution of the minimization problem:

M in σ ′ ∈C { ( σ + C : ∆ε -σ ′ ) : ( C ) -1 : ( σ + C : ∆ε -σ ′ ) } (2.44)
That is equivalent to the minimization problem under constraints:

     M in σ ′ ( σ + C : ∆ε -σ ′
) :

( C

) -1 :

( σ + C : ∆ε -σ ′ ) f ( σ ′ ) ≤ 0 (2.45)
Then, by classically introducing the Lagrange function L:

L ( σ ′ , λ ) = - ( σ + C : ∆ε -σ ′ ) : ( C ) -1 : ( σ + C : ∆ε -σ ′ ) -λf ( σ ′ ) , λ ≥ 0 (2.
46) The solution σ p.a. of the problem (2.45) satisfies the following conditions under the Kuhn and Tucker theorem [START_REF] Kuhn | Nonlinear programming[END_REF]:

   ∂L ∂σ ′ ( σ p.a. , λ * ) = 0 ∂L ∂λ ( σ p.a. , λ * ) = 0 ⇔      - ( C
) _1 : ( σ + C : ∆ε -σ p.a. ) -λ * ∂f ∂σ ′ ( σ p.a. ) = 0 f ( σ p.a. ) = 0 and λ * > 0 (2.47)
where

σ p.a. = σ + C : ∆ε -C : λ ∂f ∂σ ( σ p.a.
) , f

( σ p.a. ) = 0 and λ * > 0 (2.48)
By comparison with (2.41):

{ σ p.a. = σ + ∆σ λ * = ∆η (2.49)
it then emerges that:

σ + ∆σ = proj. C { σ + C : ∆ε } (2.50)
Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation

Implicit plasticity algorithm

The solution of the elastoplastic evolution problem returns, as demonstrated in previous section, to the resolution of a succession of the elastic problems with prescribed inelastic (plastic) strains. The difficulty consists of the evaluations of inelastic strains which are gradually corrected until convergence (refer to [4] and [START_REF] Ghazi | Modélization multiphasique pour le calcul des ouvgages renforcés par inclusions rigides[END_REF]).

Beginning the iterative procedure (k = 0) with zero plastic strains, these strains are then estimated at the iteration (k + 1) from their values in the iteration (k), as follows:

1. Solution calculation by solving an elasticity problem with imposed inelastic strain, This strain is plastic strain estimated at the iteration (k):

{ ∆σ ∆ε } (k) = Elas [ ∆Q; ∆ε p (k) ] (2.51)
2. Projection of the matrix stress field on the elasticity respective convex,

σ p.a. (k) = proj. C { σ (k) + C : ∆ε (k) } (2.52)
3. Plastic strain calculation at the iteration (k + 1)

∆ε p (k + 1) = ∆ε (k) + ( C ) -1 : ( σ (k) -σ p.a. (k) ) = ∆ε p (k) + ( C ) -1 : ( σ (k) + ∆σ (k) -σ p.a. (k) ) (2.53)
4. This procedure should be continued until convergence which corresponds to the fact that the plastically admissible fields σ p.a. (k) and the statically admissible fields { σ + ∆σ (k) } converge towards the solution in constraint:

lim i→+∞ σ p.a. = lim i→+∞ ( σ + ∆σ (k) ) = σ + ∆σ (2.54)
It should be noted that step 3 and step 4 are relative to the plasticity local treatment, whereas the first step corresponds to an overall elastic calculation of the system.

2.6. Finite element method applied to the elastoplastic calculation of composite structures 2.6 Finite element method applied to the elastoplastic calculation of composite structures

Principle of minimum potential energy

By denoting { ξ ′ } , a kinematic displacement field, and introducing a potential energy function of this field, we can define this by:

W ( ξ ) = ∫ Ω 1 2 [ ε : C : ε ] dΩ - ∫ Ω [ ε : C : ε in. ] dΩ - ∫ Ω (ρF ).ξdΩ - ∫ ∂Ω T T pd ξdS (2.55)
Generalizing the definition in the case of a natural initial state to the situation where a field of inelastic strain

{ ε in.
} must be taken into account. The functional potential energy thus defined satisfies the following minimum principle: Let ξ be the displacement solution of the elasticity problem imposed on inelastic strain, then:

{ ( ξ ′ ) ∈ C.A E ( ξ ) ≤ E ( ξ ′ ) (2.56)
where C.A denotes a kinematically admissible displacement field of the problem.

Let δξ is the deviations of the displacement fields ξ:

δξ = ξ ′ -ξ (2.57)
Deviation of the potential energy of the displacement field ( ξ ′ ) solution is given by:

δE = E ( ξ ′ ) -E ( ξ ) = E ( ξ + δξ ) -E ( ξ ) (2.58)
Taking account of (2.55):

δE = ∫ Ω 1 2 [ ( ε + δε ) : C : ( ε + δε ) ] dΩ - ∫ Ω 1 2 [ ε : C : ε ] dΩ - ∫ Ω [ ( ε + δε ) : C : ε in. ] dΩ + ∫ Ω [ ε : C : ε in. ] dΩ - ∫ Ω (ρF ). ( ξ + δξ ) dΩ - ∫ ∂Ω T T pd . ( ξ + δξ ) dS + ∫ Ω (ρF ) .ξdΩ + ∫ ∂Ω T T pd .ξdS (2.59)
Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation and therefore, after simplification we have:

δE = ∫ Ω 1 2 [ 2ε : C : δε + δε : C : δε ] dΩ - ∫ Ω [ δε : C : ε in. ] dΩ - ∫ Ω (ρF ) δξd - ∫ ∂Ω T T pd .δξdS (2.60)
where δε denotes a strain field associated with a displacement field δξ.

Introducing the solution of the problem, the stress field σ is associated with the strain fields ε by the elastoplastic behavior law as follows:

σ = C : ( ε -ε in. ) (2.61)
Expression (2.58) becomes:

δE = ∫ Ω [ σ : δε ] dΩ - ∫ Ω (ρF ) .δξdΩ - ∫ ∂Ω T T p d.δξdS + ∫ Ω 1 2 [ δε : C : δε ] dΩ (2.62)
where three first terms simplify with the use of the virtual work principle, whereas the last term is positive, which archive the demonstration of the minimum principle (2.56).

Variational formulation

Finding the solution ∆ξ (i) of the incremental problem (2.51) returns to solve the following problem by the application of the principle of the minimum potential energy established above:

{ Calculate ( ∆ξ ) ∈ C.A such as: ∀ ( ∆ξ ′ ) ∈ C.A.: E ( ∆ξ ) ≤ E ( ∆ξ ′ ) (2.63)
This is equivalent to the following variation problem:

{ Calculate ( ∆ξ ) ∈ C.A such as: ∀ ( ∆ξ ′ ) ∈ C.A.: a (( ∆ξ ) , ( ∆ξ ′ )) = L ( ∆ξ ′ ) (2.64)
where a (., .) is the bilinear form associated with the quadratic form defined by the functional potential energy which is written:

a ( ∆ξ, ∆ξ ′ ) = ∫ Ω [( ∆ε : C : ∆ε ′ )] dΩ (2.65)
2.6. Finite element method applied to the elastoplastic calculation of composite structures and L (., .) is a linear form given by:

L (∆ξ ′ ) = Φ (∆ξ ′ ) + Φ in. (∆ξ ′ ) (2.66)
Φ (., .) being the linear form expressing the work of the external forces:

Φ ( ∆ξ ′ ) = ∫ Ω (ρ∆F ) .∆ξ ′ dΩ + ∫ ∂Ω T ∆T pd .∆ξ ′ dS (2.67)
and Φ in (., .) denotes the linear form expressing the potential of the imposed inelastic strain:

Φ in. ( ∆ξ ′ ) = ∫ Ω [ ∆ε : C : ε in.
] dΩ (2.68)

Discretized form of the minimum principle

Considering a spatial discretization of the composite structure into elements, the discretized expression of the strain energy, according to [25], is given by:

W ({∆ξ}) = 1 t 2 {∆ξ} . [K] . {∆ξ} (2.69)
where [K] is the global stiffness matrix and {∆ξ} is the node displacement increments vector of the composite structure.

The potential of the external forces is written in a similar way in matrix form:

Φ ({∆ξ}) = t {∆F } . {∆ξ} (2.70)
Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation

While the potential of the inelastic strain given by (2.68), can be expressed in the form:

Φ in. ({∆E}) = ∫ Ω t {∆E in. } . [D] . {∆E} dΩ = Ne ∑ i=1 ∫ (νe)i t {∆E in. }. [D] . {∆E} (x, y) dxdy = Ne ∑ i=1 ∫ νR t {∆E in. }. [D] . {∆E} (ξ, η) .Jdξdη = Ne ∑ i=1 ∫ νR t {∆E in. } i . [D] . {B} (∆ξ) i .Jdξd = Ne ∑ i=1   t {∆E in. } i . ∫ νR [D].[B] e .Jdξdη   t {∆F in. } i .{∆ξ} i = t {∆F in. } i . {∆ξ} (2.71)
We thus obtain the discretized expression of the potential energy:

E ({∆ξ}) = W ({∆ξ}) -Φ ({∆ξ}) = 1 2 t {∆ξ} . [K] . {∆ξ} -t {∆F } . {∆ξ} -t {∆F in. } . {∆ξ} (2.72)
So the solution for the minimum leads to the following linear system:

[K] . {∆ξ} = ({∆F } + {∆F in. }) (2.73)

Numerical calculation of the strength domain based on elastoplastic calculations

This section aims to develop a FEM program based on elastoplastic FEM method, specially designed to handle an elliptic criterion matrix, for both studying the macroscopic strength estimation and studying the effect of the microstructure morphology on the homogenized strength properties.

Yield strength properties of the matrix

The matrix obeys to an isotropic elliptic strength criterion which can be expressed as follows:

2.7. Numerical calculation of the strength domain based on elastoplastic calculations

f ( σ ) = ( σ m -c a ) 2 + ( σ d b ) 2 -1 ≤ 0 (2.74)
where a, b and c depend on inter-granular interfacial strength properties and on porosity of the matrix. σ m and σ d are the mean stress and deviatoric stress which are defined as:

σ m = 1 3 tr ( σ ) ; σ d = √ σ d : σ d ; σ d = σ -σ m 1 (2.75)
In the plane (σ m , σ d ), the boundary of this resistance criterion is an ellipse with axes a 2 along σ m and b 2 along σ d .

REV morphologies

Far from addressing the problem in the general case, we will illustrate through numerical calculations, the resolution of the problem of elastoplastic computation posed on a REV in explicit cases where they are made of a homogeneous matrix and centered pores or rigid cores.

The considered REV is formed of a homogeneous matrix and a network of periodically pores or rigid cores. Such that the REV is a cube of the unit side.

A FEM code, specially designed to handle the elliptic plasticity criterion (2.74) of the matrix, is developed and used to construct point by point the macroscopic strength domain in the three following cases:

• Spherical rigid Cores (SC) are in the axisymmetric condition: The REV is a sphere rigid core surrounded by an elastic perfectly plastic matrix (figure 2.7.(a)).

• Rigid cores or porous inclusions are placed into the matrix following a simple periodic Primitive Cubic (PC) arrangement (figure 2.7.(b)).

• Rigid cores or porous inclusions are placed into the matrix with respect to a Face-Centered Cubic (FCC) arrangement (figure 2.7.(c)).

For the PC and FCC morphologies, two cases are considered: the heterogenous material is either similar to a granular material related to porous one; either similar to a matrix-inclusion composite related to matrix reinforced with rigid cores.

Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation

(b) (c) (a)

Boundary conditions

For the first morphology, SC, axisymmetric uniform strain boundary conditions are applied to the REV, whereas periodic boundary conditions with three-dimensional loading are applied to the REV in the two other cases, PC and FCC.

Loading modes

For the first configuration, the macroscopic strain tensor ∆ is axisymmetric and given by:

∆ (α) = cos α∆ 1 + sin α∆ 2 (2.76)
where

∆ 1 = 1 3 1, ∆ 2 = (e r ⊗ e r + e θ ⊗ e θ -2e z ⊗ e z ) (2.77)
α is varied from 0 to π.

For the second and third configurations, the following expression is considered:

∆ = cosα∆ 1 + sin α ( cos βs 1 + sin βs 2 ) (2.78)
where

s 1 = 1 √ 6 (2e 3 ⊗ e 3 -(e 1 ⊗ e 1 + e 2 ⊗ e 2 )) s 2 = 1 √ 2 (e 1 ⊗ e 1 -e 2 ⊗ e 2 ) (2.79)
The angular stride on the angles α and β involved in the macroscopic strain rate ∆ is set to π 180 . is the inclusion volume fraction.

Uniform strain rate boundary conditions applied on the outer surface of the hollow sphere are:

ξ = λ∆.x (2.80)
where the macroscopic strain rate ∆ follows equation (2.76).

The obtained macroscopic strength domains of the sphere made of a rigid core surrounded by an elastic perfectly plastic matrix are displayed in figures 2.9, 2.10 and 2.11 for three values of the aspect ratio of the matrix elliptic criterion a 2 b 2 =1, 10 and 100, respectively. For each value of the aspect ratio of the matrix elliptic criterion, the different values of the volume fraction of the inclusion ( ra r b = 0.3, 0.5, 0.7 and 0.8) are performed.

For the different simulations, c (in 2.74) is kept equal to zero.

Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation

In general cases, if c is not equal to zero (c ̸ = 0), macroscopic strength domains can be obtained by a translation of those obtained and presented in figures 2.9, 2.10 and 2.11 with a distance of c along the Σ m axis.

It should be emphasized that the macroscopic strength domain properties of the first configuration depends explicitly on the third stress invariant since the macroscopic strength domains are not symmetric with respect to q = 0 axis, with q = √ 2 3 (σ rr -σ zz ).

Second configuration: the PC model

The inclusions are displayed in the matrix as a periodic pattern so according to the symmetry of the problem, only the eighth of the REV is considered (see figure 2.12.(a)).

The boundary condition follows equation (2.80) are applied on the outer surface of the REV as follows: Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation

ξ + -ξ -= λ∆ ( x + -x -) (2.81)
The macroscopic strain rate ∆ follows equation (2.77). Indeed, as pointed out in [START_REF] Samir | Détermination du critère de résistance macroscopique d'un matériau hétérogène à structure périodique : apporche numérique[END_REF], due to the symmetry properties of the considered problem, the periodic conditions imply uniform strain boundary conditions. The adopted mesh of 375 20-node hexahedral three-dimensional FEM-elements and 1991 nodes is represented in figure 2.12. b. More details of the properties of 20-node hexahedral FEM element can be founded in [35].

The results obtained for macroscopic strength criteria of the periodic cubic The macroscopic strain rate ∆ still follows equation (2.77).

(a) (b)
The results obtained for macroscopic strength criteria of the FCC volume The same conclusions can be drawn for the PC model and FCC volume concerning the dependency of the strength domain properties on the third stress invariant. These effects increase with the inclusion volume fraction (ρ).

Conclusions

The determination of the macroscopic strength domain properties of composite materials displaying rigid cores or porous inclusions placed in the matrix is studied from the elastoplastic solution of an elastoplastic problem attached to the representative elementary volume. Chapter 2. A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation

The resulting properties of strength domains represented in the space of stresses show a clear dependency on the stress third invariant, especially for the spherical inclusion placed into the matrix following a periodic cubic arrangement (PC model). The numerical results clearly show that these effects increase with the inclusion volume fraction (ρ) and with the matrix criterion aspect ratio ( a b ). The obtained strength domains are compared to those estimated by means of the FFT-based method. The agreement is excellent in all studied cases, more details can be found in [START_REF] Bignonnet | Fourier-based strength homogenization of porous media[END_REF].

The FEM numerical development tool for the macroscopic strength domain calculation of composite materials in this chapter will be used to calculate and compare to those estimated by means of the static and kinematic approach of limit analysis for Callovo-Oxfordian claystone in the next steps of the study (chapter 3 and chapter 4

L

imit analysis using a strength homogenization approach is an innovative micromechanics technique which was first introduced by [START_REF] Suquet | Analyse limite et homogénéisation[END_REF], [START_REF] Suquet | Elements of homogenization theory for inelastic solid mechanics[END_REF] and was developed in the 1980s in [START_REF] De | Approche fondamentale du calcul à la rupture des ouvrages en sols renforcés[END_REF], [START_REF] Salençon | Yield strength of reinforced soils as anisotropic media[END_REF], [START_REF] Salençon | Stability analysis of reinforced soil structures through a homogenization method (in french)[END_REF], [START_REF] De | A homogenization approach to the yield strength of composite materials[END_REF], [START_REF] Salençon | An introduction to the yield design theory and its applications to soil mechanics[END_REF], [START_REF] De Buhan | Yield design of reinforced earth walls by a homogenization method[END_REF]. It is widely applied in the modeling of multi-layer heterogeneous material calculations nowadays, especially for geomaterials and composite materials such as frictional porous in [START_REF] Barthélémy | Détermination du critère de rupture macroscopique d'un milieu poreux par homogénéisation non linéaire[END_REF], [START_REF] Bignonnet | Fourier-based strength homogenization of porous media[END_REF], [START_REF] De | A homogenization approach to the yield strength of composite materials[END_REF], [START_REF] Bignonnet | Strength of a matrix with elliptic criterion reinforced by rigid inclusions with imperfect interfaces[END_REF] or (and) reinforced media in [START_REF] Barthélémy | A micromechanical approach to the strength criterion of drucker-prager materials reinforced by rigid inclusions[END_REF], [START_REF] Maghous | Micromechanical approach to the strength properties of frictional geomaterials[END_REF], [91], clay soils in [START_REF] Li | Limit analysis of 2-d and 3-d structures based on an ellipsoid yield criterion[END_REF], [START_REF] Li | Limit analysis of composite materials based on an ellipsoid yield criterion[END_REF], [START_REF] Li | Limit analysis of composite materials with anisotropic microstructures: A homogenization approach[END_REF], stone column reinforced soils in [START_REF] Hassen | A homogenization approach for assessing the yield strength properties of stone column reinforced soils[END_REF], [START_REF] Gueguin | A homogenization approach for assessing the behavior of soil structures reinforced by columns or trenches[END_REF], [START_REF] Gueguin | A homogenization approach for evaluating the longitudinal shear stiffness of reinforced soils: column versus cross trench configuration[END_REF], [START_REF] Gueguin | Numerical assessment of the macroscopic strength criterion of reinforced soils using semidefinite programming[END_REF], [START_REF] Gueguin | Ultimate bearing capacity of a foundation reinforced by columns or cross trenches under inclined loads: a homogenization approach[END_REF], composite frames in [START_REF] Bleyer | Lower bound static approach for the yield design of thick plates[END_REF], [START_REF] Bleyer | Yield surface approximation for lower and upper bound yield design of 3d composite frame structures[END_REF], [START_REF] Bleyer | Numerical methods for the yield design of civil engineering structures[END_REF], and heterogeneous plates in [START_REF] Bleyer | Yield Design of Periodically Heterogeneous Plates[END_REF], [START_REF] Bleyer | Yield design computations on homogenized periodic plates[END_REF], [START_REF] Bleyer | Breakage of non-Newtonian character in flow through a porous medium: Evidence from numerical simulation[END_REF], [START_REF] Bleyer | A computational homogenization approach for the yield design of periodic thin plates. Part II : Upper bound yield design calculation of the homogenized structure[END_REF], etc.

In this chapter, the developed method by means of two classical limit analysis theorems, both static and kinematic and the homogenization theory, for the macroscopic strength domain determination by nonlinear mathematical programming associated with an elliptic resistance criterion is considered. Indeed, a generalized study of the microstructure as a representative volume element selected from a periodic heterogeneous material to obtain a lower and an upper bound of their strength domain will be studied.

To start with strength properties description, a background of the homogenization method (briefly recalled in section 2.2.1) for the studied case, is introduced with the developed formulation of the support function applied to studied elliptic criterion in the kinematic approach of limit analysis (section 3.2.1).

Next, numerical approaches applied to the strength properties calculation of the materials with the elliptic criterion and periodic boundary conditions are performed in section 3.3.

Finally, limit analysis problems are numerically formulated in section 3.3.3 for the case of the elliptic criterion by solving second-order cone programming (SCOP) problems.

Chapter 3. Limit analysis of heterogeneous materials with an elliptic

resistance criterion matrix based on the homogenization approach 3.2 Limit analysis using the periodic strength homogenization approach

Strength properties description

The local resistance criterion can be expressed by an inequality as follows:

f ( σ ) 0 (3.1)
where σ is the stress tensor.

Local elliptic resistance criterion

Heterogeneous materials (3.1) with an elliptic resistance criterion can be mathematically expressed by an elliptic (second-order polynomial) equation which can be applied for most current yield criteria ( [START_REF] Li | Limit analysis of 2-d and 3-d structures based on an ellipsoid yield criterion[END_REF], [START_REF] Li | Limit analysis of composite materials with anisotropic microstructures: A homogenization approach[END_REF] and [START_REF] Li | Limit analysis of composite materials based on an ellipsoid yield criterion[END_REF]) as follows:

σ T Mσ + σ T N ≤ 1 (3.2)
Where M and N are coefficient matrices related to the strength properties of the material and can be directly defined from its yield criterion. For two-dimensional problems, the stress vector:

σ = (σ 11 , σ 22 , σ 12 ) T (3.3)
or three-dimensional problems, the stress vector:

σ = (σ 11 , σ 22 , σ 33 , σ 23 , σ 31 , σ 12 ) T (3.4)
are presented by the symmetrical stress tensor σ by regrouping its independent components, respectively. If matrix N is equal to vector zeros:

N = ( 0 0 0 0 0 0 ) T (3.5)
then equation (3.2) will be shortened as follows:

f (σ) = σ T Mσ -1 ≤ 0 (3.6)
3.2. Limit analysis using the periodic strength homogenization approach 3.2.1.2 Static description G(x), assumed convex and must include the null stress σ (x) = 0, is an admissible set which can be determined the local strength domain at point x in a continuum Ω, defined by the following general form:

G(x) := { ∀x ∈ Ω, f ( x, σ (x) ) 0 } (3.7)
We get the local strength domain of the elliptic resistance criterion by applying f ( x, σ (x) ) defined by equation (3.6) to equation (3.7):

G(x) := { σ(x) T Mσ (x) -1 ≤ 0; ∀x ∈ Ω } (3.8)
Recall that, the function f is convex which means ∀σ 1 (x) , σ 2 (x) ∈ G(x) and ∀λ ∈ [0, 1] we have:

f ( λσ 1 (x) + (1 -λ) σ 2 (x) ) λf ( σ 1 (x) ) + (1 -λ) f ( σ 2 (x) ) (3.9)
This characteristic of the local strength domain G(x) is illustrated in figure 3.1.

s kl s ij G( ) x Convexity s 1 s 2

Kinematic description

Beside the construction of stress fields σ (x) that satisfy the condition (3.8), an other equivalent description of the strength domain of the material G(x), leads to an resistance criterion matrix based on the homogenization approach external approach based on the principle of virtual rates of work with the construction of kinematically admissible virtual velocity field Û and the calculation of the support π-function. Details of this method are introduced in [27], [START_REF] Salençon | An introduction to the yield design theory and its applications to soil mechanics[END_REF] and [START_REF] Pecker A Salençon | Seismic bearing capacity of circular footings: a yield deisgn approach[END_REF].

The support π-function is the maximum resisting rate of work densities related to G(x) and defined as follows:

π ( d ) = sup σ(x)∈G(x) ( σ : d ) (3.10)
where d is a symmetric second order tensor.

Physically, π-function is a dissipation power and d is a strain rate tensor which is associated with the velocity field Û at all points by the relation:

d (x) = 1 2
( grad Û T + grad Û ) (x) (3.11)
The admissible stress field σ (x) can be determined by equation (3.1) or equivalently characterized by equation as follows:

σ (x) ∈ G(x) ⇔ ∀d , σ : d π ( d ) (3.12)
If the behavior of the velocity field is discontinuity across the surface S (with unit normal vector n). Denoting

[ Û ] (x) = Û 2 (x) - Û 1 (x)
is the velocity jump across the surface S, the expression of the support π-function of the discontinuity velocity field (see figure 3.2) is:

π ( n; [ Û ]) = sup σ(x)∈G(x) { ( σ.n ) . [ Û ]} on S (3.13)
The support π-function of the elliptic resistance criterion is described in the section 3.2.2.5.

Homogenization of strength properties

The replacement of the equivalent homogeneous material with the heterogeneous one in the formulation of yield design problems is an intuitive methodology in the homogenization method, the strength properties of which being specified by means of a macroscopic strength criterion. The macroscopic strength condition and the corresponding domain can be obtained from the solution of a yield design problem attached to the REV (refer to [START_REF] De | A homogenization approach to the yield strength of composite materials[END_REF], [START_REF] Salençon | An introduction to the yield design theory and its applications to soil mechanics[END_REF], [START_REF] Hassen | A homogenization approach for assessing the yield strength properties of stone column reinforced soils[END_REF]).

REV periodic boundary condition

Periodic homogenization implies that periodic boundary conditions are respected on the boundary of the REV : If any periodic material is chosen, this means that the stress field and the velocity field of two points facing each other on two opposite faces must satisfy conditions of anti-periodicity and periodicity, respectively. More details can be found in section 2.2.1.2.

Auxiliary yield design problems and boundary conditions

As pointed out in section 2.2.3, the auxiliary problem and boundary conditions are recalled here. At each point x of Ω described by a local strength domain G(x) which is characterized equivalently by the strength criterion (3.6), the set of statically admissible stress fields σ(x) with a macroscopic stress Σ is denoted S ( Σ ) is defined by the following conditions:

⋄ In absence of body forces, the equilibrium state (of the stress tensor σ (x) in Ω)

is governed by the equilibrium equation:

divσ (x) = 0 ∀x ∈ Ω (3.14)
Chapter 3. Limit analysis of heterogeneous materials with an elliptic resistance criterion matrix based on the homogenization approach ⋄ The stress vector remain continuous across any possible discontinuity surfaces of the stress field:

[ σ ] (x) . n = 0 (3.15)
where

[ σ ] (x)
is the stress jump of the stress tensor σ (x) across such a surface following its unit normal n (see figure 2.5).

⋄ σ.n is anti-periodic ∀x on ∂Ω ⋄ The macroscopic stress Σ is equal to volume average of the microscopic stress field σ over the REV :

Σ = 1 |Ω| ∫ Ω σ dΩ = ⟨ σ ⟩ (3.16)
Similarly, the set of velocity fields Û (x) which are kinematically admissible with a macroscopic strain D which is denoted C ( D ) are defined by the following conditions:

⋄ Û (x)
d ( Û ) = { d ( Û )} + [ Û ] s ⊗ nδ S (3.18)
where

{ d ( Û )} is the regular part of d ( Û )
and δ S is Dirac distribution of the discontinuity surface S.

Limit analysis using the periodic strength homogenization approach

Static definition of homogenized strength

The static definition of the homogenized strength domain G hom of the heterogeneous material with a local resistance criterion (3.1) is a set of macroscopic stress Σ for which it is possible to find a microscopic stress field satisfying the static conditions with Σ at any point x of the REV [START_REF] Salençon | An introduction to the yield design theory and its applications to soil mechanics[END_REF], [START_REF] De | A homogenization approach to the yield strength of composite materials[END_REF]. The macroscopic, so-called F (Σ) can be determined as follows:

Σ ∈ G hom ⇔ F (Σ) ≤ 0 ⇔ ∃σ that { σ(x) ∈ S ( Σ ) f ( σ(x) ) = σ T (x) Mσ (x) -1 ≤ 0, ∀x ∈ Ω (3.19)

Kinematic definition of homogenized strength

The kinematic definition of the homogenized strength domain G hom of a heterogeneous material is the determination of the macroscopic support Π-function (figure 3.3), similarly to (3.12) and (3.13), which is defined by:

Σ ∈ G hom ⇔ ∀D, Σ : D Π hom ( D ) = sup Σ∈G hom Σ : D (3.20)
The same as the static definition, the macroscopic strength domain of the heteroge-

Π hom ( D ) = min Û ∈C(D)    ⟨ π ( d )⟩ = 1 |Ω| ∫ Ω π ( d ) dΩ    (3.21)
If the velocity filed is discontinuous and expressed as equation (3.18), the average of the support π-function could be written:

⟨ π ( d )⟩ = 1 |Ω|   ∫ Ω π ( d ) dΩ + ∫ S π ( n, Û ) dS   (3.22)

Support function of the elliptic resistance criterion

If the local resistance criterion is elliptical, the matrix M in the equation (3.2) is symmetric and positive semidefinite. An alternative statement of this definition can be expressed as follows:

f (x, σ (x)) = σ(x) T Q T Qσ (x) -1 ≤ 0 (3.23) ⇔ G (x) = { σ (x) s.t ∥Qσ (x)∥ 2 1 } (3.24)
where matrix Q is Cholesky factorization of matrix M, σ is a column vector (definition 3.3 and 3.4) of the the stress components of the tensor σ and 'subject to' is abbreviated to 's.t'. By introducing a new variable σ (x) = Qσ (x) (3.25) equation (3.24) can be rewritten:

G (x) = { σ (x) s.t ∥ σ (x)∥ 2 1 } (3.26)
• In case of using continuous velocity fields As concerns the strain tensor, its components are gathered in column vector forms as follows: Following the definition of support π-function (equation 3.10), it can be obtained:

π ( d ) = sup σ∈G ( σ T d ) = sup ∥ σ∥ 1 { σ T Q d } (3.29) where Q = ( Q -1 ) T (3.30)
Denoting

Π d = QZ d (3.31)
with

Z d =          
1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 2 0 0 0 0 0 0 2

          (3.32)
finally the support π-function of the elliptic resistance criterion in case of using continuous velocity fields is:

π ( d ) = Q d = ∥Π d d∥ (3.33)
• In case of using discontinuous velocity fields As pointed out in [START_REF] Makrodimopoulos | Upper bound limit analysis using discontinuous quadratic displacement fields[END_REF], the support π-function of the discontinuity

[ Û ]
of the velocity field across an interior surface S with unit normal vector n can be similarly determined by equation (3.33) by replacing the tensor d with the symmetric tensor χ defined as:

χ = 1 2 ( n ⊗ [ Û ] + [ Û ] ⊗ n ) (3.34)
Truly:

π ( n; [ Û ]) = sup σ(x)∈G(x) { ( σ.n ) . [ Û ]} = sup σ(x)∈G(x) { σ : χ } = π ( χ ) (3.35) ⇔ π ( n; [ Û ]) = Π d χ (3.36)
where χ is a vector of the symmetric tensor χ components:

χ = ( χ 11 χ 22 χ 33 χ 23 χ 31 χ 12 ) T (3.37)
Denoting by n 1 , n 2 and n 3 , the three components of the unit vector n:

n = ( n 1 n 2 n 3 ) (3.38)
Chapter 3. Limit analysis of heterogeneous materials with an elliptic resistance criterion matrix based on the homogenization approach and defining the matrix Z n as:

Z n = 1 2           2n 1 0 0 0 2n 2 0 0 0 2n 3 0 n 3 n 2 n 3 0 n 1 n 2 n 1 0           (3.39)
then the relation of the vector χ and the discontinuity velocity field

[ Û ]
by means of the matrix Z n writes:

χ = Z n [ Û ] (3.40)
The support π-function of the discontinuity velocity field

[ Û ] can be expressed as follows ⇔ π ( n, [ Û ]) = Π [ Û ] [ Û ] (3.41)
where

Π [ Û ] = Π d Z n (3.42)
Depending on the value of strain rate tensor d and

( n; [ Û ])
, the values of the support π-function of the given domain G(x) are either finite or infinite so the kinematic admissible virtual velocity filed Û must be chosen in order that the value of the π-function remains finite. For example, the support π-function of von Mises criterion will be infinite if the condition tr(d) = 0 is not satisfied.

Static and kinematic approaches to the homogenized strength domain

Generally, except some simply particular cases, it is impossible to determine an exact strength domain of a heterogeneous materials by solving the periodic auxiliary problems attached to the representative element volume but can be approximated by optimizing the static (3.19) and kinematic (3.20) problems. The formulations of two classical method so-called interior static and exterior kinematic methods of the yield design theory ( [START_REF] De | Approche fondamentale du calcul à la rupture des ouvrages en sols renforcés[END_REF], [START_REF] Salençon | An introduction to the yield design theory and its applications to soil mechanics[END_REF]) can be considered as optimization problems in which the 3.2. Limit analysis using the periodic strength homogenization approach results (lower bound and upper bound) will be rigorous bounds of the strength domain.

Static approach

The main purpose of the static approach is determining a lower bound of the homogenized strength domain G hom by finding a relevant statically admissible stress field σ(x) ∈ S ( Σ ) based on the static definition of homogenized strength G hom (3.19) which satisfied the local resistance criterion G(x) that allow the macroscopic stress Σ pertains to the strength domain G hom .

Indeed, the first step to determine the approximated lower bound is optimizing the problem (3.19) of a giving load direction Σi in the macroscopic strain space to find a stress field which is statically admissible with that load. Clearly, the maximization of the macroscopic strain Σ i = λ i Σi coincides with the stress field σ i = λ i σi which satisfying the local resistance criterion G(x) at all the points.

S kl S ij G hom G stat ¶G hom l S + i I l S stat i ~i i S convex hull
Denoting by Σ + i = λ + i Σi the macroscopic stress tensor located at the boundary of the strength domain, λ stat i Σ i is a lower bound of Σ + i and λ stat i is defined by:

λ stat i = sup { λ i |∀x ∈ Ω , λ i σ (x) ∈ G(x) } (3.43)
constitutes a lower bound of the extreme value λ + i .

λ stat i ≤ λ + i (3.44) 57
Chapter 3. Limit analysis of heterogeneous materials with an elliptic resistance criterion matrix based on the homogenization approach

The repetition of above process with different load directions in the macroscopic strain space will provide an interior approximation of the strength boundary ∂G hom . This convex hull is called G stat (see figure 3.4).

G stat ⊆ G hom (3.45)
In reality, the choice of the statically admissible stress field and the optimization of the problem (3.43) can be solved by analytical or numerical analysis.

Kinematic approach

The starting point of the kinematic approach is to find a relevant velocity field Û which is kinematically admissible with a macroscopic strain D based on the direct definition (3.20).

In practice, the upper bound approximation of the strength domain G hom can be found first by optimizing equation (3.20) of a given direction of the macroscopic strain D i :

π cine i ( D i ) = min Û ∈C(D i ) { P rm = |Ω| ⟨ π ( d )⟩ ; ∀x ∈ Ω, π ( d ) = sup σ(x)∈G(x) ( σ : d ) } (3.46) Π hom ( D i ) ≤ π cine i ( D i ) (3.47)
The gap between the π cine i

( D i
) -function and the macroscopic support function

Π hom ( D i
) is majoring and depends on the choice of the velocity field. The repetition of above process with different macroscopic strain rate directions in the macroscopic strain space will provide a convex domain shaped by the intersection of all the related half-spaces and corresponding to an exterior approximation G cine of the macroscopic strength domain G hom (see figure 3.5).

3D numerical formulation of the auxiliary problems

3D numerical formulation of the auxiliary problems

Finite element implementation of the static approach

In the static approach, the stress-based finite element method will be used to perform a numerical calculation of the lower bound of strength domain approximation.

The REV is first discretized into n e finite elements Ω (e) , ( Ω = ∪Ω (e) , e = 1, ..., n e )

and a three-dimensional element chosen for use in the static approach is the four-node tetrahedron, called TET4-element [35], allowing to build a linear stress field.

Although the TET4-element is one of the simplest solid elements whose geometry is a right sides tetrahedron and which needs no complex numerical formulation such as integration to construct its properties, it provides enough performance for the stress analysis, especially in this static approach (viz [START_REF] Makrodimopoulos | Upper bound limit analysis using simplex strain elements and second-order cone programming[END_REF] and [START_REF] Makrodimopoulos | Upper bound limit analysis using discontinuous quadratic displacement fields[END_REF]).

Figure 3.6 shows a typical TET4-element. Clearly, its geometry can be definitively defined by four vertex coordinates with respect to the global coordinate system { x y z } :

Chapter 3. Limit analysis of heterogeneous materials with an elliptic resistance criterion matrix based on the homogenization approach

x i , y i , z i (i = 1, ..., 4) (3.48)
In practice, the position of a tetrahedron vertex may be expressed by its global coordinates system (3.48) or by its simply tetrahedral coordinates [35] as follows

{ζ i , i = 1, ..., 4} (3.49)
Recall an interesting property of the TET4-element that any function linear in

{

x y z } , say F (x, y, z), that takes the value {F i , i = 1, ..., 4} at the vertex of the TET4-element could be interpolated in terms of the tetrahedral coordinates:

F (ζ 1 , ζ 2 , ζ 3 , ζ 4 ) = 4 ∑ i=1 F i ζ i (3.50)
Thus each element has its own set of stress fields { σ } (e) i=1,...,4 at four evaluation nodes corresponding with its four vertices; the stress field of the element can be clearly interpolated in terms of those nodal stress fields. Where { σ } (e) i=1,...,4 is the stress tensor of i th node of the e th element, defined as:

{ σ } (e) i =    σ (e)

3D numerical formulation of the auxiliary problems

As mentioned in the previous section, we adopt a column vector by gathering the six independent components of the stress tensor { σ } (e) i to correspond to a stress vector, denoted {σ} (e) i , as follows: The stress field of the e th element is obtained from the stress vector at each evaluation point and from the shape function matrix:

{σ} (e) i = { σ (e)
{σ} (e) (x) = 4 ∑ i=1 N (e) i (x) {σ} (e) i =[N] (e) (x)          {σ} (e) 1 {σ} (e) 2 {σ} (e) 3 {σ} (e) 4          (3.53)
where [N] (e) (x) is the shape function matrix of the element e which is formed:

[N] (e) (x) = [ N (e) 1 (x) [I 6 ] N (e) 2 (x) [I 6 ] N (e) 3 (x) [I 6 ] N (e) 4 (x) [I 6 ] ] (3.54)
The matrix [I 6 ] denotes an identity matrix of size 6 and note that the linear shape functions of the element N The convex criterion (3.6) is satisfied at any point of a TET4-element unless this condition is satisfied at the four vertices of the element.

Indeed, at node i th of the element e th , the elliptic criterion of the stress field is described by:

f ( {σ} (e) i ) ≤ 0 ⇔ {σ} (e) i 2 ≤ 1 (3.57)
where

{σ} (e) i = [Q] (e) i {σ} (e) i (3.58) Matrix [Q] (e)
i is a Cholesky factorization of the coefficient matrix [M] (e) i related to the strength properties of the material at node i of the element e that can be directly Chapter 3. Limit analysis of heterogeneous materials with an elliptic resistance criterion matrix based on the homogenization approach defined from the yield criterion of the material.

The condition (3.58) can be rewritten in the formulation of the global stress vector {σ} and {σ} for all the vertices i (i = 1, ..., 4) of all elements e (e = 1, ..., n e ) as follows:

[Q] {σ} -{σ} = {0} (3.59)
where

{σ} =        . . . {σ} (e) i . . .        , {σ} =        . . . {σ} (e) i . . .        , i = 1, ..., 4; e = 1, ..., n e (3.60)
and

[Q] =     . . . [Q] (e) i . . .     , i = 1, ..., 4; e = 1, ..., n e (3.61)
A statically admissible stress field is built so that the equilibrium and continuity requirements inside any element and between two neighboring elements are satisfied. A body load field over the REV such as gravity or centrifugal forces is ignored so the divergence of the stress field must be equal to zero, we get: Combining this equation with equation 3.53, we obtain the equation of the equilibrium with no body force in the REV :

divσ = 0 ⇔   
[E] {σ} =     . . . [ [E] (e) 1 [E] (e) 2 [E] (e) 3 [E] (e) 4 ] . . .                          . . .
                     = {0} (3.63)
The matrix [E] (e) i is called elementary equilibrium matrix at node i of the e th element and has the following definition: As pointed out in [35], the partial derivatives of the shape functions are:

∂N (e) i ∂x = a i 6Ω (e) ; ∂N (e) i ∂y = b i 6Ω (e) ; ∂N (e) i ∂z = c i 6Ω (e)
(3.65)

So matrix [E] (e)

i can be express in terms of a i , b i , c i which can be determined in terms of the vertices's coordinates.

[E] (e) i = 1 6Ω (e)    a i 0 0 0 c i b i 0 b i 0 c i 0 a i 0 0 c i b i a i 0    (3.66)
Ω (e) is present the volume of the e th element which can be directly calculated from four vertices coordinates:

Ω (e) = 1 6 det      1 1 1 1 x (e) 1
x (e) 2

x (e) 3 The other condition of equilibrium imposes that the stress vector is continuous through the surfaces of discontinuity which are the triangles separating two adjacent elements (figure 3.7). This implies that, between each element of the mesh, the jump of the stress vector is zero (i). Moreover, since periodic conditions are applied to the REV, the stress vector must be anti-periodic on two opposite faces of the REV (ii). These two conditions (i) and (ii) can be solved in a similar way from the numerical point of view. Indeed, two adjacent elements (p) and (q) have a joint-triangular face, denoted (s) (or s + and s -facing together on two opposite faces of the edge of the REV ). The normal vector n (s) of this face is:

n (s) = t { n (s) 1 n (s) 2 n (s) 3 } (3.68)
By introducing the stress vectors at the vertex of each element related to the jointtriangular face (s):

{σ} (p) i , {σ} (p) j , {σ} (p)
k for element (p) and {σ} (q) i , {σ} (q) j , {σ}

(q)
k for resistance criterion matrix based on the homogenization approach (p) element (q), respectively (figure 3.7), we have the following relation:

{ } s k (q) { } s i (p) { } s j (q) { } s j (p) (q) { } s i (q) { } s k (q) n s (s)
[C] {σ} =          . . .    [n] (s) [0] [0] -[n] (s) [0] [0] [0] [n] (s) [0] [0] -[n] (s) [0] [0] [0] [n] (s) [0] [0] -[n] (s)    . . .                                         . . . {σ} (p) i {σ} (p) j {σ} (p) k {σ} (q) i {σ} (q) j {σ} (q) k . . .                                (3.
69) The matrix [C] has the dimensions of (18n s × 36n s ) with n s is number of joint-triangle faces and depending on the normal vector matrix [n] (s) :

[n] (s) =    n (s) 1 0 0 0 n (s) 3 n (s) 2 0 n (s) 2 0 n (s) 3 0 n (s) 1 0 0 n (s) 3 n (s) 2 n (s) 1 0    (3.70)
In addition to the anti-periodicity conditions of the stress vector, the average operation of the stress field must be equal to the macroscopic stress tensor Σ, which constitutes the second specificity related to the resolution of the auxiliary problem.

Σ = ⟨ σ ⟩ = 1 |Ω| ∫ Ω σdΩ (3.71)

3D numerical formulation of the auxiliary problems

By definition, using the property that the stress field is linear on each element, the average operation of the Σ kl component is calculated as follows [START_REF] Gueguin | A homogenization approach for assessing the behavior of soil structures reinforced by columns or trenches[END_REF]:

Σ kl = 1 |Ω| ne ∑ e=1 ∫ Ωe σ (e) kl dΩ (e) = 1 |Ω| ne ∑ e=1 Ω (e) 4 4 ∑ j=1 σ (e) kl,j (3.72)
As explained in section 3.2.2, the macroscopic stress tensor could be obtained for each tensor Σi = Σ(α,β) through the maximization of the multiplication factor λ i = λ (α,β) such as λ (α,β) Σ(α,β) is obtained on the boundary of the static approach of the strength domain.

Thus from equation (3.72) we have the following relation between the unknowns related to the stress field and the macroscopic stress tensor:

λ (α,β) { Σ} (α,β) = [P] {σ} ; α = 0 ÷ 180 • ; β = 0 ÷ 360 • (3.73) where { Σ} (α,β) = { Σ(α,β) 11 Σ(α,β) 22 Σ(α,β) 33 Σ(α,β) 23 Σ(α,β) 31 Σ(α,β) 12 } (3.74)
Matrix [P] is given by the following determinant in terms of the element volumes:

[P ] = 1 |Ω| [ • • • |Ω (e) | 4 [ [I 6 ] [I 6 ] [I 6 ] [I 6 ] ] • • •
] , e = 1, ..., n e (3.75)

Finally, the static approach reduces to the optimization defined by (3.43). For the finite element method developed here, one look to maximize the loading factor subject to the different conditions related to the respect of the local resistance criterion and the static admissibility conditions. More concretely, that means:

λ + ( Σ(α,β) ) ≥ λ static ( Σ(α,β) ) = max λ (α,β) s.t.          {σ} (e) i 2 ≤ 1       {0} [Q] -[I 24ne ] {0} [E] [0] {0} [C] [0] - Σ(α,β) [P] [0]            λ (α,β) {σ} {σ}      =          {0} {0} {0} {0}          (3.77)

Chapter 3. Limit analysis of heterogeneous materials with an elliptic resistance criterion matrix based on the homogenization approach

Next steps, as specified in section 3.2.2, the repetition of the above calculation with a number of different load directions { Σ} (α,β) by variation of α (0 ÷ 180 • ) and β (0 ÷ 360 • ) in the macroscopic stress space will provide a lower bound of the strength boundary ∂G hom . The convex hull G stat which is included in the domain G hom .

Finite element implementation of the kinematic approach

A velocity-based finite element method will be used to perform the numerical calculation for the approximation of the upper bound of strength domain. Piecewise quadratic velocity fields are considered and quadratic tetrahedron TET-10 element (see figure 3.8) are used. The use of discontinuous velocity fields is also discussed. The TET10-element using is a constant metric element (abbreviation: CM ), it must be satisfied the six midsize nodes collocation, at the midpoints between adjacent corners. By mathematically expressed as follows:

x 5 = (x 1 + x 2 ) /2, x 6 = (x 2 + x 3 ) /2, ...x 10 = (x 3 + x 4 ) /2, y 5 = (x 1 + x 2 ) /2, y 6 = (x 2 + x 3 ) /2, ...y 10 = (x 3 + x 4 ) /2, z 5 = (x 1 + x 2 ) /2, z 6 = (x 2 + x 3 ) /2, ...z 10 = (x 3 + x 4 ) /2, (3.78)
The combination of a CM TET10-element with both linearly-varied jump discontinuous and continuous velocity fields are presented in the following items.

3D numerical formulation of the auxiliary problems

⋄ Linear jump discontinuous velocity fields

The REV will be discretized into n e TET10-elements, the velocity field therefore being quadratic on each element. The discontinuity of the velocity field is permitted between two adjacent elements (see figure 3.9). The number of nodes of the mesh is n p and each node has three degrees of freedom which correspond to the three components of its velocity, denoted (u, v, w).

For a given element e, the unknown velocities are gathered in the following vector:

t {U } (e) = { u (e) 1 v (e) 1 w (e) 1 • • • u (e) 10 v (e) 10
w (e) 10

} (3.79)
The strain tensor, denoted d, is linear and coincides with the choosing of velocity fields. The vector of this tensor components at each vertex i of the e th element can be written as:

t {d} (e) i = { d (e)
} ; i = 1, ..., 4 (3.80)
The relation of the strain rate filed {d} (e) i and velocity field {U } (e) of the (p)

{ } s 4 (p) { } s 2 (p) (q) { } s 4 (q) { } s 3 (q)
{ } U 4 (q) { } U 8 (q) { } U 1 (q) { } U 5 (q) { } U 2 (q) { } U 3 (q) { } U 3 (p) { } s 3 (q) { } s 8 (q) { } s 1 (q) { } s 5 (q) { } s 2 (q) { } s 8 (p) { } s 1 (p) { } s 5 (p) { } U 8 (p) { } U 1 (p) { } U 5 (p) { } U 2 (p)
[B] (e) i =               • • •               ∂N (e
∂N (e) j ∂x 0               • • •               , j = 1, ..., 10 (3.82)
The conventional (non-hierarchical) shape functions N The global vector variables {U } and {d} are defined as Related to the discontinuity velocity condition, as mentioned in [START_REF] Makrodimopoulos | Upper bound limit analysis using discontinuous quadratic displacement fields[END_REF], the expression of the linear-jump velocity field of the face s of two adjacent elements: (p) and (q) at point ′ ij ′ (see figure 3.10) can be written as follows:

{U } =        . . . {U } (e) . . .        , {d} =        . . .
[ Û ] ij = 1 2
( [ Û ] i + [ Û ] j ) (3.87)
Using the definition (3.13) of the jump velocity across the surface s at points ′ i ′ , ′ j ′ and ′ ij ′ , extending to all discontinuity surfaces in the global system,

3D numerical formulation of the auxiliary problems

(p) (q) s i k j ik ij ki { } U i (p) { } U i (q) { } U k (q) { } U j (q) { } U j (p) { } U k (p) [ ] U k [ ] U i [ ] U j [ ] U ki [ ] U ij [ ] U ik
[ Ĉd ] {U } =     . . . [ 2 [I 3 ] [I 3 ] [I 3 ] -2 [I 3 ] -[I 3 ] -[I 3 ] ] . . .                                    . . . {U } (p) ij {U } (p) i {U } (p) j {U } (q) ij {U } (q) i {U } (q) j . . .                                = {0}
(3.88) As previously explained (section 3.2.2.4), a macroscopic strain D (α,β) , completely defined by a set of angle (α, β), is applied to the REV. Then the velocity field must be kinematically admissible with this macroscopic strain tensor. This condition leads to a linking equation of the velocity at two opposite points on two opposite faces of the REV (figure 2

.2).

Let us take two opposite faces of the boundary of the REV : ∂Ω + and ∂Ω -. The components of the velocities of the six nodes which belong to one face of the other ones are grouped in an elementary vector {U } i± , such that Chapter 3. Limit analysis of heterogeneous materials with an elliptic resistance criterion matrix based on the homogenization approach

{U } i± = { u e i± v e i± w e i± } ; i = 1, ..., 6 (3.89)
where the node i + belonging to the face ∂Ω + and the opposite node i -is situated on the face ∂Ω -, respectively. The kinematic admissible condition is then expressed in the equation as follows:

{U } i + -{U } i -= D (α,β) {∆x} i (3.90)
where the vector {∆x} i expresses difference between the coordinates of related nodes:

{∆x} i =      x i + -x i - y i + -y i - z i + -z i -      (3.91)
The global matrix [X], in the global system, will then be constructed to express the kinematically admissible boundary conditions of the REV as follows:

[X] {U } =          . . .    1 0 0 -1 0 0 0 1 0 0 -1 0 0 0 1 0 0 -1    . . .                     . . . {U } i + {U } i - . . .            = { X(α,β) } (3.92) with the definition of the vector { X(α,β) } : { X(α,β) } =                  . . .      D 11 (x i + -x i -) + D 12 (y i + -y i -) + D 13 (z i + -z i -) D 21 (x i + -x i -) + D 22 (y i + -y i -) + D 23 (z i + -z i -) D 31 (x i + -x i -) + D 32 (y i + -y i -) + D 33 (z i + -z i -)      . . .                  (3.93)
On the other hand, the other boundary conditions, consisting of fixing one or more components of velocity at certain points can be taken into account by a linear relationship:

[P] {U } = {0}
(3.94)

3D numerical formulation of the auxiliary problems

Finally, it is necessary to calculate the resisting work associated with the velocity field and (expressed previously in 3.21). The mean of the support function π ( d ) (equation 3.41) can be written as follows:

⟨ π ( d )⟩ = 1 |Ω| ( ∫ Ω π ( d ) dΩ + ∫ S π ( n, Û ) dS ) = 1 |Ω| ( ne ∑ e=1 ∫ Ω (e) π ( d ) dΩ (e) + ns ∑ s=1 ∫ S (s) π ( n, Û ) dS (s) ) (3.95)
where n s is the total number of discontinuity faces S (s) .

We specifically mention the simplest numerical integration rule, the one point rule (refer to [35]), which finds application in the evaluation of the resisting work for the CM quadratic tetrahedron element and its discontinuity face in terms of the variables at their vertices, respectively. The equation (3.95) is then rewritten:

⟨ π ( d )⟩ ≈ 1 |Ω| ( ne ∑ e=1 Ω (e) 4 Π (e) d d (e) + ns ∑ s=1 S (s) 3 Π (s) [ Û ] [ Û ] (s) ) (3.96) where        d (e) = 4 ∑ r=1 [I 4 ] {d} (e) r [ Û ] (s) = 3 ∑ t=1 [I 3 ] {[ Û ]} (s) t (3.97)
with r and t are the number of vertices and faces of the quadratic tetrahedron element, respectively. Expression (3.97) can be taken as the configuration in the global system as follows:

{ d } =     . . . [ [I 4 ] [I 4 ] [I 4 ] [I 4 ] ] . . .     {d} = [R] {d} (3.98)
and

{[ Û ]} =     . . . [ [I 3 ] [I 3 ] [I 3 ] ] . . .     {[ Û ]} = [T ] {[ Û ]} (3.

99) resistance criterion matrix based on the homogenization approach

Thus, by referring to the definition (3.46) of the auxiliary problem, the result of the numerical kinematic approach is reduced to the minimization of the maximum dissipation energy under the constraints related to the kinematic admissibility and other conditions:

π kine(α,β) = min {⟨ π ( d )⟩
π hom(α,β) ( D (α,β) ) ≤ π kine(α,β) ( D (α,β) ) (3.101)
where

π kine(α,β) ( D (α,β) ) = min 1 |Ω| ( ne ∑ e=1 |Ω (e) | 4 Π (e) d d (e) + ns ∑ s=1 |S (s) | 3 Π (s) [ [ Û ] ] [ Û ] (s) ) s.t.            -[I 24ne ] [B] [0] [0] [0] [0] [X] [0] [0] [0] [0] [P ] [0] [0] [0] [0] [ Ĉd ] [0] [0] [0] [R] [0] -[I 24ne ] [0] [0] [0] [0] [0] [T ] -[I 9ns ]                             {d} {U } { d } {[ Û ]} {[ Û ]}                  =                      {0} { X(α,β) } {0} {0} {0} {0}                     
(3.102) For the next steps, similarly in the static approach, the repetition of above calculation with a number of different load directions D (α,β) by variation of α and β (α = 0 ÷ 180 • ; β = 0 ÷ 360 • ) in the macroscopic strain space will provide an upper bound of strength domain G hom . G kine is a convex hull which envelope the domain G hom .

⋄ Continuous velocity fields

In case of discontinuity absence, a continuous condition of the velocity field of all points at the surface of two adjacent elements must be added:

3D numerical formulation of the auxiliary problems

[ Ĉc ] {U } =          . . .    [I 3 ] [0] [0] -[I 3 ] [0] [0] [0] [I 3 ] [0] [0] -[I 3 ] [0] [0] [0] [I 3 ] [0] [0] -[I 3 ]    . . .                                         . . . {U } (p) ij {U } (p) i {U } (p) j {U } (q) ij {U } (q) i {U } (q) j . . .                                = {0} (3.
103) The evaluation of the dissipation energy (3.97) can be used with the elimination of the item related to discontinuity expression:

⟨ π ( d )⟩ ≈ 1 |Ω| ( ne ∑ e=1 Ω (e) 4 Π (e) d d (e) ) (3.104)
The result of the numerical kinematic approach is then reduced to the minimization of the maximum dissipation energy under the constraints related to the kinematic admissibility and other conditions:

π kine(α,β) = min {⟨ π ( d )
π hom(α,β) ( D (α,β) ) ≤ π cine(α,β) ( D α,β ) = min 1 |Ω| ne ∑ e=1 |Ω (e) | 4 Π (e) d d (e) s.t.
      -[I 24ne ] [B] [0] [0] [0] [0] [X] [0] [0] [0] [0] [P] [0] [0] [0] [0] [ Ĉc ] [0] [0] [0]            {d} {U } { d }      =            {0} { X(α,β) } {0} {0}            (3.106)

Formulation as a SOCP problem

A variety of material resistance criteria can be written as a system of the second-order cones as introduced in [START_REF] Kuo | Mittelmann. Interior point methods for second-order cone programming and or applications[END_REF] and [START_REF] Bisbos | Second-order cone and semidefinite representations of material failure criteria[END_REF]. Several efficient applications in limit analysis have been developed and widely applied recent years ( [START_REF] Canh | Upper and lower bound limit analysis of plates using fem and second-order cone programming[END_REF], [START_REF] Bleyer | Numerical methods for the yield design of civil engineering structures[END_REF], [START_REF] Makrodimopoulos | Upper bound limit analysis using simplex strain elements and second-order cone programming[END_REF], etc).

SOCP is a nonlinear convex problem which can be a linear or convex quadratic Chapter 3. Limit analysis of heterogeneous materials with an elliptic resistance criterion matrix based on the homogenization approach program and its algorithm is able to solve problems with large amounts of variables with impressive CPU times and can be guaranteed to identify overall optimal solutions ( [START_REF] Canh | Upper and lower bound limit analysis of plates using fem and second-order cone programming[END_REF]).

In this section, the development of how SOCP can be used to solve both lower-bound and upper-bound convex optimization problems is presented.

Convex quadratic optimization

A second order problem consists of optimizing an objective function subject to conic constraints has the following form:

min h T z s.t. ∥A i z + q i ∥ ≤ u T i z + t i ; i = 1, ..., Z (3.107)
where the optimization variables:

z ∈ R n , h ∈ R n , A i ∈ R mxn , q i ∈ R m and t i ∈ R n
are the problem parameters.

For limit analysis problems, the most common second-order cones are as follows:

• Quadratic cone

C r =    z ∈ R k+1 |z 1 ≥ k+1 ∑ j=2 z j 2 = z 2→k+1    (3.108) • Rotated quadratic cone C s = { z ∈ R k+2 |2z 1 z 2 ≥ k+2 ∑ j=3 z j 2 = z 3→k+2 2 ; z 1 , z 2 ≥ 0 } (3.109) m = 3 if the problem is in 3D
and in case of m = 1, the SOCP becomes a linear programming problem.

Lower bound programming

By introducing an auxiliary variable for each node i of the e th element:

{z} (e) i = { z (e) i,1 = 1 2 z (e) i,2 = 1 } (3.110)
where σ(e) k acts as z (e) j ; j = 3, ..., 9; k = 1, ..., 6, respectively. The constraint condition (3.58) can be transformed into a linear equality constraint coupled with a rotated quadratic cone (3.109):

C s , i (e) = { {z} (e) ∈ R 9 : 2z (e) 1 z (e) 2 9 ∑ j=3 ( z (e) j ) 2 ; z (e) 1 , z (e) 2 0 } (3.111)
The maximizing problem (3.77) becomes a rotated conic optimization where the objective function being a load multiplier λ (It is obviously linear, as required for the use of SOCP).

Maximize λ (α,β)

s.t.                    {z} (e) i ∈ C (e) i ; ∀i ∈ {1, ..., 4} , ∀e ∈ {1, ..., n e }         {0} [Q] -[I 24ne ] [0] {0} [E] [0] [0] {0} [C] [0] [0] - Σ(α,β) [P ] [0] [0] {0} [0] [0] [I 2ne ]                    λ (α,β) {σ} {z} = { {σ} {z} }            =                                {0} {0} {0} {0}            . . . { 1 2 1 } . . .                                           (3.112)
here {z} is global auxiliary variable defined as:

{z} =            . . . { z (e) 1 z (e) 2 } . . .            and {z} (e) i = { {z} (e) i {σ} (e) i }

Upper bound programming ⋄ Discontinuous velocity fields

To begin, note that our problem can be modeled by using the following optimization problem:

Mininize m ∑ i=1 w i ∥s i ∥ + n ∑ i=1 v i ∥t i ∥ (3.113) Chapter 3.

Limit analysis of heterogeneous materials with an elliptic resistance criterion matrix based on the homogenization approach

To convert problem (3.113) to a standard SOCP, we have to add constraints

{ ∥s i ∥ ≤ y i , i = 1, ..., m ∥t i ∥ ≤ z i , i = 1, ..., n (3.114)
Then, the problem (3.113) can be modeled by a standard SOCP as follows:

Mininize m ∑ i=1 w i y i + n ∑ i=1 v i z i s.t.      ∥s i ∥ ≤ y i , i = 1, ..., m ∥t i ∥ ≤ z i , i = 1, ..., n y i , z i ≥ 0 (3.115)
We can now apply (3.115) to formulate (3.77) as a standard SOCP problem by adding new variables

   s e = Π (e) d d (e) , ∀e = 1, ..., n e t s = Π (s) [ [ Û ] ] [ Û ] (s) , ∀s = 1, ..., n s (3.116) then π cine(α,β) = Minimize ( ne ∑ e=1 w e y e + ns ∑ s=1 v s z s ) s.t.                                    ∥s e ∥ ≤ y e , ∀e = 1, ..., n e ∥t s ∥ ≤ z s , ∀s = 1, ..., n s y e , z s ≥ 0            -[I 24ne ] [B] [0] [0] [0] [0] [X] [0] [0] [0] [0] [P ] [0] [0] [0] [0] [ Ĉd ] [0] [0] [0] [R] [0] -[I 24ne ] [0] [0] [0] [0] [0] [T ] -[I 9ns ]                             {d} {U } { d } {[ Û ]} { Û }                  =                      {0} { X(α,β) } {0} {0} {0} {0}                      (3.117) where w e = 1 4
Ω (e) Ω and v s = 1 3

S (s) Ω (3.118)

⋄ Continuous velocity fields

Simply, the problem (3.106) can be modeled by the SOCP as follows:

3.4. Conclusions π cine(α,β) = Minimize ne ∑ e=1 w e y e s.t.                      ∥s e ∥ ≤ y e , ∀e = 1, ..., n e ; y e ≥ 0       -[I 24ne ] [B] [0] [0] [0] [0] [X] [0] [0] [0] [0] [P] [0] [0] [0] [0] [ Ĉc ] [0] [0] [0]            {d} {U } { d }      =            {0} { X(α,β) } {0} {0}            (3.119)

Conclusions

This chapter has shown the development of the formulations related to two limit analysis approaches by applying them to solve the auxiliary problems defined in the context of periodic homogenization.

Thanks to the second order cone programming, it is possible to deal with general three-dimension problems corresponding to the elliptic resistance criteria for the determination of the macroscopic resistance domain of heterogeneous materials. In this chapter a FEM based numerical formulation using the statical and kinematical approaches and taking into account periodic boundary conditions has been developed. This formulation is performed to build a numerical FEM based tool to estimated the macroscopic strength domain of the Callovo-Oxfordian claystone, the obtained results are presented in the next chapter. * * *

Introduction

N

atural heterogeneous materials and many engineering composites are multi-phase and multi-scale material systems. As explained in the introduction of Chapter 1, the need to determine the mechanical properties of these composite materials is rather great. More explanation of applying numerical approaches to actual cases in order to estimate macroscopic resistance domains of their heterogeneous materials is presented here. Indeed, the macroscopic resistance criteria of such composite materials, in general, and Callovo-Oxfordian claystone, in particular, cannot be defined simply by analytical formulate; only the use of numerical calculations makes it possible to evaluate them. They are then obtained in the form of polyhedra, depending on a very large number of hyper-planes, which makes their direct use in any set of calculations very difficult, if not impossible. This chapter, by applying the methods posited in chapter 2 and chapter 3, deals with numerical resolutions for the determination of the macroscopic strength criterion of Callovo-Oxfordian argillite, a claystone which displays a clay matrix reinforced by rigid inclusions which will be studied using the limit analysis attached to the several different morphological representative elementary volumes in three-dimensional space. These estimations can be used to determine lower and upper bounds load magnitudes of various classical geotechnical works (e.g. the bearing capacity of a soil under the action of an inclined loading, analysis of the stability of damaged zones, etc.). The study of the dependency of the macroscopic strength domain on the stress third invariant, especially for the spherical inclusion placed into the matrix following a periodic cubic arrangement will be studied as well.

To start with, the micro-meso transition of the Callovo-Oxfordian claystone is presented in section 4.2. Herein, the elliptic resistance criterion is applied to results of the analysis in [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part 1-yield criteria and flow rules for porous ductile media[END_REF] with a suitable description for the strength of the clay matrix in question. Then the transition from the mesoscopic scale to the macroscopic scale which constitutes the second homogenization step to estimate the macroscopic strength domain i.e. the main study of this chapter is presented.

Then, several morphologies are described in section 4.4 after a description of the macroscopic loading in section 4.3.

To obtain lower and upper bound estimations of the macroscopic resistance criterion, in section 4.5, techniques for the approximation of the homogenized strength domain are presented for the static and kinematic approaches, respectively.

It is currently impossible to obtain an exact evaluation of the criterion of the macroscopic resistance domain (except for special cases such as multilayered materials) for a given geometrical and material configurations. We only know that the exact macroscopic resistance domain G hom is between G stat and G cine , but its exact position remains unknown (see figure 4.1).

---

S kl S ij G hom G stat ¶G hom G kine
G stat ⊆ G hom ⊆ G kine (4.1)
Finally, the macroscopic strength criteria are numerically studied and derived in section 4.7 in cases of different matrix criterion aspect ratios by a the developed elastoplastic and limit analysis based tool is presented in section 4.6.

4.2.

Transition from micro to the meso scale: matrices with elliptic criterion and support function 4.2 Transition from micro to the meso scale: matrices with elliptic criterion and support function

As described in section 1.2, at the microscopic scale, micro pores can be observed in clay matrices. A first strength homogenization step has been performed by [23] using the appropriate nonlinear homogenization technique in order to evaluate the matrix strength criteria. The meso-structural analysis of this material at some hundreds of µm scale, referred here as the mesoscopic scale, shows a clay matrix and a random distribution of mineral inclusions (quartz and calcite) which is illustrated in figure 4.2.

An assumption of a separation between the microscopic and mesoscopic scales allows performancing a two-step homogenization: an appropriate homogenization technique [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part 1-yield criteria and flow rules for porous ductile media[END_REF] is used for the first step which deals with the transition from microscopic scale to the mesoscopic scale. Here, the matrix is described as a heterogeneous material made up of pores and a solid phase of the matrix in contact through cohesive frictional interfaces. The result of this first homogenization step is the derivation of the strength properties of the matrix at the mesoscopic scale where it is considered as a homogenous material.

Green criterion

Assuming the results of the analysis in [23] provide a suitable description for the strength of the clay matrix, as introduced by [START_REF] Barthelemy | Approche micromécanique de la rupture et de la fissuration dans les géomatériaux[END_REF] and [31], the strength criterion is isotropic and can be useful to express as an ellipse, called the Green criterion in [START_REF] Bignonnet | Fourier-based strength homogenization of porous media[END_REF], in general as follows:

f G ( σ ) = ( σ m a ) 2 + ( σ d b ) 2 -1 0 (4.2)
where a, b depends on intergranular interfacial strength properties and on the porosity of the material. σ m and σ d are the mean stress and deviatoric stress which are defined as:

σ m = 1 3 tr ( σ ) ; σ d = √ σ d : σ d ; σ d = σ -σ m 1 (4.3)
In the plane (σ m , σ d ), the boundary of this resistance criterion is a half of an ellipse with axes a 2 along σ m and b 2 along σ d . Equation (4.2) can be reformulated as the form of an elliptic resistance criterion (3.6)

f G (σ) = σ(x) T M G σ (x) -1 ≤ 0 (4.4)
where σ is presented to the symmetric stress tensor σ by regrouping its independent components as defined by equation (3.4) and

M G =           k 1 k 2 k 3 0 0 0 k 2 k 5 k 4 0 0 0 k 3 k 4 k 6 0 0 0 0 0 0 k 7 0 0 0 0 0 0 k 8 0 0 0 0 0 0 k 9           (4.5)
with

     k 1 = k 5 = k 6 = 1 9a 2 + 2 3b 2 k 2 = k 3 = k 4 = 1 9a 2 -1 3b 2 k 7 = k 8 = k 9 = 2 b 2 (4.6)
In the framework of limit analysis theory [START_REF] Salençon | An introduction to the yield design theory and its applications to soil mechanics[END_REF], a dual characterization of the strength criterion f G (σ) ≤ 0 is the support function π:

π ( d ) = sup σ(x)∈G(x) ( σ : d ) (4.7)

Macroscopic loading modes on the REV

of the convex set of admissible stress states with the kinematic admissible strain field of d.

The local strength domain at point x in a continuum Ω is now defined as:

G(x) := { f G (σ) = σ(x) T M G σ (x) -1 ≤ 0; ∀x ∈ Ω } (4.
Q = ( Q -1
) T (4.9)

Problem formulation for other yield criteria

The formulations of other yield criteria, such as the Von Mises strength criterion or Hill's yield criterion which are not used here, are presented in Appendix A.

Macroscopic loading modes on the REV

The FEM simulations are carried on the unit cell, accounting of symmetry properties.

• Static approach

For the static approach, the study is restricted to the following imposed macroscopic stress: where u 1 and u 2 are orthogonal unit tensors which are defined as follows: • Kinematic approach The study of the kinematic problem is restricted to the following imposed macroscopic strain rates:

Σ (α,β) = cos (α) 3 1 + sin (α) [ cos (β) u 1 + sin (α) u 2 ] ( 4
u 1 = √ 2
D (α,β) = cos (α) 3 1 + √ 2 3 sin (α) [ -cos ( β + π 3 ) e 1 ⊗ e 1 -cos ( β -π 3 ) e 2 ⊗ e 2 + cos (β) e 3 ⊗ e 3 ] (4.13)
The direction of loading is modified by varying angles α and β (α = 0 ÷ 180 • ; β = 0 ÷ 360 • ) with successive increments of these angles.

For a unit cell displaying three planes of symmetry with normals e 1 , e 2 and e 3 , ( [START_REF] Samir | Détermination du critère de résistance macroscopique d'un matériau hétérogène à structure périodique : apporche numérique[END_REF]), the symmetry properties imply that for the loading modes given by equation (4.13) the macroscopic stress tensor Σ, located at the point of boundary of the homogenized strength domain with outward normal D, is diagonal in the (e 1 , e 2 , e 3 ) basis.

Description of the studied mesostructure

For this model, REV at the mesoscopic scale is a composite structure with reinforcing inclusions surrounded by the homogenized matrix resulting from the micro-meso transition (figure 4.3).

The rigid inclusions are distributed throughout the porous matrix following a regular pattern. The matrix is assumed homogeneous obeying a three-dimensional resistance criterion, characterized by a convex domain in the three-dimensional space. In our concern, the matrix obeys a Green elliptic criterion (4.2).

The inclusions are supposed rigid resistant and perfectly bounded to the matrix. The reinforcement volume fraction ρ (also called replacement ratio) is used and classically defined as the ratio between the volume occupied by the inclusions and the volume of

Morphologies

The structure discussed is conceptualized as a homogeneous matrix comprising a network of inclusions, distributed at a three-dimensional periodicity. The cell includes an inclusion, centered in its middle. A FEM code, specially designed to handle the Green criterion of three-dimensional problems, is developed and used to construct the macroscopic strength domain in the two following cases in the three-dimensional space:

• The rigid inclusions are placed into the matrix following a simple periodicity Primitive Cubic (PC) arrangement. Figure.

Symmetry properties of the considered three dimensional unit cell

The symmetry property of the strength domain is defined by the symmetry property of the elliptic criterion and those of the unit cell. More details of these properties can be found in [START_REF] De | Approche fondamentale du calcul à la rupture des ouvrages en sols renforcés[END_REF] and [START_REF] Wilkins A Mcelwain | Yield criterion of porous materials subjected to complex stress states[END_REF]. In our cases, these properties can be summarized as follows: • The homogenized strength domain G hom takes over the symmetry property of the local elliptic strength domain

G(x) ∀x ∈ Ω, σ ∈ G(x) ⇒ -σ ∈ G(x) ⇒ ( Σ ∈ G hom ⇒ -Σ ∈ G hom ) (4.14)
in the same way, for the macroscopic loading direction:

Π hom ( D (α, β) ) = Π hom ( D (α + π, β) ) (4.15) Σ ( Σ m , Σ d , β = π 3 ) ∈ ∂G hom ⇒ Σ (-Σ m , Σ d , β = 0) ∈ ∂G hom (4.16)
• In the kinematic approach, when the direction e 1 , e 2 and e 3 are equivalent for the REV of the unit cell (figure 4.4) and the local elliptic strength domain is isotropic, the boundary ∂G hom of the macroscopic strength domain has the periodicity of β * = 2 3 π with respect to the Lode angle β, mathematically expressed by the expression below:

Π hom ( D (α, β) ) = Π hom ( D (α, β + β * ) ) (4.
G app low ⊆ G stat ⊆ G hom ⊆ G kine ⊆ G app up (4.18) S kl S ij G hom G stat ¶G hom

Lower bound approximation

The lower bound approximation to the strength domain may be obtained directly by performing the static approach of yield design (3.19). In practice, with a number of different load directions Σ by variation of α (0 ÷ π) and β (0 ÷ 2π) in the macroscopic stress space will provides a lower bound of G hom , which is a convex hull G app low included in the domain G stat .

Upper bound approximation

The upper bound of homogenized strength domain can be practically estimated as follows (refer to [START_REF] Bignonnet | Fourier-based strength homogenization of porous media[END_REF]):

3D-FEM tool development

The computations of all systems are performed using the conic interior-point optimizer MOSEK package [START_REF] Mosek | The mosek optimization software[END_REF] solvers.

• First configuration-The rigid inclusions are placed into the matrix following a periodic cubic arrangement -PC model.

The inclusions are placed into the matrix as a periodic pattern, only an eighth of the unit cell is considered, shown as figure 4 macroscopic stress and strain rate follow equation (4.10) and (4.13), respectively. Indeed, due to the symmetry properties of the considered unit cell, the periodic boundary conditions imply uniform strain rates boundary conditions [START_REF] Samir | Détermination du critère de résistance macroscopique d'un matériau hétérogène à structure périodique : apporche numérique[END_REF].

With PC morphology, the maximum value of the reinforcement volume fraction can be modeled ρ max ≈ 52% where the radius of the sphere is equal to the size of the unit cell.

• Second configuration -The unit cell is a FCC (face-centered cubic) volume.

The unit cell of the rigid inclusions placed into the matrix following a periodic FCC cubic arrangement is shown in figure 4 ρ max ≈ 74% where the sphere's radius is equal to one half of the diagonal of a face of the unit cell.

Numerical results

The effect of mesh density on the numeric calculations

The object of this chapter is the estimation, through numerical calculations, of the convex macroscopic strength domain of the Callovo-Oxfordian claystone described as a periodic structure. Preliminarily, static limit analysis calculation is carried out for a simple morphology of a PC model with a low inclusion density ρ = 10% and the matrix phase follows the Green criterion (4.4) to approximate a lower bound of the strength domain where α and β are equal to zero in equation (4.10) and elliptical parameters to: a = 3.0; b = 5.0 (MPa). This condition is related to a triaxial tension test and the stress tensor reads in the cartesian basis:

Σ =    1/3 0 0 0 1/3 0 0 0 1/3    (4.20)
In the second step, the same material parameters and the number of elements of the PC models are applied with α = β = 0 for macroscopic strain in equation (4.13)

Numerical results

to find the approximated upper bound of the strength do main. Two cases using continuous and discontinuous displacement fields are studied and the matrix of the values of the strain tensor coinciding with a triaxial tension test is expressed hereafter: for a predetermined problem, the CPU time cost depends greatly on the mesh density, but the value of the solutions do not. Furthermore, using static approach and kinematic approach with continuous velocity fields, the solution can be obtained after only few iterations with a short length of CPU running time. On the other hand, comparing the performance of the elements in the kinematic approach, the solution of using continuous and discontinuous displacement fields are closed but using continuous displacement field gives a significantly better solution with a shorter CPU running time required for the given mesh and material. In particular, the CPU time spend is nearly 3.42 hours of optimization process to estimate an upper bound solution in case of fine mesh using discontinuous velocity fields, forty times more than when using continuous velocity fields.

D =    1/3 0 0 0 1/3 0 0 0 1/3    (

Comparison with the elastoplastic FEM

In this section, in order to compare the performance of static and kinematic approaches with elastoplastic FEM simulation, a morphology of Primitive Cubic model The macroscopic strain (2.76) is classically divided into small loading steps and applied incrementally to the unit cell. A kinematically admissible velocity field and a statical and plastic stress field associated to the velocity field are derived at each loading step. The ultimate load is characterized by the occurrence of a failure mechanism of the unit cell, for which the macroscopic stress tensor remains constant whereas the loading can be arbitrarily increased. The adopted mesh of limit analysis calculations are composed of 581 TET4 elements and 2324 nodes (static approach), 581 TET10 elements and 5500 nodes (kinematic approach), respectively ( At low inclusion density (ρ ≈ 10%), the numerical results of limit analysis approaches give the noteworthy agreement with the elastoplastic ones, however the gap between the strength properties derived by the kinematic and elastoplastic approaches exceeds the one obtained for high inclusion densities (figure 4.11).

Estimations performed on the same model properties with increasing inclusion density (ρ ≈ 40%) are shown higher overestimation of the kinematic approach result. A clear dependency on the stress third invariant of the kinematic approach estimation was also seen whose graph is totally similar with elastoplastic one (figure 4.14).

These numerical estimations also show that the sensitivity to inclusion density is highly significant, especially for the kinematic approach and elastoplastic estimations. More studying of this aspect will be presented in the following section.

Study of unit cell effects for the Green elliptic criterion matrix reinforced with inclusions

In studying of unit cell effects for the Green elliptic criterion matrices reinforced with inclusions, the same types of calculations are carried out for aspect ratios of the resistance domain of the clay matrix with a 2 b 2 =1 , 10 and 100. Two volume fractions of inclusions are studied: ρ=10% and ρ=40% reinforcement, respectively.

For the purpose of legibility, only sections with planes Σ m = const are presented. This study indicates that the macroscopic resistance domain has an increased dependence on the Lode angle when the aspect ratio a 2 b 2 of the domain of resistance of the clay matrix is small, as well as volume fractions ρ of inclusions. The same conclusion could be drawn concerning the dependency of the strength domain on the third stress invariant for the second configuration, FCC volume. More clearly for the kinematic approach, these effects increase with the inclusion volume fraction as well.

Conclusions

In this chapter, the determination of the macroscopic strength criterion of the Callovo-Oxfordian argillite, considered as a heterogeneous material with a matrix reinforced by rigid inclusions with the Green criterion, is studied from the solution of limit analysis problems attached to the representative elementary volume for the Primitive Cubic and Face Centered Cubic morphological patterns in three-dimensional space. The studied models concern heterogeneous materials consisted of rigid inclusions surrounded by a clay matrix.

For that purpose, the contribution has targeted the second step of homogenization (transition meso to macro scale) in a two steps homogenization procedure by developing a limit analysis FEM tool to numerically estimate the macroscopic strength domain whereas the first homogenization step deals with the micro-meso scale transition as the result of [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part 1-yield criteria and flow rules for porous ductile media[END_REF].

The obtained strength criteria (the upper and lower bound) of the Callovo-Oxfordian argillite has been determined by numerical graphs and compared to results estimated by means of elastoplastic FEM technique on the basis of the nonlinear homogenization approach successively applied to similar problems ([8], [START_REF] Hassen | A homogenization approach for assessing the yield strength properties of stone column reinforced soils[END_REF], [START_REF] Gueguin | A homogenization approach for assessing the behavior of soil structures reinforced by columns or trenches[END_REF], [START_REF] Gueguin | A homogenization approach for evaluating the longitudinal shear stiffness of reinforced soils: column versus cross trench configuration[END_REF], [START_REF] Gueguin | Numerical assessment of the macroscopic strength criterion of reinforced soils using semidefinite programming[END_REF], [START_REF] Gueguin | Ultimate bearing capacity of a foundation reinforced by columns or cross trenches under inclined loads: a homogenization approach[END_REF], etc.). Moreover, additional comparison of static solutions to the kinematic limit analysis one is derived in order to prove the accuracy of the strength domain estimations.

The resulting strength domains represented in the macroscopic stress space show a clear dependency on the stress third invariant, especially for the spherical inclusion placed into the matrix following a periodic cubic arrangement. It has been clearly shown from the numerical results that these effects increase with the inclusion volume fraction (ρ) and with the matrix criterion aspect ratio ( a b ).

* * * 104

Chapter 5

Conclusions and perspectives

The homogenization method is a powerful alternative direct calculation for studying the performance and the resistance of heterogeneous structures and materials. These homogeneous media provided with the macro mechanical properties at the macroscopic scale behave on average like the heterogeneous medium. It is also an interesting method for studying behaviors and properties of complex three-dimensional morphologies.

By applying the numerical formulae to the resolution of the auxiliary problem, the strength capacities of Callovo-Oxfordian claystone which is a potential host rock for the deep underground repository of high-level radioactive waste in France are investigated.

First, from the point of view of elastoplastic behavior, aiming the determination of the argillite strength domain, a variational approach based on the principle of a discretized minimum is studied to obtain numerical formulae describing the macroscopic elastoplastic behavior of the 'geo-composite' material considered. A numerical FEM tool is developed, in chapter two, in the context of the elastoplastic behavior of the matrix with the elliptic criterion of the heterogeneous material. In the numerical simulations, three morphological patterns of the representative elementary volume are considered and subjected to an incremental loading in periodic conditions until collapse occurs: Spherical rigid Cores (SC, in the axisymmetric conditions), rigid cores or porous inclusions placed into a matrix following a simple periodic Primitive Cubic (CP) arrangement and rigid cores or porous inclusions placed into a matrix with respect to a Face-Centered Cubic (FCC) arrangement. One point of the boundary of the strength domain is obtained as a result of such elastoplastic calculation. The whole boundary strength then could be reached by successive elastoplastic calculations.

Next, with regard to the strength domain of the 'geo-composite' material by an alternative direct method, the use of limit analysis theory, allowed us to treat this by applying the homogenization method to the unit cell at the microscopic scale, it Chapter 5. Conclusions and perspectives is then possible to obtain a definition of the macroscopic strength domain of these heterogeneous materials. The implementation of the static and kinematic approach of limit analysis computations led us to frame numerically, these strength domains in a rigorous way. Starting from hypotheses on the stress and velocity fields (considering both continuous and discontinuous ones), a development of the numerical formulation of the limit analysis approaches for macroscopic strength domains of a composite with elliptic criterion matrix is presented in chapter three.

Finally, desirous of being able to study the strength domain of Callovo-Oxfordian claystone, we applied the numerical tools which is introduced in the chapter two and three to the material studied. Indeed, the result of the first homogenization step (performed by [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part 1-yield criteria and flow rules for porous ductile media[END_REF]) which deals with the transition from the microscopic scale to the mesoscopic scale of the claystone is the derivation of the strength properties of the matrix at the mesoscopic scale where it is considered as a homogenous material. The contribution is targeted the second step of homogenization (transition meso to macro scale) by developing a limit analysis FEM tool to numerically estimate the macroscopic strength domain. The stress-based (static approach) and the velocity-based (kinematic approach) finite element method are used to derive a lower bound and upper bound of strength domain of Callovo-Oxfordian claystone, respectively.

Using the entire numerical approaches mentioned above, it is found that the sensitivity to inclusion density is highly significant, especially for the kinematic approach and elastoplastic estimations. The dependency on the stress third invariant is pointed out in obtained strength domains, especially for the spherical inclusion placed into the matrix following a periodic cubic arrangement. It has been also shown from the numerical results that these effects clearly depend on the inclusion volume fraction and on the matrix criterion aspect ratio.

Future work could focus on an extension of the present methods to deal with other composite materials with an elliptic local strength criterion (e.g. Hill's yield criterion, Tsai-Hill and Tsai-Wu criterion ( [START_REF] Tsai | Strength theories of filamentary structures fundamental aspects of fibre reinforced plastic composites[END_REF], [START_REF] Tsai | A general theory of strength for anisotropic materials[END_REF])), other types of the representative elementary volume morphologies and similar heterogeneous media with imperfect interfaces between constituents.

Contents 1 introduction 1 1 . 1 1 . 2 1 . 3 16 Chapter 1 introduction Contents 1 . 1 3 1. 2 7 1. 3 . 2 8 1. 4 9 1. 4 . 1 A 10 1. 4 . 2 10 1

 1111213161113273284941104210 General Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 An argillite porous geomaterial, Callovo-Oxfordian claystone . . . . . . 5 Brief introduction of the micromechanical approach of the nonlinear behavior of geomaterials . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.1 Analytical estimation using homogenization approach . . . . . . 7 1.3.2 Numerical implementation . . . . . . . . . . . . . . . . . . . . . 8 1.4 Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4.1 A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation . . . . 10 1.4.2 Limit analysis of heterogeneous materials with an elliptic resistance criterion matrix based on the homogenization approach . 10 1.4.3 Numerical estimates of the Callovo-Oxfordian claystone strength domains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2 A homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation 13 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 A homogenization approach for the macroscopic strength domain as a result of the elastoplastic procedure . . . . . . . . . . . . . . . . . . . . General Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . An argillite porous geomaterial, Callovo-Oxfordian claystone 5 1.3 Brief introduction of the micromechanical approach of the nonlinear behavior of geomaterials . . . . . . . . . . . . . . . . 7 1.3.1 Analytical estimation using homogenization approach . . . . . . Numerical implementation . . . . . . . . . . . . . . . . . . . . . . Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . homogenization approach for assessing the macroscopic strength domain -An elastoplastic numerical estimation . . . . . Limit analysis of heterogeneous materials with an elliptic resistance criterion matrix based on the homogenization approach . .

Figure 1 . 1 :

 11 Figure 1.1: Geological schema of the eastern region of France (compiled from [3] and [53]).

Figure 1 . 2 :

 12 Figure 1.2: Microstructure of Callovo-Oxfordian argillite description on different scales:

Figure 1 . 3 :

 13 Figure 1.3: The basic steps of FEM for macroscopic behavior estimation of heterogeneous materials.

Figure 2 . 1 :

 21 Figure 2.1: Microstructure and an REV of the heterogeneous materials.

Figure 2 . 2 :

 22 Figure 2.2: REV periodic boundary conditions.

Figure 2 . 3 :

 23 Figure 2.3: Homogenous boundary conditions.

Figure 2 . 4 :

 24 Figure 2.4: 'Geo-composite' material model.

Figure 2 . 5 :

 25 Figure 2.5: Stress jump field of a discontinue surface.

Figure 2 . 6 :

 26 Figure 2.6: Elastoplastic stress response to a strain controlled loading and associated limit load.

Figure 2 . 7 :

 27 Figure 2.7: REV morphologies. (a): SC in the axisymmetric condition, (b): PC model and (c): FCC volume.

2. 7 .

 7 Numerical calculation of the strength domain based on elastoplastic calculations 2.7.5 Numerical results 2.7.5.1 First configuration: the SC model

Figure 2 . 8 .

 28 Figure 2.8.(a) shows the SC model as a hollow sphere with the points located at the internal surface are fixed. The problem is modeled in axisymmetric conditions and the adopted mesh of T6 FEM-element (the properties of T6 FEM element can be founded in [45] and [35]) is represented in figure 2.8.(b).

r

  a and r b denote the internal and external radii, respectively and ρ = (

Figure 2 . 8 :

 28 Figure 2.8: Sphere rigid core: (a). REV and (b). Adopted FEM mesh.

8 -Figure 2 . 9 :

 829 Figure 2.9: Macroscopic strength domains of the SC model, a 2 b 2 = 1.

2. 7 .Figure 2 . 10 :Figure 2 . 11 :

 7210211 Figure 2.10: Macroscopic strength domains of the SC model, a 2 b 2 = 10.

Figure 2 . 12 : 2 b 2 = 1 Figure 2 . 13 : 2 b 2 = 1 Figure 2 . 14 :

 212221213221214 Figure 2.12: The rigid inclusions are placed into the matrix following a periodic cubic arrangement: (a). REV and (b). Adopted FEM mesh.

Figure 2 . 15 :

 215 Figure 2.15: FCC volume: (a). REV and (b). Adopted FEM mesh.

Chapter 2 . 2 b 2 = 1 Figure 2 . 16 : 2 b 2 = 1 Figure 2 . 17 :

 2221216221217 Figure 2.16: Homogenized strength properties of the FCC volume with ρ = 10%: octahedral plane cuts for Σm Σ max m = 0, 0.5, 0.7, 0.9

Figure 3 . 1 :

 31 Figure 3.1: Convexity of the local strength domain G(x).

3. 2 .SFigure 3 . 2 :

 232 Figure 3.2: The velocity jump [ Û ] of the discontinuous velocity field Û through a surface S.

  is a piecewise continuously differentiable set. ⋄ Û (x) -D.x is periodic ∀x on ∂Ω ⋄ d ( Û ) satisfied equation (3.11) and ⋄ D is equal to the volume average of strain field d (

Figure 3 . 3 : 53 Chapter 3 .

 33533 Figure 3.3: Kinematic definition of the macroscopic strength domain.

d = ( d 11 d 22 d 33 d 23 d 31 d 12 )d 11 d 22 d 33 2d 23 2d 31 2d 12 )

 1212

Figure 3 . 4 :

 34 Figure 3.4: Static approach of the macroscopic strength domain.

Figure 3 . 5 :

 35 Figure 3.5: Kinematic approach of the macroscopic strength domain.

Figure 3 . 6 :

 36 Figure 3.6: Typical four-node tetrahedron, TET4-element using in the numerical static approach.

i 4 ({ 1

 41 (x) are simply the tetrahedra's coordinates: N (e) i (x) = ζ i , i = 1, ..., at node i; i = 1, ..., 4 0 at other nodes. (3.56)

3. 3 .

 3 3D numerical formulation of the auxiliary problems

Figure 3 . 7 :

 37 Figure 3.7: Triangles surface s of two adjacent elements (p) and (q).

λ

  static = max {λ under conditions: (3.57), (3.59), (3.63), (3.69) and (3.73) } (3.76)

Figure 3 . 8 :

 38 Figure 3.8: Quadratic element: 10-node quadratic tetrahedron.

Figure 3 . 9 :Chapter 3 .

 393 Figure 3.9: Internal evaluation points of the strain field in the kinematic approach.

  are given by[35] in terms of the tetrahedral coordinates:N (e) j = ζ j (2ζ j -1) , j = 1, ..., 4 N (e) j = 4ζ k ζ l , j =5, ..., 10 and kl = 12, 23, 31, 14, 24, 34 respectively. (3.83) In the global system, the equation (3.81) can be written in terms of global vector {d} and {U } as follows: {d} -[B] {U } = {0} (3.84)

  ; i = 1, ..., 4; e = 1, ..., n e (3= 1, ..., 4; e = 1, ..., n e (3.86)

Figure 3 . 10 :

 310 Figure 3.10: Velocity jump at the discontinuity face of two adjacent TET10-elements.

  ⟩ under conditions (3.57), (3.88), (3.93), (3.94), (3.98) and (3.103) } (3.105) And in the explicit form:

Figure 4 . 1 :

 41 Figure 4.1: Relative position of the exact macroscopic resistance domain G hom between G stat and G kine

Figure 4 . 2 :

 42 Figure 4.2: Different scales of the argillite Callovo-Oxfordian. (a) macroscopic scale, (b) argillite at mesoscopic scale and (c) porous matrix at microscopic scale.

8 )

 8 Recall section (3.2.2.5) applied for the domain G(x), the support function π ( d ) of the domain G(x) for the Green criterion (4.4) is: • For continuous velocity fields: The support function of the domain (4.8) is follow equation (3.33) and • For discontinuous velocity fields: Clearly, the matrix M G is symmetry and positive semi definitive; it is possible to obtain its Cholesky factorization: Q. The support function of the domain (4.8) in case of using discontinuous velocity fields is defined by equation (3.41) with the matrix Π [ Û ] is defined by equation (3.42), and Q is defined as:

. 10 ) 4 .

 104 Chapter Numerical estimates of the Callovo-Oxfordian claystone strength domains

4. 4 .Figure 4 . 3 :

 443 Figure 4.3: Rigid inclusion reinforced porous matrix and microstructure of the unit cell Ω.

•

  4.4.(a). The unit cell is a Face-Centered Cubic (FCC) volume. Figure. 4.4.(b).

Figure 4 . 4 :

 44 Figure 4.4: Considered microstructure of the unit cell. (a): PC model and (b): FCC model.

17

 17

  )

4. 5

 5 Approximation of the homogenized strength domainBy numerically performing the static and kinematic approaches of the limit analysis, it was possible to determine a lower bound G stat and a upper bound G kine of the macroscopic domain. The definition of these (numerical) resistance criteria depends 4.5. Approximation of the homogenized strength domain on the discretization retained for the loading space.The approximation of the boundary ∂G hom of the homogenized strength domain G hom is gained by the hyper-planes Σ :D = Π app up ( D )in the kinematic approach and point by point in the static approach (figure 4.5).

Figure 4 . 5 :

 45 Figure 4.5: Representation of the different evaluations of a macroscopic strength domain G hom

Figure 4 . 6 :

 46 Figure 4.6: Periodic cubic arrangement model (PC). (a): REV and (b): Adopted mesh of the one eighth model.

Figure 4 . 7 :

 47 Figure 4.7: Face-Centered Cubic volume. (a): REV and (b): Adopted mesh of the one eighth model.

Chapter 4 .

 4 Numerical estimates of the Callovo-Oxfordian claystone strength domains with a inclusion density of ρ ≈ 10% and ρ ≈ 40% is chosen. The simulations are carried out on the one eighth of the PC unit cell, accounting for symmetry properties.The matrix phase follows the Green criterion (4.4) with parameters: a = 7.31; b = 5.78 (MPa) which are roughly corresponding to Callovo-Oxfordian argillite (refer to[START_REF] François | Caractérisation expérimentale et modélisation micromécanique de la perméabilité et la résistance de roches argileuses[END_REF]). The adopted constructed mesh of the elastoplastic FEM is composed of 375 20-node hexahedral elements and 1991 nodes (figure 4.8.(a)).

  figure 4.8.(b)). The angular stride on the angles α and β involved in the macroscopic loads (equation

Figure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Adopted mesh of one eighth PC unit cell. (a): for elastoplastic FEM performance and (b): for limit analysis FEM performance.

Figure 4 . 10 :

 410 Figure 4.10: Homogenized strength properties of the PC model, low inclusion density ρ ≈ 10%, in planes of Σ m = const with different values of Σ m = 0, 0.5, 0.7 and 0.9 of the maximum macroscopic isotropic tensile strength Σ max m . (a): lower bound, (b): elastoplastic FEM and (c): upper bound estimation.

Figure 4 . 11 :Figure 4 . 12 :

 411412 Figure 4.11: The comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (a) and the plane of the equation Σ 11 + Σ 22 = 0 (b) of the corresponding method for the PC model, low inclusion density ρ ≈ 10%.

Figure 4 . 13 :

 413 Figure 4.13: Homogenized strength properties of the PC model, high inclusion density ρ ≈ 40%, in planes of Σ m = const with different values of Σ m = 0, 0.5, 0.7, and 0.9 of the maximum macroscopic isotropic tensile strength Σ max m . (a): lower bound, (b): elastoplastic FEM and (c): upper bound estimation.

Figure 4 . 14 :

 414 Figure 4.14: The comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (a) and the plane of the equation Σ 11 + Σ 22 = 0 (b) of the corresponding method of the PC model, high inclusion density ρ ≈ 40%.

4. 7 . 3 . 1 2 b 2 . 4

 731224 First configuration-PC model • With ρ = 10%, the results obtained for lower bound and upper bound macroscopic strength criteria of the periodic cubic arrangement cell are displayed in the figure 4.15 for different values of the aspect ratio of the matrix elliptic criterion a More results of the static and kinematic approaches are shown in figure B.1, figure B.2 and figure. B.3 (Appendix B) for a 2 b 2 =1, 10 and 100, respectively. • With ρ = 40%, the obtained results for lower bound and upper bound macroscopic strength criteria of the periodic cubic arrangement cell are displayed in the figure. 4.16 for different values of the aspect ratio of the matrix elliptic criterion 98

Figure 4 . 15 :Chapter 4 . 2 b 2 .

 415422 Figure 4.15: PC model with ρ = 10%. Lower (a) and upper (b) bound homogenized strength properties in planes of Σ m = const: the different values of Σ m = 0, 0.5, 0.7, and0.9 of the maximum macroscopic isotropic tensile strength Σ max m .

Figure 4 . 16 :

 416 Figure 4.16: PC model with ρ = 40%. Lower (a) and upper (b) bound homogenized strength properties in planes of Σ m = const: the different values of Σ m = 0, 0.5, 0.7, and0.9 of the maximum macroscopic isotropic tensile strength Σ max m .

Figure 4 . 17 :

 417 Figure 4.17: FCC model with ρ = 10%. Lower (a) and upper (b) bound homogenized strength properties in planes of Σ m = const: the different values of Σ m = 0, 0.5, 0.7, and0.9 of the maximum macroscopic isotropic tensile strength Σ max m .

Figure 4 . 18 :

 418 Figure 4.18: FCC model with ρ = 40%. Lower (a) and upper (b) bound homogenized strength properties in planes of Σ m = const: the different values of Σ m = 0, 0.5, 0.7, and0.9 of the maximum macroscopic isotropic tensile strength Σ max m .

  (a) Macroscopic strain domain G hom low in the space (Σ 11 , Σ 22 , Σ 33 ).

  Macroscopic strain domain G hom up in the space (Σ 11 , Σ 22 , Σ 33 ).

Figure B. 1 : 2 b 2 = 1 .

 1221 Figure B.1: PC model with ρ = 10%, a 2 b 2 = 1. Lower (a) and upper (b) bound homogenized strength and comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (c) and the plane of the equation Σ 11 + Σ 22 = 0 (d).

B. 2 .

 2 Second configuration, FCC volume (a) Macroscopic strain domain G hom low in the space (Σ 11 , Σ 22 , Σ 33 ).

  Macroscopic strain domain G hom up in the space (Σ 11 , Σ 22 , Σ 33 ).

Figure B. 2 : 2 b 2

 222 Figure B.2: PC model with ρ = 10%, a 2 b 2 = 10. Lower (a) and upper (b) bound homogenized strength and comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (c) and the plane of the equation Σ 11 + Σ 22 = 0 (d).

  (a) Macroscopic strain domain G hom low in the space (Σ 11 , Σ 22 , Σ 33 ).

  Macroscopic strain domain G hom up in the space (Σ 11 , Σ 22 , Σ 33 ).

Figure B. 3 : 2 b 2

 322 Figure B.3: PC model with ρ = 10%, a 2 b 2 = 100. Lower (a) and upper (b) bound homogenized strength and comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (c) and the plane of the equation Σ 11 + Σ 22 = 0 (d).

B. 2 .

 2 Second configuration, FCC volume (a) Macroscopic strain domain G hom low in the space (Σ 11 , Σ 22 , Σ 33 ).

  Macroscopic strain domain G hom up in the space (Σ 11 , Σ 22 , Σ 33 ).

Figure B. 4 : 2 b 2 = 1 .

 4221 Figure B.4: PC model with ρ = 40%, a 2 b 2 = 1. Lower (a) and upper (b) bound homogenized strength and comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (c) and the plane of the equation Σ 11 + Σ 22 = 0 (d).

  (a) Macroscopic strain domain G hom low in the space (Σ 11 , Σ 22 , Σ 33 ).

  Macroscopic strain domain G hom up in the space (Σ 11 , Σ 22 , Σ 33 ).

Figure B. 5 : 2 b 2

 522 Figure B.5: PC model with ρ = 40%, a 2 b 2 = 10. Lower (a) and upper (b) bound homogenized strength and comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (c) and the plane of the equation Σ 11 + Σ 22 = 0 (d).

B. 2 .

 2 Second configuration, FCC volume (a) Macroscopic strain domain G hom low in the space (Σ 11 , Σ 22 , Σ 33 ).

  Macroscopic strain domain G hom up in the space (Σ 11 , Σ 22 , Σ 33 ).

Figure B. 6 : 2 b 2

 622 Figure B.6: PC model with ρ = 40%, a 2 b 2 = 100. Lower (a) and upper (b) bound homogenized strength and comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (c) and the plane of the equation Σ 11 + Σ 22 = 0 (d).

  (a) Macroscopic strain domain G hom low in the space (Σ 11 , Σ 22 , Σ 33 ).

  Macroscopic strain domain G hom up in the space (Σ 11 , Σ 22 , Σ 33 ).

Figure B. 7 : 2 b 2 = 1 .

 7221 Figure B.7: FCC model with ρ = 10%, a 2 b 2 = 1. Lower (a) and upper (b) bound homogenized strength and comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (c) and the plane of the equation Σ 11 + Σ 22 = 0 (d).

B. 2 .

 2 Second configuration, FCC volume (a) Macroscopic strain domain G hom low in the space (Σ 11 , Σ 22 , Σ 33 ).

  Macroscopic strain domain G hom up in the space (Σ 11 , Σ 22 , Σ 33 ).

Figure B. 8 : 2 b 2

 822 Figure B.8: FCC model with ρ = 10%, a 2 b 2 = 10. Lower (a) and upper (b) bound homogenized strength and comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (c) and the plane of the equation Σ 11 + Σ 22 = 0 (d).

  (a) Macroscopic strain domain G hom low in the space (Σ 11 , Σ 22 , Σ 33 ).

  Macroscopic strain domain G hom up in the space (Σ 11 , Σ 22 , Σ 33 ).

Figure B. 9 : 2 b 2

 922 Figure B.9: FCC model with ρ = 10%, a 2 b 2 = 100. Lower (a) and upper (b) bound homogenized strength and comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (c) and the plane of the equation Σ 11 + Σ 22 = 0 (d).

B. 2 .

 2 Second configuration, FCC volume (a) Macroscopic strain domain G hom low in the space (Σ 11 , Σ 22 , Σ 33 ).

  Macroscopic strain domain G hom up in the space (Σ 11 , Σ 22 , Σ 33 ).

Figure B. 10 : 2 b 2 = 1 .

 10221 Figure B.10: FCC model with ρ = 40%, a 2 b 2 = 1. Lower (a) and upper (b) bound homogenized strength and comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (c) and the plane of the equation Σ 11 + Σ 22 = 0 (d).

  (a) Macroscopic strain domain G hom low in the space (Σ 11 , Σ 22 , Σ 33 ).

  Macroscopic strain domain G hom up in the space (Σ 11 , Σ 22 , Σ 33 ).

Figure B. 11 : 2 b 2 B. 2 .

 11222 Figure B.11: FCC model with ρ = 40%, a 2 b 2 = 10. Lower (a) and upper (b) bound homogenized strength and comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (c) and the plane of the equation Σ 11 + Σ 22 = 0 (d).

  Macroscopic strain domain G hom up in the space (Σ 11 , Σ 22 , Σ 33 ).

Figure B. 12 : 2 b 2

 1222 Figure B.12: FCC model with ρ = 40%, a 2 b 2 = 100. Lower (a) and upper (b) bound homogenized strength and comparison of the strength domain boundary cuts including the deviatoric plane at Σ m = 0 (c) and the plane of the equation Σ 11 + Σ 22 = 0 (d).

  ).

		Chapter 3
	Limit analysis of heterogeneous
	materials with an elliptic resistance
	criterion matrix based on the
	homogenization approach
	*	*
	*

Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2 Limit analysis using the periodic strength homogenization approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2.1 Strength properties description . . . . . . . . . . . . . . . . . . . 48 3.2.2 Homogenization of strength properties . . . . . . . . . . . . . . . 50 3.2.3 Static and kinematic approaches to the homogenized strength domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.3 3D numerical formulation of the auxiliary problems . . . . . . 59 3.3.1 Finite element implementation of the static approach . . . . . . . 59 3.3.2 Finite element implementation of the kinematic approach . . . . 66 3.3.3 Formulation as a SOCP problem . . . . . . . . . . . . . . . . . . 73 3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 3.1 Introduction

  4.21) All simulations were performed on a 2.4 GHz Core i7 PC with 8GB RAM under Microsoft windows 10. The results of evaluation are shown in table 4.1.

	Clearly from the obtained results, on the one hand, it can be concluded that
		Number of iterations		CPU time (s)		Value
	Model mesh	Coarse Medium	Fine	Coarse Medium Fine	Coarse Medium Fine
	λ stat	10	9	9	2.23	28.97	101.80	9.000	9.000	8.999
	Π kine using V con	13	12	15	2.13	90.25	309.70	3.497	3.511	3.516
	Π kine using V dis	31	39	56	18.56	871.31	12302	3.497	3.511	3.516
	TABLEAU 4.1: Mesh density effect on numerical calculations

•

  In case of high inclusion density, ρ ≈ 40%The approximated strength domains G hom are presented in (Σ 11 , Σ 22 , Σ 33 ) space (figure 4.12) and in the deviatoric plane (figure4.13) for different values of Σ m : Σ m = 0, 0.5, 0.7 and 0.9 of the maximum macroscopic isotropic tensile strength Σ max m .

	Chapter 4. Numerical estimates of the Callovo-Oxfordian claystone
														strength domains
					(a)								(b)
	6							Lower Bound							Lower Bound
								Elasto-Plastic						Elasto-Plastic
								Upper Bound	10						Upper Bound
	4
	2								5
	0								0
	-2								-5
	-4
									-10
	-6
	-8	-6	-4	-2	0	2	4	6	8	-15	-10	-5	0	5	10	15

(a) (b)

) , the direction of the macroscopic strain rate D will be changed. The result of each simulation is the approximated macroscopic support function Π app up (

D

) .

• Definition of a hyper-plane Σ :

in the macroscopic stress states space whose normal is D and whose distance to the origin (or free stress state) is Π app up (

D

) .

• The interior convex hull of the hyper-planes obtained for all studied directions D then provides an exterior approach to the boundary ∂G hom of the homogenized strength domain G hom .

Under consideration of the symmetry property of the strength domain and the symmetry property of the elliptic criterion and the morphologies, following equation (4.17):

3D-FEM tool development

In this section, the proposed numerical method (presented in Chapter 3, section 3.3) is applied to evaluate the macroscopic strength domain. Three-dimensional microstructures are considered and each model is solved first the using static approach (with TET4-element) and then the kinematic approach, in both TET10-elements combined with discontinuities and continuous velocity fields, respectively.

In the static approach, with the number of elements n e and the number of nodes of each element n node , the system is solved with

× n e number of optimization variables and n node × n e of SCOP cones.

The system in the kinematic approach is solved with = +∞ (more details can be founded in [START_REF] Hill | A theory of the yielding and plastic flow of anisotropic metals[END_REF]):

As a particular case of (4.2), Equation (A.1) can be also reformulated as an elliptic resistance criterion with the material coefficient matrix M V in the form of (4.5) with:

Note that, if σ s is uniaxial yield stress of the material then k s = 1 2σ 2 s , we get:

The support function of the Von Mises criterion can be founded in [START_REF] Bignonnet | Fourier-based strength homogenization of porous media[END_REF]:

Where d υ = tr ( d ) = 0, pertinent strain rates must be incompressible.

The the function π relative to a discontinuity of velocity as follows:

Appendix A. Yield criteria formulated as ecliptic resistance criteria and support functions π

.n = 0 +∞ else (A.5)

A.2 Hill's yield criterion

Hill's yield criterion is frequently used for anisotropic materials and it takes the form:

) -1 0 (A.6) The ellipsoid equation (3.6) can be extended of using to define Hill's yield criterion by setting matrix M H with:

F , G , H , L , M , N are constants characteristic of the current state of anisotropy. If X , Y , Z are the tensile yield stresses in the principal anisotropic direction, we have:

and if R, S, T are the yield stresses in shear with respect to the principal axes of anisotropy, then:

Appendix B

Lower bound and upper bound estimation comparisons

The comparison of lower bound and upper bound properties of the Green elliptic criterion matrix reinforced with inclusions (section 4.