Keywords: Micromécanique, Tenseur de contribution, FEM, Élasticité, Conduction thermique, Propriétés effectives, Multiéchelle, Homogénéisation, Pore concave, Phases transversalement isotropes, RSA, RVEs périodiques Micromechanics, Contribution tensor, FEM, Elasticity, Thermal conduction, Effective properties, Multiscale, Homogenization, Concave pore, Transversely isotropic phases, RSA, Periodic RVEs

de Mori-Tanaka-Benveniste et celui de Maxwell. Les propriétés effectives des composites à pores réguliers sont également estimées à l'aide de l'approche dite champ complet par simulations numériques, puis comparées à la modélisation micromécanique. L'effet de microstructure complexes est étudié en considérant des Volumes Elémentaires Représentatifs (VERs) périodiques contenant des arrangements aléatoires des pores noyés dans des matrices isotropes transverses.

There are many who helped me along the way on this research.

With immense pleasure and deep sense of gratitude, I wish to express my sincere thanks to my supervisor Prof. Albert Giraud and co-supervisor Dr. Long Cheng. Without their invaluable guidance, endless motivation and continuous encouragement, this research would not have been successfully completed. Their brilliant, skilful supervision enriched this study higher than my expectation.

I would like to say big thank you for Prof. Igor Sevostianov and Dr. Jean-Franc ¸ois Barthélémy for their collaborative effort.

I express my sincere thanks to Prof. Djimédo Kondo, for his kind attendance of my CSI these years and thoughtful feedback , as well as the referee and support for my future post-doctoral application and as the president of the defence committee.

I like to acknowledge the support rendered by Dr. Stella Brach in several ways throughout my research work and the attendance of my CSI these years.

I would like to express my appreciation to my committee members, Prof. Sébastien Brisard, Prof. Oana Cazacu, and Dr. Kostas Danas, who so generously took time out of their schedules to attend my defense. Their insightful suggestions and valuable remarks enriched my knowledge and encouraged me for a further exploration.

From the bottom of my heart, I would like to thank Prof. Shaoqiang Tang and Prof. Chuang Liu for his invaluable guidance, endless motivation and continuous kindness.

I am grateful to all the GeoRessources research group members, my kind colleagues Ph.D. students & Postdoc, for their energy, understanding and help throughout these three years. I wish to extend my profound sense of gratitude to my parents and sister for all the sacrifices they made during my research and also providing me with moral support and encouragement whenever required.

Last but not the least, I could not have completed this dissertation without the support of my family group and friends, who provided stimulating discussions as well as happy distractions to rest my mind outside of my research.

Résumé

Les propriétés mécaniques des géomatériaux hétérogènes sont évaluées en prenant simultanément en compte l'anisotropie microstructurale ainsi que celle du matériau matriciel. A cet effet, l'anisotropie de la microstructure est représentée par la complexité de forme poreuse et/ou d'inclusion qui est considérée dans le présent travail comme concave ou convexe en portant nos attentions particulières aux pores supersphériques et supersphéroïdaux axisymétriques. Les tenseurs de concentration et de contribution sont calculés numériquement à l'aide de la méthode des élément finis (FEM), qui est ensuite utilisés au niveau de la modélisation semi-analytique pour l'objectif d'évaluer des propriétés effectives associées, telles que des réponses effectives élastiques et celles de conductivité.

Plus précisément, afin de résoudre le 2ème problème d'Eshelby [START_REF] Eshelby | Elastic inclusions and inhomogeneities[END_REF]) dans le cas d'inhomogénéités 3D et non ellipsoïdales, nous utilisons conditions aux limites adaptées récemment développées par [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF] basé sur une solution en champ lointain [START_REF] Sevostianov | Elastic fields generated by inhomogeneities: Far-field asymptotics, its shape dependence and relation to the effective elastic properties[END_REF]) pour intégrer l'anisotropie matricielle et la correction du biais induit par le caractère borné du domaine du maillage, ce qui permet d'accélérer la convergence du calcul sans sacrifier sa précision. En adoptant la technique d'homogénéisation numérique, les tenseurs de contribution compliance/résistivité sont calculés pour différentes formes de pores (attention particulière des pores supersphéroïdaux et supersphériques) noyés dans une matrice isotrope transverse. La méthode numérique proposée s'avère efficace et précise après un grand nombre d'estimations et leurs validations. Dans certains cas particuliers, ces validations s'effectuent avec des comparaisons entre les résultats analytiques et numériques disponibles dans littérature. En prenant en compte les résultats numériques obtenus pour des microstructures tridimensionnelles (3D) considérées, les tenseurs de contribution dans les deux cas d'inclusions/pores concaves indiqués ci-dessus, supersphère et supersphéroïde axisymétrique, sont développées dans les contextes des problèmes élastiques et thermiques. Notons ici que la forme d'inclusion/pore sphérique (i.e. paramètre de concavité p = 1) ainsi que celle de fissure circulaire (i.e. rapport d'aspect γ → 0), qui peuvent être considérés comme deux cas particuliers, sont également étudiés. Cela permet d'évaluer et de valider la méthode proposée dans le présent travail. De plus, dans le cadre de l'homogénéisation, une application aux géomatériaux poreux à matrice isotrope transverse, tels que les roches argileuses, est présentée pour illustrer l'impact du paramètre de concavité et celui de l'anisotropie de la matrice sur les propriétés globales à travers plusieurs schémas d'homogénéisation micromécanique, tels que l'approche basée sur l'approximation de non-interaction (i.e.NIA: Non-Interaction Approximation), schéma iii 2.9 Geometry and mesh of 1/8 spherical inclusion in the center of spherical matrix . . . . . . . . . . . . . 2.10 Numerical estimations of H i jkl for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/5 and different scale ratio a/L ∈ [4,18]. . . . . . . . . . . . . . . . . . .

2.11

Relative errors of numerical computations of H i jkl with respect to the analytical results [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/5 and different scale ratios a/L ∈ [4, 18]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 Geometry and mesh for a cubic model comprising a superspherical pore with p = 0.4. . . . . . . . . .

2.13

Geometry and mesh of cubic model comprising a superspheroidal pore with p = 0.4. . . . . . . . . .

Evolution of components

2.18

Left: ratios V se (p)/V 0 and V so (p)/V 0 functions of concavity parameter p, right: ratio V se (p)/V so (p) . .

2.19

The 6 independent components H i jkl of the tetragonal compliance contribution tensor of a superspherical pore embedded in TI matrix, as a function of concavity parameter p. Comparison between FEM results (dashed lines) and approximate relations (plain lines). . . . . . . . . . . . . . .

2.20

The 5 independent components H i jkl of the T I compliance contribution tensor of an axisymmetric superspheroidal pore embedded in TI matrix, as a function of concavity parameter p. Comparison between FEM results (dashed lines) and approximate relations (plain lines). Note that H 1212 is used to check accuracy of transverse isotropy by comparing to (H 1111 -H 1122 )/2.. . . . . . . . . . . . . . .

H

The two independent components H se

11 and H se 33 of a superspherical pore embedded in TI matrix, as a function of concavity parameter p and anisotropic degree of host matrix κ. Comparison between FEM results (star points) and approximate relations (plain surface). . . . . . . . . . . . . . . . . . . . . . .

3.9

The two independent components H so 11 and HH so 33 of a superspheroidal pore embedded in TI matrix, as a function of concavity parameter p and anisotropic degree of host matrix κ. Comparison between FEM results (star points) and approximate relations (plain surface). . . . . . . . . . . . . . . . . . . . A.3 Relative errors of numerical computations of A i jkl with respect to the analytical results [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/5 and different scale ratios a/L ∈ [4,18] A.5 Relative errors of numerical computations of H i jkl with respect to the analytical results [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] for the spherical pore γ = 1 embedded in a transversely isotropic matrix with different scale ratios A.7 Relative errors of numerical computations of H i jkl with respect to the analytical results [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/2 and different scale ratios a/L ∈ [4,18] [START_REF] Rodin | Eshelby's inclusion problem for polygons and polyhedra[END_REF] pour les polygones et les polyédriques, les fissures planes (Sevostianov and Kachanov (2002b)) ou non [START_REF] Mear | Elastic compliances of non-flat cracks[END_REF]) et les fissures planes sécantes [START_REF] Grechka | The influence of crack shape on the effective elasticity of fractured rocks[END_REF]).

Effective conductivity λ MTB

A la connaissance de l'auteur, des résultats analytiques purs ne sont disponibles que dans le cas d'inhomogénéités ellipsoïdales noyées dans une matrice isotrope, ou d'inhomogénéités ellipsoïdales alignées selon la direction d'une matrice transversale isotrope [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF]; [START_REF] Sevostianov | Effective elastic properties of matrix composites with transversely-isotropic phases[END_REF]; [START_REF] Levin | Elastic properties of inhomogeneous transversely isotropic rocks[END_REF] (2016) a comparé les deux approches pour prédire les propriétés élastiques effectives des solides avec des pores de forme régulière et irrégulière. Les images MEB de Grgic et al. (2013) ont illustré les pores concaves existant entre les grains de calcite dans le calcaire oolithique (Lavoux, France). Ce type de formes pourrait être décrit mathématiquement par un superellipsoïde [START_REF] Sevostianov | On computation of the compliance and stiffness contribution tensors of non ellipsoidal inhomogeneities[END_REF]). L'effet du facteur de concavité des supersphères et des pores concaves axisymétriques a été analysé dans les travaux de Sevostianov et al. (2008); [START_REF] Sevostianov | Generalization of Maxwell homogenization scheme for elastic material containing inhomogeneities of diverse shape[END_REF]; [START_REF] Chen | Combined effect of pores concavity and aspect ratio on the elastic properties of a porous material[END_REF]; Sevostianov et al. (2016a). Les auteurs ont complété la modélisation par éléments finis avec des approximations analytiques pour les tenseurs de contribution de conformité des pores de telles formes.

Ces résultats ont été utilisés pour calculer les propriétés élastiques globales de matériaux à pores concaves multiples dans une matrice isotrope : roche oolithique ( • soit par une approche semi-analytique qui combine les simulations numériques d'une inhomogénéité isolée dans un volume élémentaire et la méthode d'homogénéisation analytique. ( Roberts andGarboczi (2000, 2001); [START_REF] Arns | Computation of linear elastic properties from micrographic images: Methodology and agreement between theory and experiment[END_REF]; Garboczi and Douglas (2012))

• soit directement par l'approche des champs complets visant à numériquement calculer les champs de contraintes et de déformations ainsi que celui de déplacements d'une microstructure complexe représenté par un VER avec des inclusions/pores inclut dans la phase matricielle. [START_REF] Kushch | Conductivity of a periodic particle composite with transversely isotropic phases[END_REF]; [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF]; [START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF])

Cependant, leurs résultats correspondants ont eu tendance à se concentrer sur les matériaux à matrice isotrope. Un petit nombre de résultats explicites ont été réalisés sur les propriétés élastiques de matériaux hétérogènes à matrice anisotrope. A titre des exemples, les propriétés piézoélectriques de matériaux isotropes transverses contenant des fibres circulaires alignées avec les axes de symétrie de la matrice ont été calculées en utilisant diverses techniques d'homogénéisation par [START_REF] Dunn | Analysis of piezoelectric composite materials containing ellipsoidal inhomogeneities[END_REF]; [START_REF] Chen | Piezoelectric properties of multiphase fibrous composites: Some theoretical results[END_REF]; [START_REF] Tungyang | Micromechanical estimates of the overall thermoelectroelastic moduli of multiphase fibrous composites[END_REF]; [START_REF] Sevostianov | On the modeling and design of piezocomposites with prescribed properties[END_REF]). Sevostianov et al. (2005) a étudié le tenseur de contribution à la compliance pour une inhomogénéité sphéroïdale de rapport d'aspect arbitraire intégré dans un matériau transversalement isotrope. [START_REF] Levin | Elastic properties of inhomogeneous transversely isotropic rocks[END_REF] a calculé les propriétés élastodynamiques effectives de roches transversalement isotropes contenant des pores sphériques alignés et sphéroïdaux fortement aplatis. Les propriétés porothermoélastiques effectives des roches isotropes transverses telles que les mudstones, les argilites, les shales avec des inclusions ellipsoïdales arbitrairement orientées ont été étudiées dans le cadre de la théorie des médias effectifs (EMT), par [START_REF] Giraud | Effective poroelastic properties of transversely isotropic rocks-like composites with arbitrarily oriented ellipsoidal inclusions[END_REF][START_REF] Giraud | Effective porothermoelastic properties of transversely isotropic rock-like composites[END_REF]. Les recherches ont montré que le comportement mécanique de tels matériaux nécessite une caractérisation quantitative non seulement de la microstructure qui peut être complexe, mais aussi des propriétés matérielles de la matrice, qui retiennent notre attention dans nos travaux.

Objectifs

Les matériaux ciblés dans cette thèse seront en particulier les géomatériaux anisotropes (e.g. roches argileuses telles que marnes et argilites, etc). Ces roches sont présentées comme roches de stockage dans le cas des recherches liées au projet de stockage profond de déchets radioactifs (e.g.

Méthodologie

Cette thèse a été réalisée par une approche semi-analytique qui combine la méthode type champ moyen (e.g.

approche par tenseurs de contribution, etc.) basée sur la solution du problème de l'inclusion isolée/ du pore isolé, et celle en champ complet. A l'aide du logiciel Abaqus [START_REF] Smith | ABAQUS/Standard User's Manual[END_REF]), cette dernière s'effectue via les simulations numériques en adoptant des microstructures tridimensionnelles complexes avec, soit inclusion/pore isolé, soit un certain nombre des inclusions/pores noyés dans la phase matricielle. Ceci a été développé par des outils numériques basés sur l'approche RSA (Random Sequential Adsorption) [START_REF] Torquato | Random Heterogenous Materials[END_REF][START_REF] Lopez-Pamies | The nonlinear elastic response of suspensions of rigid inclusions in rubber: Ii -a simple explicit approximation for finite-concentration suspensions[END_REF][START_REF] Cheng | A homogenization model for porous ductile solids under cyclic loads comprising a matrix with isotropic and linear kinematic hardening[END_REF][START_REF] Zerhouni | Quantifying the effect of two-point correlations on the effective elasticity of specific classes of random porous materials with and without connectivity[END_REF] qui permettra de générer des VERs avec des microcavités et des inclusions de formes 3D diverses en relation avec les travaux récents de l'équipe HGM du laboratoire Georessources [START_REF] Kalo | Effective elastic moduli of a heterogeneous oolitic rock containing 3-d irregularly shaped pores[END_REF]; [START_REF] Grgic | The influence of CO 2 on the long-term chemo-mechanical behavior of an oolitic limestone[END_REF]). La solution numérique servira de référence pour valider et/ou calibrer les approches semi-analytiques simplifiées par champ moyen. [START_REF] Mura | Micromechanics of defects in solids[END_REF].

Méthodes d'homogénéisation

Cependant, à la connaissance de l'auteur, la solution analytique pure du deuxième Problème d'inhomogénéité n'est disponible que pour les matériaux contenant des inhomogénéités ellipsoïdales noyées dans une matrice isotrope, ou des inhomogénéités ellipsoïdales parallèles alignées dans la direction d'une matrice transversalement isotrope [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF], [START_REF] Sevostianov | Effective elastic properties of matrix composites with transversely-isotropic phases[END_REF], [START_REF] Levin | Elastic properties of inhomogeneous transversely isotropic rocks[END_REF] [START_REF] Sevostianov | On the shape of effective inclusion in the Maxwell homogenization scheme for anisotropic elastic composites[END_REF]Kachanov (2013, 2012)). Pour une fraction volumique plus élevée lorsque l'interaction n'est plus négligeable, des schémas micromécaniques plus avancés tels que le schéma MTB et Maxwell peuvent être utilisés.

Dans le schéma de Mori-Tanaka-Benveniste (MTB), toute inhomogénéité de matériau se comporte comme une inhomogénéité isolée dans la matrice. Les interactions sont considérées en supposant que chaque inhomogénéité est soumise à un champ extérieur constant qui coïncide avec le champ de contraintes moyen dans la matrice du VER [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF]; [START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF]). Le schéma d'homogénéisation de Maxwell (Maxwell (1873)) est peut-être la méthode la plus ancienne pour calculer explicitement les propriétés élastiques globales des matériaux hétérogènes. Elle est largement étudiée dans le cas où l'influence de l'interaction pourrait être prise en compte "collectivement" [START_REF] Drach | Numerical modeling of carbon/carbon composites with nanotextured matrix and 3d pores of irregular shapes[END_REF][START_REF] Sevostianov | On the shape of effective inclusion in the Maxwell homogenization scheme for anisotropic elastic composites[END_REF] One of the most important subjects in micromechanics is to describe the macroscopic behaviours by considering the evolution of the complex microstructures both in geometry and materials. It is especially for the porous composite geomaterials relating to the propagation of the micro-cracks and the growth of the micro-voids, which enormously affect their damage process, which generally introduces the deterioration of the mechanical properties and induces an softening precursor [START_REF] Mura | Micromechanics of defects in solids[END_REF]).

In the area of applied mechanics, the key importance is to be able to determine the contribution of microstructure to the effective/overall material properties of composite materials. A major objective of this kind of study is to find the statistical average material properties of the heterogeneous material through various homogenization methods [START_REF] Li | Introduction to micromechanics and nanomechanics[END_REF]). The modelling to study the mechanical behaviour of multiphase constituents generally starts with the explicit consideration of their heterogeneity at a length scale where matrix and inhomogeneity are clearly distinguishable, which is called the microscale in the following [START_REF] Bohm | Multi-Inclusion Unit Cell Studies of Reinforcement Stresses and Particle Failure in Discontinuously Reinforced Ductile Matrix Composites[END_REF]). The aim of micromechanics is just to relate the macroscopic behaviour of the heterogeneous media to the details of their microscopical constitution through the basic idea of homogenization. The conventional multiscale modelling of micromechanics is a special mathematical homogenization model that is usually not connected with any fixed length scale but require the classical length-level assumption. In continuum mechanics, we generally consider three scales, as illustrated in Fig. 1.1:

• the scale of structure or macroscale, of characteristic dimension D

• the representative volume element (RVE) scale L

• the characteristic size of the heterogeneities d (which called microscale)

1

Chapter 1

1.1

The behaviour of the RVE defines Equivalent Homogeneous Behaviour, or even effective behaviour. The homogenization step leads to the connection between the average stress/strain fields at micro and macro to constitute an overall behaviour of the RVE. Although various choices of a RVE's size are available in the literature, a rule is that the typical dimension of RVE (L) should be much larger than the microscopic one (d), and be much smaller than the typical dimension of the macrostructure (D): d ≪ L ≪ D.

Figure 1.1: The multiple scales ("macro-meso-micro") of homogenization [START_REF] Temizer | A numerical method for homogenization in non-linear elasticity[END_REF] There are a variety of methods for homogenization of porous solids, including micromechanical modelling schemes based on the contribution solution of individual inhomogeneity, direct simulations of representative volume elements (RVEs) of the material by using either the finite element method (FEM) or the fast Fourier transformedbased numerical scheme (FFT) , and establishment of variational bounds of Hashin-Shtrikman type (Hashin andShtrikman (1962b, 1963)). The present work discusses the first two approaches by using the FEM but the former one of micromechanical schemes is the main tool used in the present work.

Eshelby's solutions on the ellipsoidal inhomogeneity problem in an infinite isotropic matrix [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF])

pave the way to the development of a vast number of continuum micromechanics based homogenization schemes for evaluation of effective properties. Some extensions have then been proposed in literature. The most important ones might be the consideration of different forms of the inhomogeneities such as [START_REF] Rodin | Eshelby's inclusion problem for polygons and polyhedra[END_REF] for the polygons and polyhedral ones, the planar cracks (Sevostianov and Kachanov (2002b)) or not [START_REF] Mear | Elastic compliances of non-flat cracks[END_REF]) and intersecting planar cracks [START_REF] Grechka | The influence of crack shape on the effective elasticity of fractured rocks[END_REF]). To the author's knowledge, pure analytical results are only available in the case of ellipsoidal inhomogeneities embedded in an isotropic matrix, or aligned ellipsoidal inhomogeneities aligned in the direction of a transversely isotropic matrix [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF], [START_REF] Sevostianov | Effective elastic properties of matrix composites with transversely-isotropic phases[END_REF], [START_REF] Levin | Elastic properties of inhomogeneous transversely isotropic rocks[END_REF]). Due to mathematical difficulties, the application of such pure analytical homogenization to the materials with complex microstructure is shown limited especially for the treatment of the non-ellipsoidal/non-convex shape of the inhomogeneity in the homogenization procedure. In three-dimensional cases, the problem of irregular (nonellipsoidal) inhomogeneities reduces to integral equations and generally requires computational approaches (although, in some cases, solution can be obtained in the form of infinite series, see, for example [START_REF] Krasnitskii | Effect of a rigid toroidal inhomogeneity onthe elastic properties of a composite[END_REF]) . [START_REF] Böhm | Effects of particle shape on the thermoelastoplastic behavior of particle reinforced composites[END_REF]; [START_REF] Rasool | Effects of particle shape on the macroscopic and microscopic linear behaviors of particle reinforced composites[END_REF] analysed shape effects on the effective elastic and thermal properties of the composites containing randomly oriented and distributed spherical, octahedral, cubical and tetrahedral particles. [START_REF] Drach | Prediction of the effective elastic moduli of materials with irregularly-shaped pores based on the pore projected areas[END_REF] proposed to evaluate effect of pores of irregular shape on the overall elastic moduli using pore projected areas. This approach works well for prediction of the overall Young's moduli in different directions. [START_REF] Drach | Comparison of full field and single pore approaches to homogenization of linearly elastic materials with pores of regular and irregular shapes[END_REF] compared the two approaches to predict the effective elastic properties of solids with regular and irregular shaped pores. The SEM images of [START_REF] Grgic | Impact of chemical weathering on micro/macro-mechanical properties of oolithic iron ore[END_REF] illustrated the concave pores existing between calcite grains in the oolitic limestone (Lavoux, France). This kind of shapes could be mathematically described by superellipsoid [START_REF] Sevostianov | On computation of the compliance and stiffness contribution tensors of non ellipsoidal inhomogeneities[END_REF]). Effect of the concavity factor of superspheres and axisymmetric concave pores was analysed in the works of [START_REF] Sevostianov | On computation of the compliance and stiffness contribution tensors of non ellipsoidal inhomogeneities[END_REF]; [START_REF] Sevostianov | Generalization of Maxwell homogenization scheme for elastic material containing inhomogeneities of diverse shape[END_REF]; [START_REF] Chen | Combined effect of pores concavity and aspect ratio on the elastic properties of a porous material[END_REF]Sevostianov et al. (2016a). The authors supplemented finite element modelling with analytical approximations for compliance contribution tensors of pores of such shapes. These results were used to calculate overall elastic properties of materials with multiple concave pores in isotropic matrix: oolitic rock [START_REF] Kalo | Effective elastic moduli of a heterogeneous oolitic rock containing 3-d irregularly shaped pores[END_REF]) and 3 -D printed S i3N4 ceramics [START_REF] Lurie | Mechanical behavior of porous S i 3 N 4 ceramics manufactured with 3d printing technology[END_REF]). Their results have also been used by [START_REF] Sevostianov | On the compliance contribution tensor for a concave superspherical pore[END_REF] to derive approximate analytical expressions for the components of the compliance contribution tensor. This derivation have recently been corrected and extended to resistivity contribution tensor [START_REF] Chen | Evaluation of the effective elastic and conductive properties of a material containing concave pores[END_REF]; Sevostianov et al. (2016a)). The resistivity contribution tensor has been introduced by [START_REF] Kachanov | Explicit cross-property correlations for porous materials with anisotropic microstructures[END_REF] in the context of the cross property connection between elastic and conductive properties of heterogeneous materials. [START_REF] Kushch | Dipole moments, property contribution tensors and effective conductivity of anisotropic particulate composites[END_REF] developed explicit relations between resistivity contribution tensor and dipole moments. Trofimov et al. (2017b) compared predictions of overall elastic properties of composites reinforced with particles of different polyhedral shapes by direct finite element analysis approaches and micromechanical schemes. The results of the two approaches are in good agreement for volume fractions up to 30% for all studied material combinations. In the aboved mentioned works, effective properties of the concerned materials have been studied, by utilizing different approaches, for example:

• semi-analytical method that combines the theory of analytical homogenization and the technique of numerical simulations. (Roberts andGarboczi (2000, 2001); [START_REF] Arns | Computation of linear elastic properties from micrographic images: Methodology and agreement between theory and experiment[END_REF]; Garboczi and Douglas (2012))

• numerical homogenization for estimating the stress, strain and displacement fields with the complex microstructure by a RVE. [START_REF] Kushch | Conductivity of a periodic particle composite with transversely isotropic phases[END_REF]; [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF]; [START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF])

However, these studies have tended to focus on materials with isotropic matrix. A small number of explicit results has been carried out on elastic properties of heterogeneous materials with anisotropic matrix. Piezoelectric properties of transversely isotropic materials containing circular fibers aligned with the axes of symmetry of the matrix have been calculated using various homogenization techniques by [START_REF] Dunn | Analysis of piezoelectric composite materials containing ellipsoidal inhomogeneities[END_REF]; [START_REF] Chen | Piezoelectric properties of multiphase fibrous composites: Some theoretical results[END_REF]; [START_REF] Tungyang | Micromechanical estimates of the overall thermoelectroelastic moduli of multiphase fibrous composites[END_REF]; [START_REF] Sevostianov | On the modeling and design of piezocomposites with prescribed properties[END_REF]. [START_REF] Sevostianov | Effective elastic properties of matrix composites with transversely-isotropic phases[END_REF] calculated compliance contribution tensor for a spheroidal inhomogeneity of arbitrary aspect ratio embedded in a transversely-isotropic material. [START_REF] Levin | Elastic properties of inhomogeneous transversely isotropic rocks[END_REF] calculated effective elastodynamics properties of transversely isotropic rocks containing aligned spherical and strongly oblate spheroidal pores. Effective porothermoelastic properties of transversely isotropic rocks such as mudstones, argillites, shales with arbitrarily oriented ellipsoidal inclusions have been studied in the frame of Effective Media

Chapter 1

1.2

Theory (EMT), by [START_REF] Giraud | Effective poroelastic properties of transversely isotropic rocks-like composites with arbitrarily oriented ellipsoidal inclusions[END_REF][START_REF] Giraud | Effective porothermoelastic properties of transversely isotropic rock-like composites[END_REF]). Research has consistently shown that the mechanical behaviour of such materials require quantitative characterization not only for the microstructure which could be complex, but also on the material properties of the matrix, which catch our attention in our work.

Objectives

The materials concerned in this thesis will particularly be the anisotropic geomaterials (e.g. clay rocks, etc.), which are usually studied in the research areas for example, that correspond to the applications of the storage of the radioactive wastes (e.g. Fig. 1.2 Laboratory ANDRA in Meuse-Haute Marne). The effect of their complex microstructures on the mechanical properties at macroscopic scale deals with an important subject in the framework of underground storage. Concerning the clay rock, the microstructure has been studied at different scales via experimental observations [START_REF] Yven | Mineralogy, Texture and Porosity of Callovo-Oxfordian Argillites of the Meuse/Haute-Marne Region (Eastern Paris Basin)[END_REF]; [START_REF] Robinet | Minéralogie, porosité et diffusion des solutés dans l'argilite du Callovo-Oxfordien de Bure (Meuse, Haute-Marne, France) de l' échelle centimétrique à micrométrique[END_REF]). Fig. 1.3 shows a conceptual model illustrating the passage between them (i.e.

nano -micro -macro -structural scales) as well as different forms of inclusions and pores (e.g. sphere, spheroidal, irregular ones, etc) that are included and randomly distributed in its matrix phase. From the mechanical point of view, this complex microstructure introduces some anisotropic characters in the heterogeneous rocks that strongly affect their mechanical properties at the macroscopic scale. Moreover, the anisotropic response of the concerned rocks can be induced by another source due to the matrix material of the rock that depends on different orientations of the clay particles and inclusions/pores at relatively smaller scales (i.e. nano-scale, etc.). After recalling the relationship between different scales in the context of micromechanics, we will next look at the technique of homogenization that will be used in this thesis to study the effective properties of concerned geomaterials. In fact, it consists of determining the macroscopic properties by taking into account the available information of the heterogeneous material at its microscopic scale. The information is usually in relation to the corresponding microstructure that could be represented by a Representative Volume Element (RVE). This method is widely applied in the fields of mechanics and material science. However, the existing models in literature obtained by the purely analytical homogenization deal only with the "simplified" materials, such as the porous composite materials with an isotropic matrix and inclusions/pores in regular forms (i.e. convex forms). Moreover, the macroscopic models obtained from this approach have not been sufficiently applied in industry/engineering problems. These main reasons were due to the mathematical description of the complex microstructures during the theoretical homogenization process, in which even some other phenomenons such as the hydro-thermo-mechanical coupling could be observed that also depend on these complex microstructures.

This thesis can be considered as a first step to propose a homogenization model for such kind of geomaterials by overcoming the above-mentioned difficulties. The objective is hence to study, in the context of geomechanics, the macroscopic effective characters of geomaterial by considering the anisotropy introduced simultaneously by its complex microstructure and the matrix material. In this light, we are interested in the effect of the non-ellipsoidal forms of inclusions/pores (i.e. concave and/or irregular forms) as well as that of the elastic and/or thermal properties that will be also anisotropic. It is worthy to note that during the numerical simulation, the adapted boundary conditions, firstly proposed by [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF] for isotropic matrix, is extended in the context of anisotropic one for both elastic and thermal solution in our work. A major contribution is found as the sufficiently exact computation results and analysis that quantitatively illustrate the significant effects of the pore concavity and matrix anisotropy, and derive the analytical approximate expressions of effective elastic & thermal properties. More specifically, our attentions will be mainly paid to the porous materials with a transversely isotropic matrix as well as the superspherical Chapter 1

1.3 and superspheroidal pores, which is corresponding the aforementioned concave forms. They can be mathematically described as:

x a 2p + y b 2p + z c 2p = 1 (1.1) x 2 + y 2 a 2 p + z c 2p = 1 (1.2)
with p being the concavity parameter. a, b and c are the radii respectively along the e x , e y and e z directions in the Cartesian coordinates. It is worth mentioning that, as shown in Fig. 4, when p > 0.5, their shapes are defined as convex and concave if p < 0.5. The supersphere (Fig. 4a) and the axisymmetric superspheroid (Fig. 4b) both recover a sphere in the case of p = 1. While their particular cases when p → 0 are totally different: the supersphere tends to be three orthogonal needles and, the superspheroid approaches a penny crack perpendicularly crossed by a needle along the axis of symmetry e z .

Methodology

This thesis was realized by adopting a semi-analytical approach that combines the method of mean fields (i.e.

method depending on the contribution tensors and/or theoretical homogenization etc.) and that of the full fields (i.e. direct numerical estimations upon a whole RVE, and/or numerical homogenization). The latter one was carried out by adopting the three dimensional complex microstructures (i.e. RVEs) with either a single inclusion/pore or a finite number of inclusions/pores embedded in the matrix phase. The RVEs are numerically generated based on the algorithm of Random Sequential Adsorption (RSA) ( [START_REF] Torquato | Random Heterogenous Materials[END_REF][START_REF] Lopez-Pamies | The nonlinear elastic response of suspensions of rigid inclusions in rubber: Ii -a simple explicit approximation for finite-concentration suspensions[END_REF][START_REF] Cheng | A homogenization model for porous ductile solids under cyclic loads comprising a matrix with isotropic and linear kinematic hardening[END_REF][START_REF] Zerhouni | Quantifying the effect of two-point correlations on the effective elasticity of specific classes of random porous materials with and without connectivity[END_REF]) that is in relation to the recent works of the HGM team of laboratory GeoRessources [START_REF] Kalo | Effective elastic moduli of a heterogeneous oolitic rock containing 3-d irregularly shaped pores[END_REF]; [START_REF] Grgic | The influence of CO 2 on the long-term chemo-mechanical behavior of an oolitic limestone[END_REF]. The obtained numerical solutions are used to validate and/or calibrate the results obtained from the semi-analytical method.

Analytical homogenization methods

As aforementioned, Eshelby's solutions on the ellipsoidal inhomogeneity problem in an infinite isotropic matrix [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF]) have been extended in a vast number of developments and applications for evaluating the mechanical properties in the framework of the effective medium theory. The first Eshelby's problem is Inclusion (transformation) Problem [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF]), which considers an inclusion region in an infinite homogeneous linear elastic solid that is subjected to a uniform permanent deformation due to a nonelastic deformation as thermal expansion, phase transformation, etc. The prescribed deformation refers to eigenstrains without the constraint of surrounding matrix (stress-free transformation strains). It is worth to emphasize that both the inclusion and the matrix have the same elastic constants. This Eshelby's inclusion problem is of interest to solve the elastic fields both

Chapter 1 1.3
in the inclusion and in the matrix. The second one is Inhomogeneity Problem [START_REF] Eshelby | Elastic inclusions and inhomogeneities[END_REF]), in contrast with the Inclusion Problem, which has elastic modulus different from those of the infinitely surrounding linear elastic solid. The disturbed elastic fields of heterogeneous material submitted to remotely applied loading and prescribed eigenstrain are of great interests. The detailed solutions of Eshelby's problems can be found in [START_REF] Mura | Micromechanics of defects in solids[END_REF].

However, to the best of the author's knowledge, the purely analytical solutions of the second Inhomogeneity Problem are only available for materials containing ellipsoidal inhomogeneities embedded in an isotropic matrix, or parallel ellipsoidal inhomogeneities aligned in the direction of a transversely isotropic matrix [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF], [START_REF] Sevostianov | Effective elastic properties of matrix composites with transversely-isotropic phases[END_REF]; [START_REF] Levin | Elastic properties of inhomogeneous transversely isotropic rocks[END_REF]). Hence, a semi-analytical method via the micromechanical schemes is proposed and utilized in this thesis. In order to characterize the contribution of an individual pore to the effective mechanical properties of heterogeneous materials, the analytical homogenization approaches used in the present work are mainly those of the Non-interaction approximation (NIA) scheme, the Mori-Tanaka-Benveniste (MTB) and the Maxwell ones. More specifically, the NIA scheme is applicable in the case of low volume fractions, for which we suppose that the interaction between inhomogeneities in the composite could be ignored [START_REF] Sevostianov | On the shape of effective inclusion in the Maxwell homogenization scheme for anisotropic elastic composites[END_REF]Kachanov (2013, 2012)). For higher volume fractions problems, more advanced micromechanical schemes such as MTB and Maxwell schemes can be used when the interaction is no longer negligible. Concerning the MTB scheme, every inhomogeneity in the material is presented as isolated in the matrix, for which the interactions are considered by assuming that every inhomogeneity is subjected to a constant external field that is coincided with the average stress field in the matrix of the RVE [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF]; [START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF]). Maxwell homogenization scheme [START_REF] Maxwell | Treatise on Electricity and Magnetism[END_REF]) consists of explicitly calculating the overall elastic properties of heterogeneous materials. It is largely investigated in the cases when the influence of interaction could be accounted "collectively" [START_REF] Drach | Numerical modeling of carbon/carbon composites with nanotextured matrix and 3d pores of irregular shapes[END_REF][START_REF] Sevostianov | On the shape of effective inclusion in the Maxwell homogenization scheme for anisotropic elastic composites[END_REF]). Still in the Maxwell scheme, the far field induced by the presence of inhomogeneities is equivalent with respect to the one produced by a fictitious domain of certain shape with unknown effective properties.

Numerical algorithm

In contrast to the purely analytical modelling, the numerical method has less limit on the complexity of the microstructure (i.e. forms of inclusions/pores, etc.) regarding its mathematical description. Two numerical methods are adopted in this thesis, which are respectively based on the computation of the contribution given by an individual inhomogeneity and direct predictions of the RVE subjected to different mechanical and thermal solicitations. Note that these numerical computations are realized by using the software Abaqus [START_REF] Smith | ABAQUS/Standard User's Manual[END_REF]) via the Finite Element Method.

More specifically, in order to compute the contribution of an individual inhomogeneity, we make use of a recently developed method that allows redefining some adapted boundary conditions. The main advantage lies in reducing the bias of bounded domain of the RVE during the numerical computations, which has been proposed based on the far-field asymptotics [START_REF] Sevostianov | Elastic fields generated by inhomogeneities: Far-field asymptotics, its shape dependence and relation to the effective elastic properties[END_REF] of an inclusion problem. It is important to note here that in Chapter 1

1.4

the elastic cases, such a correction of boundary conditions is given as functions of the Green tensor and its gradient as dependent on the anisotropic elasticity of matrix material, which are rigorously calculated by means of the Fourier transform based integral method in particular for regularizing the singularities on the symmetric axis of the transverse isotropy. A large number of simulations have been carried out in the context of numerical homogenization to estimate the contribution tensors of the inhomogeneity, as aforementioned, by using FEM simulations in the framework of continuum mechanics. These contribution tensors can be used to calculate the extra elastic field produced by presence of the inhomogeneity in an elastic media subjected to a remotely applied linear elastic field. As already mentioned above, different linear homogenization schemes will be adopted and calibrated with some appropriate approximations.

The numerical results outlined in this work via the contribution tensors are of direct relevance to the effective properties of solids with inhomogeneities, but are obtained outside of the Eshelby's theory framework.

In the case of full field simulations, the RVE is first numerically generated by adopting the RSA algorithm that allows to randomly define the positions of the pores with the control of some physical information of the concerned microstructure, such as the porosity, the orientation of the pores, etc. This microstructure is subjected to either mechanical loadings (e.g. displacements or stresses) or thermal ones (e.g. temperature) for the study of the effective elastic response and that of the thermal conductivity, respectively. In this light, the whole numerical procedure is realized in this thesis by use of different softwares (e.g. MATLAB, ABAQUS, Python, Netgen and Fortran, etc.).

Specifically, the numerical generation of the RVEs is carried out by user-defined scripts of HGM team at laboratory GeoRessources that are implemented in MATLAB and/or Fortran language. NETGEN is adopted for the meshing process of the RVEs with appropriate type of quadratic 3D elements (e.g. C3D10 elements for elastic problem and DC3D10 ones for thermal problem) . Thanks to Abaqus/Standard and the supercomputer EXPLOR(un Ensemble de Calcul Scientifique Pour la Lorraine), recently developed at the University of Lorraine, the finite element simulations are carried out via parallel computations, for which the numerical homogenization models are implemented in Python language.

Publications

Until now, two papers have been published during this thesis:

• Du K., Cheng L., Barthélémy J.F., Sevostianov I., [START_REF] Du | Effective elastic properties of transversely isotropic materials with concave pores[END_REF]. Effective elastic properties of transversely isotropic materials with concave pores. Mechanics of Materials, Vol. 153, 103665.

• Du K., Cheng L., Barthélémy J.F., Sevostianov I., [START_REF] Du | Numerical computation of compliance contribution tensor of a concave pore embedded in a transversely isotropic matrix[END_REF]. Numerical computation of compliance contribution tensor of a concave pore embedded in a transversely isotropic matrix. International Journal of Engineering Science, Vol. 152, 103306.

whose detail scientific contents are presented in section 2.3 and 2.5.

Chapter 1

1.5

Concerning the remaining parts of this thesis, two papers have been finished and will be submitted to journals for potential publications:

• Du K., Cheng L., Barthélémy J.F., Sevostianov I., Giraud A., Adessina A. Numerical estimation of resistivity contribution tensor of a concave pore embedded in a transversely isotropic matrix

• Du K., Cheng L., Barthélémy J.F., Sevostianov I., Giraud A., Adessina A. Effective thermal properties of a composite containing transversely isotropic matrix and irregular shapes of pores

Organization of the work

This thesis is organized as follows. Chapter 1 reviews and summarizes the engineering and research background as well as the main methodologies adopted in the present work. Chapter 2 is devoted to investigating the effect of the concavity parameter of 3D non-ellipsoidal concave pores on the overall effective elastic properties in the framework of micromechanics. It begins by laying out the theoretical introduction of compliance & stiffness contribution tensors.

Namely, these tensors serve as the basic building block to evaluate effective elastic properties of heterogeneous materials. Numerical calculations with corrected boundary conditions are performed for non-ellipsoidal pores especially superspherical and axisymmetric superspheroidal ones embedded in a transversely isotropic matrix.

Based on the obtained numerical results, the compliance contribution tensors are next mathematically expressed by considering some appropriated approximations. It is finally taken into account for the semi-analytical homogenization in order to study the effective elastic properties of the concerned geomaterials (i.e. clay rocks). Chapter 3 deals with the effective thermal properties. The effects of pore shape and matrix anisotropy on the pore contribution parameters are qualitatively and quantitatively studied. After a large number of numerical computations, the resistivity contribution tensor was approximated and expressed as dependent on the concavity of the pore and the matrix anisotropy. Based on this semi-analytical development, the relating micromechanical modelling is proposed for the effective thermal conductivity prediction. Chapter 4 is devoted to the full field simulations by considering the "complex" microstructure of geomaterial, in which the numerical homogenization is applied. The numerical estimations of the corresponding effective responses are compared to the semi-analytical results obtained in the previous chapters for their assessment and validation. Finally, Chapter 5 gives a brief summary of this work as well as some potential perspectives in future research.

Chapter 2

Effective elastic properties of anisotropic solids with concave pores

As mentioned in the literature review, very little results was found on the study of mechanical properties of heterogeneous materials with concave pore. For determining the macroscopic elastic properties of the porous composite geomaterials, the first step of present part is devoted to evaluate the property contribution tensors by particularly taking into account the microstructural complexity and/or the matrix anisotropy. And then predict the overall elastic properties via micromechanical homogenization based on the properties contribution tensor (H-tensor) of single pore problem. This chapter is organized as follows. The two reference shapes, supersphere and axisymmetrical superspheroid are presented in Section 2.1. Section 2.2 begins by laying out the theory of property contribution tensors and looks at the numerical procedure with classic Hashin-type boundary condition (uniform strain or uniform stress, [START_REF] Hashin | Analysis of composite materials -a survey[END_REF]). The section 2.3 extends a recent numerical procedure of correction of Hashin-type boundary condition, in the case of an isotropic host matrix, by [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF], for the problem of an isolated inhomogeneity embedded in an infinite matrix domain, to the case of an anisotropic host matrix. This method allows to replace the reference problem with infinite matrix domain by an approximate problem with bounded matrix domain and corrected boundary condition based on a far field analysis of the Eshelby problem presented in [START_REF] Sevostianov | Elastic fields generated by inhomogeneities: Far-field asymptotics, its shape dependence and relation to the effective elastic properties[END_REF]. It is then applied in Section 2.3.4 to the numerical homogenization method for the computation of the contribution tensors. Numerical estimations respectively in the cases of isotropic matrix as well as the transversely isotropic one are then carried out and discussed. We investigate in section 2.5 the extension to transverse isotropy of approximation formula for the compliance contribution tensor of 3D pores of particular shapes previously studied on the basis of numerical results in section 2.3. The effective properties of a transversely-isotropic material containing concave pores are also discussed and illustrated on the example of porous clay matrix via different micromechanical homogenization schemes introduced in section 2.4. All of them are based on the solution for a single inhomogeneity problem and can be easily formulated for ellipsoidal inhomogeneities using Eshelby results [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Eshelby | Elastic inclusions and inhomogeneities[END_REF]; [START_REF] Wu | The effect of inclusion shape on the elastic moduli of a two-phase material[END_REF]). Note that section 2.3 and section 2.5 are two published articles.

Microstructure

In this section, we mathematically describe the microstructure containing concave pores that are common in both One purpose of this study set out to assess the effect of pore shape. For this goal, the superspherical and axisymmetric superspheroidal shapes are considered as two promising candidates to comply with the benchmarking representation of the pore concavity abovementioned, whose geometries are respectively expressed as Eq.(2.1) and Eq.(2.2) when a

= b = c. x a 2p + y b 2p + z c 2p = 1 (2.1) x 2 + y 2 a 2 p + z c 2p = 1 (2.2)
with p called the concavity factor or a parameter of deviation from ellipsoidal shape. As shown in Fig. 2.4, if p > 0.5, these shapes are convex and concave if p < 0.5. Supersphere 2.4a and axisymmetrical superspheroid 2.4b coincide with sphere in the case p = 1 but strongly differ in the limiting case p → 0: supersphere tends to three orthogonal needles along coordinates axes and superspheroid tends to a circular crack of unit radius crossed by a perpendicular needle along symmetry axis x 3 . The different shapes of related microstructures numerically generated by applying a user-defined MATLAB scripts.

A function is defined in MATLAB to generate the area of a cube that removes the superspherical or supersphroidal pore in the center. The 3D grid coordinates of this domain are returned by "meshgrid" function. A "stl" format file (information of triangular mesh) is output as the input for mesh generation by extracting the "isosurface" data from volume data.

Property contribution tensors and classically numerical estimations of elastic problem

In this section, we briefly recall the definition of the contribution tensors in the context of linear homogenization.

Numerical estimations respectively in the cases of isotropic matrix as well as the transversely isotropic one subjected to classical Hashin-type boundary conditions [START_REF] Hashin | The elastic moduli of heterogeneous materials[END_REF]) are then carried out and discussed.

Theoretical background on contribution tensors

∆ε = V I V H : σ ∞ (2.3) Chapter 2 2.2
Similarly, the extra stress ∆σ over the volume V due to the presence of the inhomogeneity subjected to remote strains ε ∞ boundary condition can be expressed as [START_REF] Sevostianov | Compliance tensor of ellipsoidal inclusion[END_REF]):

∆σ = V I V N : ε ∞ (2.4)
Following [START_REF] Sevostianov | Relations between compliances of inhomogeneities having the same shape but different elastic constants[END_REF], H and N are interrelated as:

H = -S 0 : N : S 0 , N = -C 0 : H : C 0 (2.5)
where S 0 and C 0 are respectively the compliance and stiffness tensors of the matrix material. Moreover, in the case of the ellipsoidal inhomogeneity, H and N are obtained by (Sevostianov and Kachanov (2002a)):

H = (S I -S 0 ) -1 + Q -1 N = (C I -C 0 ) -1 + P -1 (2.6)
where S I and C I are the compliance and stiffness tensors of the inhomogeneity, P and Q denote the fourth order Hill's tensors [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF] in terms of Green tensor G (see Section A.2).

Compliance contribution tensors in classic model

The RVE is introduced to describe the relatively complex microstructure of the composite material as the classical concept in homogenization strategy. It is required to be the smallest sample containing representative information of the whole composite materials. In this section, we focus on numerical procedures using FEM in a RVE to evaluate compliance contribution tensor of individual superspherical pore with the shape parameter 0.3 ≤ p ≤ 1.

Classical FEM calculation of H-tensor in isotropic cubical solid

It is proposed by [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF] a numerical procedure to calculate the H tensor for the composite medium with an isolated inhomogeneity. We are firstly interested in generating the related microstructures numerically by using the MATLAB scripts developed in HGM (Hydrogéomécanique multi-échelles) team of laboratory GeoRessources. The software NETGEN (Schöberl (1997)) is then devoted to create the 3D tetrahedral mesh (C3D10) of the obtained RVE and all finite element calculations are complied by Abaqus. Theoretically, for a dilute numerical simulation, one needs to use an extremely large domain compared to the size of inhomogeneity.

Hence, in this section, the scale of the cube is assumed to be twenty times bigger than the superspherical voids that allows to represent the infinite medium. (see Fig.

σ i j m = 1 V N l=1 σ (l) i j m • V (l) , (i, j = 1, 2, 3; m = 1, 2) (2.7)
where σ i j m is the component of volume average stress calculated from the results of the m-th load case, V and V (l) is the total and elemental volume, σ (l)

i j m
is the stress component i j at the centroid of the finite element l in m-th load case correspondingly, and N is the total number of elements in the model. N can then be calculated based on the results of numerical homogenization through:

σ i j m = C hom i jkl : ε 0 kl m , N = (C hom -C 0 )/ f 0 (2.8)
where ε 0 kl m is prescribed component of homogeneous strain and f 0 is initial porosity. Note that the prescribed strains are set to 1.0e -5 to ensure small deformations so that the volumes could be considered constant. Putting Eq.(2.8) into (2.5), the components of compliance contribution tensor H in isotropic matrix with different shapes of pore are then found. Its 3 independent components are presented in table 2.1 (with properties of matrix E = 1000MPa, ν = 0.3) for superspherical pores. In order to assess the reliability of current model, some comparisons are made in table 2.1 with Chapter 2 2.2 the solution of Gaussian Approximation given by [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF]. 2018) were also contrasted and similar. Due to the limit of the page, these comparison are not illustrated here. We can conclude by those comparisons that the numerical algorithm is assessed and validated for isotropic cases and irregular shapes.

Classical FEM calculation of H-tensor in transversely isotropic cubical solid

In this section, we discuss the influence of the concavity parameter (the coefficient of form p) on the contribution tensor H in transversely isotropic matrix. The properties of the matrix utilized for this part are: E 1 = 20.44GPa, E 3 = 11.306GPa, ν 12 = 0.1027, ν 31 = 0.1798, G 13 = 1.5851GPa. The same process of RVE modelling and mesh as built in Section 2.2.2.1 are adopted in the numerical homogenization. However, due to the anisotropy of the matrix, four simulations related to two uniaxial tension/compression and two simple shearing cases are executed. The components of the compliance contribution tensor for superspherical pores are then determined and the results are presented in the In the last line of table 2.2, the relative errors between analytical and numerical results for p = 1 show that the accuracy of the results is excellent. A more strong relationship between concave shapes and property contribution tensors has been reported in the results than the convex shapes. Due to the complexity of the void and/or inclusions, the volume size and the mesh refinement, the numerical estimations which induce in practice always consume large modelling and expensive time. A new reformulation of the classical boundary condition will be introduced in the next section in order to help reduce the computation time and improve the accuracy of the solutions. The main objective of this section is to estimate the compliance contribution tensor of the concave pore inhomogeneity surrounded by a transversely isotropic matrix. In this light, we make use of a recently developed adapted boundary conditions based Finite Elements Method to incorporate the matrix anisotropy and the correction of the bias induced by the bounded character of the mesh domain, which allows to accelerate the computation convergence without sacrificing its accuracy. The correction of the boundary conditions is given as functions of the Green tensor and its gradient as dependent on the anisotropic elasticity of the matrix material, which are rigorously calculated by means of the Fourier transform based integral method in particular for regularizing the singularities on the symmetric axis of the transverse isotropy. Simultaneously by complying with the numerical homogenization technique, the compliance contribution tensor is computed for different forms of pores (e.g. superspheroidal and superspherical ones, etc.) embedded in an transversely isotropic matrix. The proposed numerical method is shown to be efficient and accurate after several appropriate assessment and validation by comparing its predictions, in some particular cases, with analytical results and some available numerical ones. Finally, the effect of the pore concavity on the compliance contribution tensor is quantitatively illustrated.

Introduction

This section focuses on the analysis of the effect of a concave pore in a transversely-isotropic material on its overall elastic properties. For this goal we calculate compliance contribution tensors of concave pores (superspherical and superspheroidal) using adapted boundary conditions based Finite Elements Method and evaluate effect of the pore concavity. The work is motivated by multiple experimental observations on irregular character of pores shape in various natural and man-made materials. Particular attention could be payed on the SEM images of [START_REF] Grgic | The influence of CO 2 on the long-term chemo-mechanical behavior of an oolitic limestone[END_REF] showing the concave pores between the calcite grains in the oolitic limestone present (see also [START_REF] Emmanuel | Interfacial energy effects and the evolution of pore size distributions during quartz precipitation in sandstone[END_REF] for sandstones and [START_REF] Wark | Reassessment of pore shapes in microstructurally equilibrated rocks, with implications for permeability of the upper mantle[END_REF] in the case of the harzburgites), which could be described by introducing an concavity-convexity factor proposed by [START_REF] Sevostianov | On computation of the compliance and stiffness contribution tensors of non ellipsoidal inhomogeneities[END_REF].

While the pores and inhomogeneities of irregular shape are typical for materials studied by various branches of natural science, analytical modelling of the properties of materials with microstructures formed by inhomogeneities of non-ellipsoidal shape has not been well developed. The inhomogeneities are typically assumed to be ellipsoids of identical aspect ratios and analytical micromechanical approximations of effective properties are based on the classical Eshelby solution for ellipsoidal inhomogeneities [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Eshelby | Elastic inclusions and inhomogeneities[END_REF]. This approximation is largely responsible for the huge gap between methods of micromechanics and materials science applications. The reason for this lack is quite obvious: while for 2-D non-elliptical inhomogeneities the technique for evaluation of elastic fields associated with a single inhomogeneity is well developed (see [START_REF] Kachanov | Effective Moduli of Solids With Cavities of Various Shapes[END_REF]; Tsukrov andNovak (2002, 2004); Lanzoni et al.

( These results were used to calculate overall elastic properties of materials with multiple concave pores: oolitic rock [START_REF] Kalo | Effective elastic moduli of a heterogeneous oolitic rock containing 3-d irregularly shaped pores[END_REF]; [START_REF] Giraud | Micromechanical modeling of the effective elastic properties of oolitic limestone[END_REF]) and 3 -D printed S i 3 N 4 ceramics [START_REF] Lurie | Mechanical behavior of porous S i 3 N 4 ceramics manufactured with 3d printing technology[END_REF]. In all of the above mentioned works, it was assumed that the representative elementary volume is an infinite domain (similarly to the Eshelby's hypotheses), which induces in practice, especially for the numerical estimations, some expensive time consuming cost due to the volume size and the mesh refinement. This can be overcame by a recently proposed adapted boundary condition method (Adessina et al. ( 2017)) dedicated to the numerical resolution of the arbitrary shaped inhomogeneity problem. Consequently, the resulting predictions can be found to converge for a relatively small matrix domain and the process is shown to be less time consuming by holding a sufficiently accurate precision.

The corrected boundary conditions in this method depend on the elastic properties of the matrix and the method, initially formulated for isotropic matrix only, is extended in the present paper to the case of a transversely isotropic matrix.

This section is organized as follows. In Section 2.3.2, the classical Eshelby problem is reformulated for an inhomogeneity embedded in a finite transversely isotropic matrix by introducing the Green tensor based correction of boundary conditions. It is then applied in Section 2.3.4 to the numerical homogenization method for the computation of the contribution tensors. Next, we assess and validate the proposed method in Section 2.3.5 by comparing its predictions with the analytical and available numerical results to systematically justify its efficiency and accuracy with respect to the material anisotropy and the pore concavity. The whole procedure leads to some numerical estimations, as presented in Section 2.3.5, in the cases of the superspheroidal and superspherical voids planted in the transversely isotropic matrix. Particular attention should be payed to the significant combined effect of the material anisotropy and the shape of pore especially when it is concave. We finally present some concluding remarks in Section 2.3.7.

Green tensor based correction of boundary conditions

In Eshelby's footsteps [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF], let us consider an infinite domain Ω comprising a matrix surrounding an inhomogeneity E of arbitrary shape. The matrix is linear elastic of homogeneous stiffness tensor denoted by C 0 and the inhomogeneity is also linear elastic but not necessarily homogeneous. Before specifying the work to a transversely isotropic behaviour of the matrix embedding a porous concave domain, it is worth noticing that the main results of this section remain theoretically valid in the most general case of anisotropy of the matrix as well as arbitrary shape and internal heterogeneity of E. The infinite domain Ω is submitted to the Hashin-type boundary condition:

ξ(x) ∼ ∥x∥→∞ E.x (2.9)
where ξ is the displacement field at the position x and E denotes the remote homogeneous strain. The above mentioned problem is described as:

(P) unbounded                              div (σ(x)) = 0 (x ∈ Ω) σ(x) = C(x) : ε(x) (x ∈ Ω) ε = 1 2 gradξ + t gradξ (x ∈ Ω) ξ(x) = E • x (x ∈ ∂Ω) (2.10)
By introducing the polarization tensor

p = (C(x) -C 0 ) : ε(x) (2.11)
which is non-zero only in the inhomogeneity E, one has the displacement solution of Eq.(2.10) [START_REF] Sevostianov | Elastic fields generated by inhomogeneities: Far-field asymptotics, its shape dependence and relation to the effective elastic properties[END_REF]; [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF]; [START_REF] Barthélémy | Effective properties of ageing linear viscoelastic media with spheroidal inhomogeneities[END_REF]:

ξ(x) = E • x + x ′ ∈E gradG 0 (x -x ′ ) : p(x ′ ) dΩ x ′ (2.12)
where G 0 is the second-order Green tensor of the infinite matrix of elasticity C 0 . The first term in the r.h.s. of Eq.(2.12) represents the remote displacement field and the second one corresponds to the disturbance caused by the inhomogeneity.

The idea then is to derive from Eq.(2.12) a Taylor expansion of the displacement for remote values of x that could be eventually used at relatively large but finite distance (for instance at a mesh boundary in a Finite Element Chapter 2

2.3 computation) where the sole higher order of the asymptotic behaviour (i.e. E • x) may lack of accuracy, as shown in 2.7. For this purpose it can first be noticed that the following approximation is relevant when ∥x∥ ≫ ∥x ′ ∥:

G 0 (x -x ′ ) ∼ ∥x∥→∞ G 0 (x) ∀x ′ ∈ E (2.13)
Consequently, Eq.(2.12) can be recast into (see also [START_REF] Sevostianov | Elastic fields generated by inhomogeneities: Far-field asymptotics, its shape dependence and relation to the effective elastic properties[END_REF]):

ξ(x) = E • x+ | E | gradG 0 (x) : P (2.14)
where P defines the average polarisation tensor inside the inhomogeneity which is given as:

P = ⟨σ⟩ E -C 0 : ⟨ε⟩ E (2.15) with ⟨ε⟩ E = 1 | E | x ′ ∈E ε(x ′ ) dΩ x ′ (2.16)
being the average strain field of the inhomogeneity.

Correction of boundary condition for a finite domain

We focus henceforward on a finite domain D with a matrix containing an inhomogeneity E. Note once again that the following developments do not require any limitation on the material symmetry of the matrix nor on the shape or content of the inhomogeneity By taking into account Eq.(2.14), one has the so called bounded problem expressed as follows:

(P) bounded                              div (σ(x)) = 0 (D) σ(x) = C(x) : ε(x) (D) ε = 1 2 gradξ + t gradξ (D) ξ(x) = E • x+ | E | gradG 0 (x) : P (∂D)
(2.17)

It can be considered as the superposition of two elementary linear elastic problems with different boundary conditions, which can be respectively described as:

(P) E bounded                              div (σ(x)) = 0 (D) σ(x) = C(x) : ε(x) (D) ε = 1 2 gradξ + t gradξ (D) ξ(x) = E • x (∂D)
(2.18) and (P)

P bounded                              div (σ(x)) = 0 (D) σ(x) = C(x) : ε(x) (D) ε = 1 2 gradξ + t gradξ (D) ξ(x) =| E | gradG 0 (x) : P (∂D) (2.19)
By separately solving the subproblems (P) E bounded and (P) P bounded , one has:

⟨ε⟩ E = A E : E, ⟨σ⟩ E = B E : E (2.20) ⟨ε⟩ P = A p : P, ⟨σ⟩ P = B p : P (2.21)
where A E , B E , A p and B p are the strain and stress concentration tensors in the (P) E bounded and (P) P bounded problems, respectively.

Due to their linearity, the solution of the initial problem (P) bounded (i.e. Eq.(2.17)) can be obtained by the following superposition:

         ⟨ε⟩ E = A E : E + A p : P ⟨σ⟩ E = B E : E + B p : P (2.22)
Next, by inserting Eq.(2.22) into (2.15), one obtains:

P = D : E (2.23) with D = I -B p + C 0 : A p -1 : B E -C 0 : A E (2.24)
Consequently, Eq.(2.22) can be rewritten as:

         ⟨ε⟩ E = A E 0 : E, A E 0 = A E + A p : D ⟨σ⟩ E = B E 0 : E, B E 0 = B E + B p : D (2.25)
where A E 0 and B E 0 are respectively the average strain and stress concentration tensors of the bounded problem (P) bounded (see also Eq.(2.17)).

Determination of elastic compliance and stiffness contribution tensors

As recalled in [START_REF] Kachanov | Effective Moduli of Solids With Cavities of Various Shapes[END_REF]; [START_REF] Sevostianov | Compliance tensor of ellipsoidal inclusion[END_REF]; [START_REF] Sevostianov | Generalization of Maxwell homogenization scheme for elastic material containing inhomogeneities of diverse shape[END_REF], the compliance and stiffness contribution tensors denoted respectively by H E 0 and N E 0 allow to calculate the extra strain and stress induced by the presence of the inhomogeneity in a dilute situation such that

∆ε = f H E 0 : Σ, ∆σ = f N E 0 : E with f = | E | | Ω | (2.26)
where Σ is the remotely applied stress and as aforementioned, E is the remotely applied strain. Moreover, by applying the consistency laws ensuring that Σ and E are also the average stress and strain within a representative elementary volume, it is useful to notice that the extra stress can be expressed by means of the averages of stress and strain within the inhomogeneity phase in the second term of the r.h.s. of the following decomposition

Σ = C 0 : E + f ⟨σ⟩ E -C 0 : ⟨ε⟩ E (2.27)
Eq.(2.27) casts a new light on the definition of the average polarisation tensor in (2.24) and consequently on the stiffness contribution tensor N E 0 which is approximated here by no other than D introduced in Eq.(2.23) and related to the partial concentration tensors in Eq.(2.24)

N E 0 = D = I -B p + C 0 : A p -1 : B E -C 0 : A E (2.28)
In the general case of non ellipsoidal shapes, contribution and concentration tensors related to an inhomogeneity need to be calculated numerically as it is presented in this paper (see also in [START_REF] Eroshkin | On micromechanical modeling of particulate composites with inclusions of various shapes[END_REF] details on such calculations). Here whereas the average concentration tensors A E 0 and B E 0 are estimated by Eq.(2.25), the contribution ones H E 0 and N E 0 can be interrelated as:

H E 0 = -S 0 : N E 0 : S 0 , N E 0 = -C 0 : H E 0 : C 0 (2.29)

Case of an homogeneous inhomogeneity

The previous developments, leading to the construction of concentration tensors in Eq.(2.25) and contribution tensors in Eq.(2.28) and (2.29) from the partial tensors in Eq.(2.20) and (2.21), apply to an arbitrary inhomogeneity in terms of shape or content that is possible made of heterogeneous material. However, it may be interesting for practical implementation to examine how these tensors write in the case of an inhomogeneity of uniform stiffness tensor C E .

First it is clear that the partial concentration tensors are related by B E = C E : A E and B p = C E : A p . It follows that Eq.(2.28) becomes:

N E 0 = D = (C E -C 0 ) -1 -A p -1 : A E (2.30)
The relationships between the compliance and stiffness contribution tensors as expressed in Eq.(2.29) still hold as well as the concentration tensors A E 0 and B E 0 in Eq.(2.25) with D = N E 0 given by Eq.(2.30). However they can alternatively be written here:

A E 0 = (C E -C 0 ) -1 : N E 0 , B E 0 = C E : A E 0 (2.31)

Case of an ellipsoidal homogeneous inhomogeneity

The ellipsoidal homogeneous inhomogeneity is of particular interest in the present since analytical expressions of contribution and concentration tensors are available and can then further be compared to the numerical ones to validate the methodology in Section 2.3.5.1.In the particular case of an ellipsoidal inhomogeneity E embedded in an infinite matrix 0 of stiffness C 0 and compliance S 0 tensors, compliance H E 0 and stiffness N E 0 contribution tensors write (see [START_REF] Kachanov | Explicit cross-property correlations for porous materials with anisotropic microstructures[END_REF]; [START_REF] Kachanov | Micromechanics of materials, with applications[END_REF] for details) :

H E 0 = (S E -S 0 ) -1 + Q E 0 -1 , N E 0 = (C E -C 0 ) -1 + P E 0 -1
(2.32)

where P E 0 and Q E 0 denote the fourth order Hill's tensors [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF] of the inhomogeneity. Strain concentration tensor of the ellipsoidal inhomogeneity writes

A E 0 = I + P E 0 : (C E -C 0 ) -1
(2.33)

Green tensor of a transversely isotropic matrix

The practical implementation of the reasoning presented in the previous section 2.3.2 relies on the calculation of the Green tensor and its gradient in the problem of Eq.(2.19). The general anisotropic case of Green tensor is briefly recalled in A.2. In the sequel a transversely isotropic matrix is particularly considered. It is worthy to point out that the three dimensional (3D) solution of G 0 as well as its gradient gradG 0 in the case of transversely isotropy could be analytically calculated by using the results published in literatures [START_REF] Elliott | Three-dimensional stress distributions in hexagonal aeolotropic crystals[END_REF]; [START_REF] Yoo | Elastic interaction of small defects and defect clusters in hexagonal crystals[END_REF]; [START_REF] Pan | Point force solution for an infinite transversely isotropic solid[END_REF]; [START_REF] Mura | Micromechanics of defects in solids[END_REF]; [START_REF] Karapetian | Green's functions for the isotropic or transversely isotropic space containing a circular crack[END_REF]; [START_REF] Pouya | Green's function solution and displacement potentials for transformed transversely isotropic materials[END_REF][START_REF] Pouya | Ellipsoidal anisotropy in linear elasticity: Approximation models and analytical solutions[END_REF]; [START_REF] Kachanov | Micromechanics of materials, with applications[END_REF]. However, some of them might present accuracy problems. The misprint and an apparent degenerate case, as quoted by many authors, in Pan and Chou solution [START_REF] Pan | Point force solution for an infinite transversely isotropic solid[END_REF] were corrected by [START_REF] Pouya | Green's function solution and displacement potentials for transformed transversely isotropic materials[END_REF][START_REF] Pouya | Ellipsoidal anisotropy in linear elasticity: Approximation models and analytical solutions[END_REF]. It can also be observed that the gradient gradG 0 obtained by [START_REF] Pan | Point force solution for an infinite transversely isotropic solid[END_REF]; [START_REF] Pouya | Green's function solution and displacement potentials for transformed transversely isotropic materials[END_REF] are singular on the symmetry axis e 3 . In this paper, the Green functions G 0 i j as well as their gradients G 0 i j,k , are calculated on the one hand by using the solution given in [START_REF] Pouya | Green's function solution and displacement potentials for transformed transversely isotropic materials[END_REF] that is also briefly recalled in A.2.2, and on the other hand by applying a Fourier transform based integration in order to overcome the singularity problem on the symmetry axis e 3 . The latter is described in 2.3.3.1 and validated by comparison with Pouya's solution [START_REF] Pouya | Green's function solution and displacement potentials for transformed transversely isotropic materials[END_REF] 1 . Note that the Gaussian integration rule is adopted in the proposed Fourier transform integral based method that is shown to be very efficient and robust in the corresponding numerical implementations and applications. It must be emphasized that the Fourier transform based solution, for the Green tensor and its gradient is valid in the general anisotropic case including all others classes of symmetry.

Exact results on symmetry axis in transversely isotropic case obtained by using Fourier transform solution

In the particular case of transversely isotropic matrix stiffness tensor with symmetry axis x 3 , see relations (A.7-A.8-A.9)

C 0 = 6 i=1 c i E i (2.34)
One considers a position vector located on symmetry axis

x = x 3 e 3 , ∥x∥ = r = |x 3 |, v = x ∥x∥ = e 3 , u ψ = cos ψe 1 + sin ψe 2 (2.35)
Inverse of acoustic tensor writes

K -1 =           c 2 (1-cos(2ψ))+2c 5 c 5 (c 2 +c 5 ) -c 2 sin(2ψ) c 5 (c 2 +c 5 ) 0 -c 2 sin(2ψ) c 5 (c 2 +c 5 ) c 2 (1+cos(2ψ))+2c 5 c 5 (c 2 +c 5 ) 0 0 0 2 c 6           (2.36)
and corresponding Green function

G 0 i j (x = x 3 e 3 ) = 1 8 π 2 ∥x∥ 2 π 0 K -1 i j (u ψ ) dψ (2.37)
Non zero components write

G 0 11 (x = x 3 e 3 ) = G 0 22 (x = x 3 e 3 ) = 1 4 π ∥x∥ c 2 + 2c 5 c 5 (c 2 + c 5 ) , G 0 33 (x = x 3 e 3 ) = 1 2 π ∥x∥ c 6 , ∥x∥ = |x 3 | (2.38)
and it coincides with relations (A.44-A.45).

K -1 . v.C 0 .u ψ + u ψ .C 0 .v .K -1 = 2 √ 2c 3 + c 6 (c 2 + c 5 ) c 6           0 0 cos(ψ) 0 0 sin(ψ) cos(ψ) sin(ψ) 0           (2.39) Chapter 2 2.3 G 0 i j,k (x = x 3 e 3 ) = 1 8 π 2 ∥x∥ 2 2 π 0 K -1 . v.C 0 .u ψ + u ψ .C 0 .v .K -1 i j u ψ k -K -1 i j v k dψ, ∥x∥ = |x 3 | (2.40)
Non zero components write

G 0 11,3 (x = x 3 e 3 ) = G 0 22,3 (x = x 3 e 3 ) = - 1 4 π ∥x∥ 2 c 2 + 2c 5 c 5 (c 2 + c 5 ) , G 0 33,3 (x = x 3 e 3 ) = - 1 2 π ∥x∥ 2 c 6 (2.41) G 0 32,2 (x = x 3 e 3 ) = G 0 31,1 (x = x 3 e 3 ) = G 0 13,1 (x = x 3 e 3 ) = G 0 23,2 (x = x 3 e 3 ) = 1 4 π ∥x∥ 2 √ 2c 3 + c 6 (c 2 + c 5 ) c 6 (2.42)
or, in terms of usual components

G 0 11 (x = x 3 e 3 ) = 1 8 π ∥x∥ 3C 0 1111 -C 0 1122 C 0 1111 C 0 1111 -C 0 1122 , G 0 33 (x = x 3 e 3 ) = 1 4 π ∥x∥ C 0 2323 (2.43) G 0 11,3 (x = x 3 e 3 ) = - 1 8 π ∥x∥ 2 3C 0 1111 -C 0 1122 C 0 1111 C 0 1111 -C 0 1122 , G 0 33,3 (x = x 3 e 3 ) = - 1 4 π ∥x∥ 2 C 0 2323 (2.44) G 0 23,2 (x = x 3 e 3 ) = 1 8 π ∥x∥ 2 C 0 1133 + C 0 2323 C 0 1111 C 0 2323
(2.45)

Numerical framework of the compliance contribution tensor estimate

This section deals with the numerical procedure for the computation of the compliance contribution tensors H E 0 . For the sake of keeping this work focused and concise, we especially consider the inhomogeneity in the case of pore (denoted also by E) that is surrounded by a bounded domain D with a transversely isotropic matrix D M = D\E. Note again that the matrix anisotropy is defined around the axis e 3 in Cartesian system. Different shapes of the pore, in particular the concave ones, will be considered in the present work whose 3D geometries are realized by adopting a user-defined Matlab script. Moreover, as will be detailed in the following part of this section, the bounded domain will be represented either by a spherical model or a cubic one. They are meshed by utilizing the Netgen software (Schöberl, 1997) with quadratic 3D elements (C3D10) that are compatible with the Finite Elements computations via Abaqus/Standard software [START_REF] Smith | ABAQUS/Standard User's Manual[END_REF]. Furthermore, the boundary conditions on the external surface ∂D are given by Eq.(2.17) that could be decoupled in Eqs.(2.18) and (2.19). As a consequence, the numerical procedure will be Chapter 2 2.3 carried out simultaneously in two parts and through 8 different numerical computations. More specifically,

• for (P) E bounded problem, different boundary strain fields are respectively defined for two simple tension loadings in the directions of e 1 (i.e. E = E 0 e 1 ⊗ e 1 ) and e 3 (i.e. E = E 0 e 3 ⊗ e 3 ), and two simple shear loadings in the planes of e 1e 2 (i.e. E = E 0 (e 1 ⊗ e 2 + e 2 ⊗ e 1 )) and e 1e 3 (i.e. E = E 0 (e 1 ⊗ e 3 + e 3 ⊗ e 1 )) with an arbitrary small constant amplitude E 0 = 1.0e -5 . The displacement field subjected on the external boundary is calculated by

ξ E i = E i j x j (2.46)
• for (P) p bounded problem, four analogical computations are realized: two simple tractions with polarization tensor P = P 0 e 1 ⊗ e 1 and P = P 0 e 3 ⊗ e 3 , and two simple shear loadings with P = P 0 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ) and P = P 0 (e 1 ⊗ e 3 + e 3 ⊗ e 1 ). A normalized value of P 0 is taken as 1GPa in the numerical computations. The displacements subjected on the external boundary are given by

ξ P i =| E | ∂G i j ∂x k P k j (2.47)
It is important to emphasize that the proposed method can be applied to any form of the bounded media with transversely isotropic matrix2 . We display in Fig. 2.8a the spherical model comprising the ellipsoidal pore represented by its 1/8 geometry, which is described as

x 1 a 2 + x 2 a 2 + x 3 c 2 = 1 s.t. x 1 ≥ 0, x 2 ≥ 0, x 3 ≥ 0 (2.48)
with a and c being the major and minor radii, respectively. It can be readily obtained that the aspect ratio γ = c/a.

Moreover, the superspherical and axisymmetric superspheroidal shapes are considered as two promising candidates to comply with the benchmarking representation of the pore concavity, whose geometries are respectively expressed as:

x 1 a 2p + x 2 a 2p + x 3 a 2p = 1 (2.49) and x 2 1 + x 2 2 a 2 p + x 3 a 2p = 1 (2.50)
where p is the concavity parameter, a and c are the radii of the principal axes. Note that both the superspherical and axisymmetric superspheroidal pores are concave when p < 0.5 and convex if p > 0.5.

Note that for the sake of efficient mesh generation and convenient numerical homogenization that will be introduced in next sections, we adopt the whole cubic model in the cases of concave pores (see Figs.2.8b-2.8c).

Numerical average method

As the considered media being defined in a finite domain (i.e. porous media), the numerical homogenization method is adopted in order to estimate the compliance contribution tensor H. Following Eq.(2.25), the strain concentration tensor A E 0 will be first computed. In the framework of homogenization, one has:

⟨ε⟩ E D = E = (1 -f ) ⟨ε⟩ E D M + f ⟨ε⟩ E E ⟨ε⟩ P D = (1 -f ) ⟨ε⟩ P D M + f ⟨ε⟩ P E (2.51) where f =| E | / | D |
is the porosity, ⟨ε⟩ E D and ⟨ε⟩ P D denote the subjected macroscopic strain respectively for the (P) E bounded and (P) P bounded problems, ⟨ε⟩ E E , ⟨ε⟩ P E , ⟨ε⟩ E D M and ⟨ε⟩ P D M are the corresponding average strains in the porous phase and those in the matrix phase.

Following the Gauss theorem, ⟨ε⟩ P D M can be computed as an integral over the external boundary ∂D:

ε i j P D = 1 2 | D | M m=1 ξ P i n j + n i ξ P j (m) S (m) (2.52)
where M is the total number of the surface elements and S (m) is the area of the m-th one.

ε i j E D M and ε i j P D M can be obtained, respectively for the (P) E bounded and (P) P bounded problems, by averaging the matrix strains of the elements as:

ε i j D M = 1 | D M | N n=1 ε i j V n (2.53)
Chapter 2

2.3

with N being the total number of the volume elements and V (m) giving the volume of the n-th one.

Having in hand the above computed quantities, the average strain field in the porous phase can be obtained as:

⟨ε⟩ E E = E -(1 -f ) ⟨ε⟩ E D M f ⟨ε⟩ P E = ⟨ε⟩ P D -(1 -f ) ⟨ε⟩ P D M f (2.54) 2.3.4.2

. Strain concentration tensor and compliance contribution tensor

In this section, the strain concentration tensor A E 0 and the compliance contribution one H will be calculated based on the numerical computation as described in Section 2.3.4.1. To this end, we aim first at computing the A E and A p that both have 7 independent non-zero components. More specifically, for each of them, the components denoted as A 1111 , A 3333 , A 1122 , A 1133 and A 3311 are calculated from 2 tension and traction loadings and those of A 1212 and A 1313 can be obtained from 2 shear ones. Note that the particular consideration of the pore inhomogeneity leads to the fact that the stress concentration tensors B E and B p both vanished. Consequently, the stiffness contribution tensor N E 0 (see also Eq.(2.30)) and the compliance contribution tensor H E 0 can be simplified as:

N E 0 = -S 0 + A p -1 : A E , H E 0 = S 0 : S 0 + A p -1 : A E : S 0 (2.55)
2.3.5. Assessment and validation of the proposed numerical procedure

In this section, the estimations of compliance contribution tensor H obtained from the proposed numerical procedure will be assessed and validated by comparison with the avalaible analytical and numerical results in literatures. In this light, we systematically consider the spheroidal pores embedded in a transversely isotropic matrix as well as the superspherical ones surrounded by an isotropic matrix to respectively justify its accuracy on the anisotropy of the matrix material and that on the concavity of the pores. More specifically, in the spheroid case, comparison will be made with respect to the analytical results proposed by [START_REF] Mura | Micromechanics of defects in solids[END_REF]. Whereas in the concave superspherical cases, due to the lack of analytical results, we compare the numerical predictions with the available (

) 2018
). The transverse isotropy of the elastic matrix material is described by the parameters as shown in Table2.3. We consider a spherically bounded model comprising an ellipsoidal pore, due to its symmetry, whose 1/8 geometry and mesh are illustrated in Fig. 2.9. The corrected boundary conditions (2.17) and the uncorrected ones (2.18) will be respectively subjected to the same mesh to study the efficiency and the accuracy of the proposed numerical procedure. Note that in order to evaluate the efficiency of the numerical method with corrected boundary conditions, the computations will be carried out by adopting different spherical models with a fixed aspect ratio γ by varying the the scale ratio a/L between the major radius of the ellipsoid a and the the radius of the spherical model L. Each independent components of the compliance contribution tensor H i jkl will be numerically computed and compared with the analytical solution [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF]. Additionarilly, we show in Fig. 2.11 the relative errors of the numerical estimations with respect to the analytical solution [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] (recalled in appendix A.2). First of all, it can be readily observed that when a/L is approximately in the range of [4,8], the relative errors obtained from the corrected boundary conditions are shown to be very small (as around the value of 10 -4 ), whereas those obtained from the classical modeling is relatively significant. It is logical and evident that the corrections of the boundary conditions allows to accelerate the numerical convergence without degrading the computation accuracy. It may be noticed that non monotonous evolutions of relative errors as function of ratio a/L may be observed for small values of a/L and components (H 3333 , H 1122 , H 1133 ), in the case of corrected boundary condition. In the paper [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF], this aspect of the curve is assumed to be probably due to the error compensations, and is interpreted as an compensation between the effect of the mesh refinement and that of the boundary condition correction. By contrast, when a/L ∈ [8, 18] , the relative errors asymptotically converge and attaining a value around 10 -3 . This can be interpreted as the fact that when the scale ratio a/L is sufficiently big, the representative bounded model tends to be an infinite one such that the correction of the boundary conditions is hence neither efficient nor useful in the numerical modeling. For completeness, we provide in section A.4 the computation results of the strain contribution tensor A in the case of γ = 1/5 as well as some supplementary ones of the compliance contribution tensor H for some other ellipsoidal pores such as γ = 1 (i.e. spherical pore), γ = 1/2 and γ = 1/10. By considering all of the above mentioned computation results, the scale ratio a/L = 8 is supposed to be fixed in the next part of this work. This is of course an approximation but will be shown as sufficiently accurate in the following numerical estimations even in the case of cubically bounded matrix comprising a concave pore.

Properties of material

E 1 (GPa) E 3 (GPa) ν 12 ν 31 G 13 (

Superspherical pore in the cubically bounded model

This section deals with the assessment and validation of the numerical procedure by paying particular attention to the concavity of the pore shape. In this light, we switch off the matrix anisotropy and consider that the pore is in a superspherical form (see also Eq.(2.49)). In the case of isotropic elasticity of the matrix material, the Young modulus and Poisson's ratio are respectively supposed as E = 1GPa, ν = 0.3. Note that the computation result of the contribution tensor does not depend on the Young modulus but on the Poisson's ratio. As aforementioned, the bounded domain is represented by a cubic model with the scale ratio a/L = 8 for the convenient mesh generation and numerical homogenization. Different geometries and meshes are realized by varying the value of the concavity p. We show in Fig. 2.12 an example in the case of p = 0.4 and its mesh that is sufficiently refined in the transition zone between the matrix and porous phases. Moreover, the bounded cubic model is subjected only to the corrected boundary conditions. The predictions are compared with the FEM results that have been published in [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF]. Chapter 2 2.3

Numerical estimation in the case of concave pores

In this section, we propose to carry out the study of the compliance contribution tensor H in the case of transversely isotropic matrix comprising concave pore to understand in more detail their combined effect. The superspheroidal and superspherical shapes of the pores will be respectively considered by varying the concavity parameter p in a relatively large interval such as [0.2, 5]. Again, we restrict the study, particularly in the superspheroidal case (i.e. shape verifying the symmetry of revolution), to the assumption that the directions of the symmetry between the matrix anisotropy and that of the pore are aligned on the same direction. In Figs.2.13 and 2.12 , we show the cubic geometries that comprising the superspheroidal and superspherical pores as well as the corresponding surface meshes on internal boundaries. Moreover, for the sake of prediction accuracy, the numerical computation will be carried out based on the sufficiently refined meshes, for which the number of nodes and elements is detailed in A.3. Note that we keep the material parameters fixed as previously introduced in Table .2.3.

(a) (b)

Figure 2.13: Geometry and mesh of cubic model comprising a superspheroidal pore with p = 0.4.

Superspheroidal pores

We provide in Chapter 2 2.3 Fig. 2.15 shows the evolution of the H i jkl components on the concavity parameter p. As same as the superspheroidal cases, it can also be finally concluded that the concavity of the superspherical pore significantly affect its compliance contribution tensor H.

-1 -0. For completeness, the numerical evaluation of the strain concentration tensor A are reported in A.5.

Concluding remarks

In this section, we have numerically evaluated the compliance contribution tensor of the concave pore inhomogeneity embedded in a transversely isotropic matrix. This has been realized by use of an original developed numerical method complying with the adapted boundary conditions based method recently formulated by [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF]. The proposed numerical procedure was carried out for an arbitrarily bounded representative elementary volume and is shown to be efficient and accurate in the numerical modeling. By paying particular attentions to the pore concavity and the matrix anisotropy, a major contribution of this work is found as the sufficiently exact computation results and analysis that illustrate the significant effect of the pore concavity on the elastic properties of the matrix.

Specifically, the adapted boundary conditions based method was extended in the context of the matrix anisotropy thanks to the Green functions and their gradients applied in the correction of the boundary conditions, which have been rigorously reformulated via the Fourier transform based integral method in particular by solving the singularity Chapter 2

2.3

problem on the axis of the transverse isotropy. Moreover, the numerical homogenization method has been utilized in the proposed numerical procedure that has been firstly assessed and validated by comparing its predictions with the analytical and existing numerical results in particular cases. It is then used to the investigate the effect of the concave pore on the transversely isotropic matrix from the quantitative estimates of the compliance contribution tensor, which was found to be of critical importance especially in the case of the concavity being less than 0.5.

Last but not least, the proposed numerical method is able to deal with any general anisotropy of the matrix material but such a study has not been attempted here for the sake of keeping the work focused and concise. In the perspective point of view, effective properties such as those predicted from semi-analytical homogenization models could be developed based on the obtained numerical computations. Alternatively, new formulations of the "adapted boundary conditions" as well as their applications in the context of non-linear homogenization may also constitute a challenging extension.

Micromechanical modelling schemes based on individual inhomogeneities

To characterize contribution of an individual pore to the effective elastic properties of heterogeneous materials, the compliance contribution tensor H tensor is introduced [START_REF] Sevostianov | On the effect of interactions of inhomogeneities on the overall elastic and conductive properties[END_REF]; [START_REF] Kachanov | Micromechanics of materials, with applications[END_REF]).

The extra strain or stress induced by the presence of inhomogeneities in a RVE can be evaluated using the contribution tensor. The effective compliance tensor of a composite with inhomogeneities is then presented as:

S = S 0 + H RV E (2.56)
where S 0 is the compliance tensor of the matrix material and H RV E is the entire contribution tensor of RVE. With different assumptions, the implementation of H RV E is different.

Non-interaction approximation

For dilute distribution of inhomogeneities where interaction between inhomogeneities in the composite is ignored, each inhomogeneities can be supposed to be applied by the same remotely applied stress or strain. The non-interaction approximation (NIA) can then be used in low volume fraction due to its simplicity ( [START_REF] Sevostianov | On the shape of effective inclusion in the Maxwell homogenization scheme for anisotropic elastic composites[END_REF]Kachanov (2013, 2012))), and H RV E is found by direct summation of contributions from total inhomogeneities in the RVE:

H NIA RV E = H (i) (2.57)
where H (i) is the compliance contribution tensor of the i-th inhomogeneities [START_REF] Kachanov | Effective Moduli of Solids With Cavities of Various Shapes[END_REF]; [START_REF] Eroshkin | On micromechanical modeling of particulate composites with inclusions of various shapes[END_REF]).

Mori-Tanaka-Benveniste scheme

For higher volume fraction when the interaction is no longer negligible, more advanced micromechanical schemes can be used. In Mori-Tanaka-Benveniste scheme (MTB), every inhomogeneity in material behaves as isolated one in the matrix. The interactions are considered by assuming that every inhomogeneity is subjected to a constant external field that is coincided with the average stress field in the matrix of the RVE [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF]; [START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF]). As shown in [START_REF] Eroshkin | On micromechanical modeling of particulate composites with inclusions of various shapes[END_REF], the predictions for the effective elastic compliance tensor via the MTB scheme can be easily obtained when the NIA is given:

H MT B RV E = H NIA RV E : f (S I -S 0 ) + H NIA RV E -1 : (S I -S 0 ) (2.58) Chapter 2 2.4
where f is the volume fraction of inhomogeneities and S I is therir compliance tensor. For the limiting case when the inhomogeneities are pores, the corresponding contribution tensor turns into:

H MT B RV E = H NIA RV E 1 -f (2.59)
2.4.3. Maxwell scheme

Maxwell homogenization scheme [START_REF] Maxwell | Treatise on Electricity and Magnetism[END_REF]) may be the oldest method to explicitly calculate the overall elastic properties of heterogeneous materials. It is largely investigated in the case where the influence of interaction could be accounted "collectively" [START_REF] Drach | Numerical modeling of carbon/carbon composites with nanotextured matrix and 3d pores of irregular shapes[END_REF]; [START_REF] Sevostianov | On the shape of effective inclusion in the Maxwell homogenization scheme for anisotropic elastic composites[END_REF]). In Maxwell scheme, the far field induced by the presence of inhomogeneities is equated with the far field produced by a fictitious domain of certain shape with unknown effective properties. [START_REF] Sevostianov | Generalization of Maxwell homogenization scheme for elastic material containing inhomogeneities of diverse shape[END_REF]; [START_REF] Sevostianov | Maxwell's methodology of estimating effective properties: Alive and well[END_REF] rewrite its expressions on the basis of the compliance and stiffness contribution tensors:

H Maxwell RV E = H NIA RV E -1 -Q Ω -1 (2.60)
where Q Ω is Hill's tensor [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF]) for the effective inclusion of chosen shape Ω.

Elastic properties of transversely isotropic materials with concave pores Du et al. (2021)

The aim of this section is to extend recent works devoted to the study of the effect of 3D pores of concave shape embedded in isotropic matrix to the case of transversely isotropic (TI) matrix. In the first part of the section, approximate relations for the compliance contribution tensor of pores of two reference shapes, supersphere and axisymmetrical superspheroid, are developed on the basis of 3D Finite Element Modelling in 2.3, and known exact solutions for the limiting cases of spherical pores and circular crack. In the second part, application to effective elastic coefficients of transversely isotropic materials such as clay rocks, in the frame of homogenization theory is presented to illustrate the impact of concavity parameter on overall properties.

Introduction

In the present section, we analyse the effect of the concavity of pores on the overall elastic properties of a porous material with transversely-isotropic solid phase. For this goal, we use two homogenization techniques: Mori-Tanaka-Benveniste scheme and Maxwell scheme. Both of them are based on the solution for a single inhomogeneity problem and can be easily formulated for ellipsoidal inhomogeneities using Eshelby results [START_REF] Eshelby | Elastic inclusions and inhomogeneities[END_REF]. Nonellipsoidal shapes of the inhomogeneities are not so well studied and most of the results are obtained in 2 -D by conformal mapping [START_REF] Kachanov | Effective Moduli of Solids With Cavities of Various Shapes[END_REF]. For three-dimensional case, the problem of irregular (non-ellipsoidal) inhomogeneities reduces to integral equations and generally requires computational approaches although, in some cases, solution can be obtained in the form of infinite series, see, for example [START_REF] Krasnitskii | Effect of a rigid toroidal inhomogeneity onthe elastic properties of a composite[END_REF] . They can be generally subdivided onto two groups: (i) direct computation of stress and strain fields for a given (deterministic)

microstructure by discretizing the domain and using the FEM, and then post-processing the averages of the stress and strain fields (see, for example, [START_REF] Garboczi | Elastic moduli of composites containing a low concentration of complex-shaped particles having a general property contrast with the matrix[END_REF]) and (ii) computation of the contribution of one isolated inhomogeneity into the effective elastic properties as a function of its shape. The latter results constitute basic building blocks for theoretical models that cover diverse orientation distributions and concentrations of inhomogeneities. [START_REF] Böhm | Effects of particle shape on the thermoelastoplastic behavior of particle reinforced composites[END_REF] analysed shape effects on the effective elastic and thermal properties of the composites containing randomly oriented and distributed spherical, octahedral, cubical and tetrahedral particles. [START_REF] Drach | Prediction of the effective elastic moduli of materials with irregularly-shaped pores based on the pore projected areas[END_REF] proposed to evaluate effect of pores of irregular shape on the overall elastic moduli using pore projected areas.

This approach works well for prediction of the overall Young's moduli in different directions. [START_REF] Drach | Comparison of full field and single pore approaches to homogenization of linearly elastic materials with pores of regular and irregular shapes[END_REF] performed comprehensive numerical analysis of the pore shape on the overall properties of solids with porosity levels up to 25%. Trofimov et al. (2017b) compared predictions of overall elastic properties of composites reinforced with particles of a different polyhedral shapes by FEM and micromechanical schemes. The results of the two approaches are in good agreement for volume fractions up to 30% for all studied material combinations. The inverse problem -design of material microstructure has been has been done by [START_REF] Zohdi | Genetic design of solids possessing a random-particulate microstructure[END_REF], who determine optimal geometrical and mechanical properties of inhomogeneities for prescribed overall elastic moduli.

Effect of the concavity factor of superspheres and axisymmetric concave pores was analyzed in the works of Chen used approach developed by [START_REF] Guerrero | On a possible approximation of changes in elastic properties of a transversely isotropic material due to an arbitrarily oriented crack[END_REF] to calculate overall properties of a transversely-isotropic material containing arbitrarily oriented cracks. This result was used by [START_REF] Seyedkavoosi | Modeling of the overall elastic behavior of a transversely isotropic material reinforced with arbitrarily oriented transversely isotropic platelets[END_REF] to estimate properties of a transversely-isotropic material with multiple arbitrarily oriented oblate inhomogeneities and applied to calculation of the overall properties of dentine. In the text to follow, we use the recent numerical results of Du et al.

(2020) where compliance contribution tensors of concave pores in a transversely-isotropic material are obtained. We suggest an analytical approximation of the components of these tensor using approach of [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF] and evaluate overall elastic properties of transversely-isotropic materials containing such pores. The results are illustrated by example of shale rock containing concave pores. • superspheroidal pore • axisymmetrical superspheroidal pore In what follows we will only consider ς = 1, the first shape is then a supersphere, and the second shape obtained by a rotation about symmetry axis x 3 . Supersphere and axisymmetrical superspheroid coincide with sphere in the case p = 1 but strongly differ in the limiting case p → 0:

Compliance contribution tensor of a concave pore

| x 1 a | 2p + | x 2 a | 2p + | x 3 ςa | 2p = 1 (2.61) p = 0.
x 2 1 + x 2 2 a 2 p + | x 3 ςa | 2p = 1 (2.
supersphere tends to three orthogonal needles along coordinates axes and superspheroid tends to a circular crack of unit radius crossed by a perpendicular needle along symmetry axis x 3 .

Compliance contribution tensor H E 0 of a superspherical pore aligned with the direction of a TI matrix, with symmetry axis x 3 , respects tetragonal symmetry. With three orthogonal planes of symmetry (with normal e i ) and equivalence between x 1 and x 2 axes, its 6 independent elastic coefficients are

H E 0 = 6 i=1 h i E i (2.63)
Average compliance contribution tensors H E 0 of superspherical and axisymmetrical superspheroidal pores have been numerically calculated by using 3D Finite Element Method (FEM), see [START_REF] Du | Numerical computation of compliance contribution tensor of a concave pore embedded in a transversely isotropic matrix[END_REF] considering TI matrix with elastic properties of a shale (see Tables 2.7-2.8). 2020) are summarized in Tables (2.6-2.5). In the particular case of an ellipsoidal pore E embedded in an infinite matrix 0 of stiffness C 0 and compliance S 0 tensors, compliance H E 0 and stiffness N E 0 contribution tensors are analytical and write (see [START_REF] Kachanov | Micromechanics of materials, with applications[END_REF] for details) :

H E 0 = (S E -S 0 ) -1 + Q E 0 -1 , N E 0 = (C E -C 0 ) -1 + P E 0 -1
(2.64)

where P E 0 and Q E 0 denote the fourth order Hill's tensors [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF] of the inhomogeneity, related by relation

Q E 0 = C 0 : I -P E 0 : C 0 (2.65)
Strain Hill tensor P E 0 of a spheroidal inhomogeneity aligned in a TI matrix may be found in [START_REF] Sevostianov | Effective elastic properties of matrix composites with transversely-isotropic phases[END_REF]; [START_REF] Barthélémy | Simplified approach to the derivation of the relationship between hill polarization tensors of transformed problems and applications[END_REF] and it is recalled for convenience in appendix A.6. In the next Section we approximate these results analytically and then use them to calculate overall elastic properties of transversely isotropic matrix containing multiple concave pores.

Approximation formula for compliance contribution tensor of a superspherical or axisymmetrical superspheroidal pore embedded in a transversely isotropic host matrix

We investigate in this section the extension to transverse isotropy of approximation formula for the compliance contribution tensor of 3D pores of particular shapes previously presented. We restrict the study to the following assumptions • same directions of symmetry between matrix and pore inclusion (aligned case)

• study is focused on the concavity parameter p • in the case of superspherical pore, compliance contribution tensor respects tetragonal symmetry (6 independent components) but we will consider for applications random orientation distributions in the isotropic plane x 1x 2 (x 3 denotes the symmetry axis) on the one hand, in 3D space on the other hand. Related compliance contribution tensors respects transversely isotropic symmetry.

Volume and surface area of superspherical and axisymmetrical superspheroidal pores

Approximation formula may be obtained by using basic geometric information related to the considered reference shapes, supersphere and axisymmetrical superspheroid, defined in relations (2.61-2.62), with aspect ratio ς = 1. These informations are volume, total surface area and projected areas onto planes 0x i x 3 (with i = 1, 2).. and corresponding volumes write (Γ denotes Euler Gamma function, see [START_REF] Chen | Evaluation of the effective elastic and conductive properties of a material containing concave pores[END_REF]; Sevostianov et al. (2016a);[START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF] for details)

V se (p) = 2 3 Γ 1 2p 3 p 2 Γ 3 2p , V so (p) = 4 π 3 Γ 1+2p 2p Γ 1 p Γ 3 2p (2.66)
where superscripts se and so respectively refer to supersphere and axisymmetric superspheroid. In the range 0 < p ≤ 1, supersphere and axisymmetrical superspheroid with unit semi-length a are superscribed by unit sphere of volume V 0 = 4π/3. Ratios V se (p)/V 0 , V so (p)/V 0 and V se (p)/V so (p) are presented in figure (2.18).. Supersphere and axisymmetrical superspheroid coincide with sphere in the case p = 1 but strongly differ in the limiting case p → 0:

supersphere tends to three orthogonal needles along coordinates axes and superspheroid tends to a circular crack of unit radius crossed by a perpendicular needle along symmetry axis x 3 . . Except for some particular values of concavity parameter (p = 1 4 , 1 2 , 1), the total surface area needs to be calculated by numerical integration. As in [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF], we use the surface area of the supersphere A se (p) given by [START_REF] Trott | The Mathematica Guidebook for Numerics[END_REF], and the surface area A so (p) of the axisymmetrical superspheroid is given by the single integral accounting for symmetry of revolution

A so (p) = 1 0 1 -x 2p 1 2p 1 + x -2(1+2p) 1 -x 2p 1-2p p 1 2 dx (2.67)
The projection area S proj (p) of both 3D shapes onto planes x i x 3 (i = 1, 2, x 3 denotes symmetry axis of the axisymmetrical superspheroid) writes (Beta denotes Euler Beta function) /V so V se (p)/V so (p)

S proj (p) = 2 p Beta 1 + 1 2p , 1
Figure 2.18: Left: ratios V se (p)/V 0 and V so (p)/V 0 functions of concavity parameter p, right: ratio V se (p)/V so (p)

Approximation formula for superspherical pore

We restrict this study to the range 0.2 ≤ p ≤ 1, on the basis of the non-dimensional factors proposed by [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF], we propose approximation formula

H i jkl (p) =          S proj (p)/(V se (p)) 2/3
S octa proj /(V se octa ) 2/3 f se-a i jkl (p) H octa i jkl 0.2 ≤ p < 0.5 , no sum over i and j

A se (p)/(V se (p)) 2/3 A octa se /(V se octa ) 2/3 f se-b i jkl (p) H octa i jkl 0.5 ≤ p ≤ 1 , no sum over i and j

(2.69)

where octa denotes octahedron (particular case of supersphere at p = 0.5 represents an octahedron). Functions f se-a i jkl (p) and f se-b i jkl (p) are given in A.7. A quadratic fit has been considered for components H ii j j (no sum over i and j) whereas a fourth degree polynomial has been necessary to fit shear components

Approximation formula for axisymmetrical superspheroidal pore

Similar approximation formula are proposed for components H 1111 , H 1122 , H 1133 but total surface area has been used instead of projection area, in the concave range p < 0.5.

H 11ii (p) =         
A so (p)/(V so (p)) 2/3 A octa so /(V so octa ) 2/3 f so-a 11ii (p) H octa 11ii 0.2 ≤ p < 0.5 , no sum over i, i ∈ [1, 2, 3]

A so (p)/(V so (p)) 2/3 A octa so /(V so octa ) 2/3 f so-b 11ii (p) H octa 11ii 0.5 ≤ p ≤ 1 , no sum over i, i ∈ [1, 2, 3] (2.70)
where octa corresponds to the case (p = 0.5) which is not an octahedron but a double-conical shape. Functions f so-a 11ii (p) and f so-b 11ii (p) are given in A.7. Polynomials of degree 4 have been considered for both concave and convex domains in the range 0.2 < p < 1. It may be noticed that semi-analytical approximations using the limiting cases of aligned circular crack p → 0 and sphere p → 1 may be used for components H 3333 and H 2323 (with x 3 symmetry axis of axisymmetrical superspheroid and TI matrix ). Approximate solutions writes

H 3333 (p) = V sphere V so (p) 1 -p 1 -α o H c 3333 + p -α o 1 -α o H sphere 3333 , α o = 0.19, 0.2 ≤ p ≤ 1 (2.71) H 2323 (p) =          f 2323 (p) 0.2 ≤ p < 0.5 H sphere 2323 0.5 ≤ p ≤ 1 (2.72) f 2323 (p) = A + B exp (-ω(p -p 0 )) , p 0 = 0.2 (2.73)
Constants A and B are determined by using analytical solutions for the limiting cases of circular crack and sphere, which are imposed at p = 0.2 and p = 0.5 (difference between solutions of axisymmetrical octahedron p = 0.5 and sphere are not significant for H 2323 component).

f 2323 (p = 0.2) = H c 2323 V sphere V so (0.2) , f 2323 (p = 0.5) = H sphere 2323
(2.74) and then

A + B = H c 2323 V sphere V so (0.2) , A + B exp (-0.3 ω) =

Evaluation of the effective elastic properties of materials with transversely isotropic matrices

In this section, we calculate effective elastic properties using three homogenization techniques: Non Interaction Approximation, Mori Tanaka-Benveniste and Maxwell schemes (respectively referred with superscripts NIA, MTB and MX), see [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF]; [START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF]; [START_REF] Kachanov | Micromechanics of materials, with applications[END_REF])

S NIA = S 0 + φH E 0 , S MTB = S 0 + φ 1 -φ H E 0 , S MX = S 0 + 1 φ H E 0 -1 -Q Ω 0 -1 (2.79)
where φ denotes the porosity. Q Ω 0 denotes the second Hill tensor of the effective inclusion of the Maxwell scheme, which is supposed of spheroidal shape (with aspect ratio γ Ω ) and aligned with the directions of the TI host matrix. Q Ω 0 is related to the strain Hill tensor P Ω 0 by the relation (see appendix A.6 for details):

Q Ω 0 = C 0 : I -P Ω 0 : C 0 (2.80)
For numerical examples, we use elastic constants of shale and mudstone given in Tables 2.7 and 2.8 [START_REF] Giraud | Effective porothermoelastic properties of transversely isotropic rock-like composites[END_REF]; [START_REF] Cosenza | Effect of the local clay distribution on the effective elastic properties of shales[END_REF]. In what follows, effective elastic properties of porous clay matrix at mesoscopic scale are estimated by homogenizing micropores. We do not consider solid mineral inclusions of calcite and quartz which would need to be added for the transition from mesoscopic to the macroscopic scale (the centimeter scale of standard geomechanical laboratory tests).

The porosity of clay matrix φ (denoted f I p in [START_REF] Giraud | Effective porothermoelastic properties of transversely isotropic rock-like composites[END_REF]) is comprised in the range φ ≤ 0.30 which has been considered for the sensitivity study. It must be emphasized that most of the existing homogenization results do not account for anisotropy of the host matrix attributing the overall anisotropy to the microstructure of the pore space.

This assumption is invalid for shale rock, in particular.

For geomaterials, the superspherical shape of pores is more realistic than the axisymmetric one since it approximately represents intergranular pores. We consider only a random orientation distribution of superspherical pores which does not violate the orientation of the symmetry axes of the transversely isotropic matrix. as a function of porosity φ, MTB approximation, aligned axisymmetrical superspheroidal pores. App 1 (plain lines): approximation formula (relations 2.70-2.71-2.72) for axisymmetrical spheroidal pores, App 2 (dashed lines) : approximation oblate spheroid with same volume (relations 2.78-2.83). Note that the two approximations coincide only for normal Young's modulus E MTB 3 . 0.5. It is expected as the supersphere tends to three orthogonal needles (with zero volume and zero surface) when p tends to zero, whereas the axisymmetrical superspheroid tends to a circular crack with one central orthogonal needle (the latter having zero volume but non zero surface).. As previously indicated, the most relevant shape compared to microstructure of porous materials is certainly supersphere.

• a significant anisotropic degree in the case of aligned axisymmetrical superspheroidal pores, in the limit p → 0. as a fonction of porosity φ, MTB approximation, superspherical (plain lines) and aligned axisymmetrical superspheroidal (dashed lines) pores.

Comparisons between NIA, MTB and Maxwell homogenization schemes

Effective elastic properties predicted by Maxwell, MTB, NIA are presented in figures (2.25-2.26) for respectively for axisymmetrical superspheroidal and superspherical pores randomly oriented in transverse plane. The shape of the effective inclusion of the Maxwell scheme is still an open issue when host matrix is anisotropic (see [START_REF] Sevostianov | On the shape of effective inclusion in the Maxwell homogenization scheme for anisotropic elastic composites[END_REF]; [START_REF] Giraud | Effective electrical conductivity of transversely isotropic rocks with arbitrarily oriented ellipsoidal inclusions[END_REF]). The sensitivity study on the shape of the effective inclusion, and oblate spheroid of aspect ratio γ Ω = 0.5 -1 confirms that it is a parameter of major importance.

Concluding remarks

In the present work, effective properties of a transversely-isotropic material containing concave pores are discussed and illustrated on the example of porous clay matrix. For this goal we used NIA, MTB and Maxwell homogenization schemes. All techniques require the explicit analytical representation of the compliance contribution tensor for a single pore. These tensors were calculated for the set of superspherical and axisymmetrical superspheroidal pores with concavity parameter p in the range 0.2 ≤ p ≤ 1 which covers both concave (0.2 ≤ p < 0.5) and convex shapes (0.5 ≤ p ≤ 1) using FEM. Based on the numerical solution for octahedron (p = 0.5) and analytical solution for sphere we built analytical approximations of the compliance contribution tensor in terms of the pore concavity parameter. The accuracy of this approximation is better than 5% for all the tensor components. We show that the concavity parameter p is a parameter of major importance on the overall elastic behavior of transversely isotropic materials containing such pore shapes. It is impossible to match effect of concave pores by oblate spheroidal pores except in the specific case of normal Young's modulus and axisymmetrical shape. The main novelty of this study is the account of concavity effects related to 3D shapes embedded in an anisotropic matrix while previous studies where done in the case of isotropic matrix. On the basis of presented results, it appears that it is not possible to separate effect of anisotropy from the effect of concavity to extend previous results obtained in the isotropic case. Accordingly, it is not possible to take effect of concavity from isotropic matrix and put it into TI matrix.

Chapter 3

Effective thermal properties of anisotropic solids with concave pores Heat transfer through porous media has caught increasingly interest in diverse fields. Its widespread application in biomedicine, chemistry, electronics, mechanic, geology and nuclear engineering has triggered a wave of research from different disciplines. For these reasons, our work in this chapter focuses on the determination of the effective linear thermal conductivity for anisotropic porous media. Property contribution tensors are also utilized here to compute the effective conductivity properties of random heterogeneous materials, in the framework of micromechanical modelling.

The resistivity contribution tensor has been introduced by Sevostianov and Kachanov (2002a) in the context of the cross-property connection between elastic and conductive properties of heterogeneous materials. They have also reported the fact that inclusion shapes affect the elasticity and the conductivity differently. In contrast with the work that focus on linear elasticity in chapter 2, this part derives the resistivity contribution tensor of linear conductivity problem. In this sense, the present work can be viewed as an extension of chapter 2 for the concave pores in anisotropic conductivity media. The thermal problem is of lower order which allows us to consider much more materials. That's why we conduct the thermal problem. This chapter seeks to explain the effect of pore shape and matrix anisotropy on the conductivity properties of the anisotropic porous materials (containing irregular pores and diversely anisotropic degree of matrix). Moreover, the superspherical and axisymmetric superspheroidal shapes are still two typical shapes of pores which we'll pay particular attention in this part. The corrected boundary condition, initially proposed by [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF], extended to the anisotropic case of elastic problem in section 2.3.2, are adopted to the problem of heat conduction during the numerical procedure in this chapter. The contents of this chapter are organized into three sections. The first section 3.1 reviews the basic and introductory knowledge of property contribution tensors in thermal issue. Section 3.2 outlines the procedure of simulation and the subsequent finite element (FE) realization. It then quantitatively presents the contribution of individual inhomogeneity due to its shapes and matrix anisotropy to the resistivity contribution tensor via the FE based numerical homogenization. Note that the content of this section is being submitted as a separate article. Section 3.3, in the first step, is devoted to the approximate expressions of property contribution tensor in terms of numerical results in section 3.2. The effective thermal properties of composite materials are detailed afterwards via various implementations such as NIA, MTB and Maxwell schemes with random or aligned microstructure.

Property contribution tensors of conductivity problem

Contrarily to the elastic problem and related fourth order tensors, the thermal conductivity problem and the corresponding second-order contribution tensors are considered in this chapter (which is mathematically equivalent to the electric conductivity problem). Table 3.1 summarizes the equivalences between thermal conductivity properties and elasticity. linear elasticity linear thermal conduction displacement vector ξ temperature T strain tensor thermal gradient

ε = 1 2 grad(ξ) + t grad(ξ) ε = gradT ε i j = 1 2 ξ i, j + ξ j,i ε i = T ,i
stress tensor heat flux vector Assuming a linear conduction law (linear relation between the temperature gradient ε and the remotely applied heat flux vector Σ, same goes for the heat flux σ and the remotely applied thermal gradient E), we then have a linear relation for the change due to the inhomogeneity in a dilute situation:

σ = C : ε σ = -λ • ε σ i j = C i jkl ε lk σ i = -λ i j ε j elastic stiffness tensor thermal conductivity tensor C λ σ = C : ε σ = -λ • ε elastic compliance tensor thermal resistivity tensor S = C -1 r = λ -1 ε = S : σ ε = -r • σ
∆ε = -f H E 0 • Σ, ∆σ = -f N E 0 • E with f = | E | | Ω | (3.1)
where ∆ε and ∆σ respectively are the extra thermal gradient and heat flux vector due to the presence of inhomogeneity.

The symmetric second rank tensors H E 0 and N E 0 can be called the resistivity and conductivity contribution tensors. The first segment of formulas Eqs.(3.1) can be interpreted as the change in temperature gradient that is required to maintain the same heat flux after the inhomogeneity has been introduced. The property contribution tensors then turn into conductivity N and resistivity H contribution tensors in contrast with the elasticity. The contribution tensors H E 0 and N E 0 of a given inhomogeneity of any shape are interrelated as follows [START_REF] Sevostianov | On computation of the compliance and stiffness contribution tensors of non ellipsoidal inhomogeneities[END_REF]):

H E 0 = -r 0 • N E 0 • r 0 , N E 0 = -λ 0 • H E 0 • λ 0 (3.2)
where r 0 , λ 0 are resistivity and conductivity tensor of host matrix respectively.

Computational resistivity contribution tensor of heterogeneous materials with concave pores and transversely isotropic matrix

In this section, we focus on the effect of non-ellipsoidal concave pore on thermal conduction properties of porous media with an infinite transversely isotropic matrix. This effect is described by the resistivity contribution tensors that will be computed via the Finite Elements (FE) based numerical homogenization. The FE computations will be carried out with some adapted and bounded boundary conditions that are formulated as the gradient of the Green function for the three dimensional Poisson's equation in infinite anisotropic medium. It allows to incorporate the matrix anisotropy and the correction of the bias induced by the bounded character of the mesh domain. The boundary conditions are constructed and applied in such a way that they accelerate the convergence of numerical computations, and therefore preserve the accuracy of estimations. This is proved after several appropriate assessment and validation by comparing its predictions, in some particular cases, with analytical results and some available numerical ones. Finally, the effect of the pore concavity as well as that of the matrix anisotropy on the resistivity contribution tensor are quantitatively illustrated.

This section is organized as follows. In Section 3.2.1, the classical concerned problem with an infinite transversely isotropic matrix is reformulated for an inhomogeneity embedded in a finite one by introducing the Green tensor based correction of boundary conditions. It is then applied in Section 3.2.2 to the reformulation of the contribution tensors.

Next, in section 3.2.3, a numerical framework based on the adapted boundary conditions is proposed by adopting the numerical homogenization method, which is also assessed and validated by comparing its predictions with some analytical and available numerical results to systematically justify its efficiency and accuracy with respect to the pore concavity and the material anisotropy. The whole procedure leads to some numerical estimations, as presented in Section 3.2.5, in the cases of the superspheroidal and superspherical voids planted in the transversely isotropic matrix.

Particular attention should be payed to the significant combined effect of the material anisotropy and the shape of pore especially when it is concave. We finally present some concluding remarks in Section 3.2.6.

Green function based correction of boundary conditions

We consider an infinite domain Ω comprising an inhomogeneity E of arbitrary shape surrounded by a matrix. The matrix is thermally conductive obeying the Fourier law:

σ = -λ • ε (3.3) Chapter 3 3.2
where σ, λ and ε respectively denote the heat flux vector, the thermal conductivity tensor and thermal gradient with T being the temperature field. Assuming that the infinite domain Ω is submitted to remote Hashin-type boundary condition:

T (x) ∼ ∥x∥→∞ E • x (3.4)
with T (x) being the temperature field at the position x and E denoting the remote homogeneous thermal gradient, the above mentioned problem is described as:

(P) unbounded                              div (σ(x)) = 0 (x ∈ Ω) σ(x) = -λ • ε(x) (x ∈ Ω) ε = gradT (x ∈ Ω) T (x) = E.x (x ∈ ∂Ω) (3.5)
Borrowing ideas from [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF][START_REF] Adessina | Micromechanical model for the diffusion properties of materials embedding complex structures[END_REF] (see also [START_REF] Du | Numerical computation of compliance contribution tensor of a concave pore embedded in a transversely isotropic matrix[END_REF] or section 2.3), the temperature solution of Eq.(3.5) can be calculated as:

T (x) = E • x + x ′ ∈E gradG 0 (x -x ′ ) • p(x ′ ) dΩ x ′ (3.6)
with the polarization vector [START_REF] Ammari | Polarization and Moment Tensors: With Applications to Inverse Problems and Effective Medium Theory[END_REF]):

p(x) = -λ(x) -λ 0 • ε(x) (3.7)
which is non-zero only in the inhomogeneity E, and the Green function G 0 for the three dimensional Poisson equation of the infinite medium with thermal conductivity λ 0 . The expression of G 0 as well as its gradient gradG 0 are briefly recalled in A.8. Note that the first term in the r.h.s. of Eq.(3.6) represents the remote temperature field and the second one corresponds to the disturbance caused by the inhomogeneity.

We perform an approximation by assuming that:

G 0 (x -x ′ ) ∼ ∥x∥→∞ G 0 (x) ∀x ′ ∈ E (3.8)
Eq.(3.6) can then be rewritten as:

T (x) = E • x+ | E | gradG 0 (x) • P (3.9) 66 Chapter 3 3.2
in which P is the average polarisation vector inside the inhomogeneity, reading as:

P = p E = 1 | E | x ′ ∈E p(x ′ ) dΩ x ′ = 1 | E | x ′ ∈E -λ(x ′ ) • ε(x ′ ) dΩ x ′ ⟨σ⟩ E +λ 0 • 1 | E | x ′ ∈E ε(x ′ ) dΩ x ′ ⟨ε⟩ E (3.10)
with ⟨σ⟩ E and ⟨ε⟩ E being the averages heat flux and thermal gradient vector of the inhomogeneity, respectively.

In this context, it is crucial to remark that a careful attention should be paid on the approximation given by Eq.(3.8) which is theoretically true when || x ||≫|| x ′ ||. In fact, by developing the Taylor expansion of the remote temperature at x, the higher order of the asymptotic behavior (i.e. E • x) may lack of accuracy. Nevertheless, in the perspective of the FE computations, it is convenient to define an appropriate finite mesh scale that we will show in the next sections this approximation delivers very accurate results and thus is used here as well.

Based on Eq.(3.9), the approach of the corrected boundary conditions consist in reformulating the Eshelby like problem Eq.(3.5) into the one on a finite domain D, which can be expediently written as:

(P) bounded                              div (σ(x)) = 0 (D) σ(x) = -λ(x) • ε(x) (D) ε = gradT (D) T (x) = E • x + |E| gradG 0 (x) • P (∂D) (3.11)
At this point, it is emphasized that the following developments do not require any limitation on the material symmetry of the matrix nor on the shape or content of the inhomogeneity. See for instance section 2.3 in the context of linear anisotropic elasticity.

Due to the linearity of Eq.(3.11), it can be considered as the superposition of two elementary linear thermal problems with different boundary conditions: one is composed by the remote Hashin-type boundary condition T (x) = E • x, denoted as the (P) E bounded problem:

(P) E bounded                              div (σ(x)) = 0 (D) σ(x) = -λ(x) • ε(x) (D) ε = gradT (D) T (x) = E • x (∂D) (3.12)
and according to the polarization field, another one depends on the gradient of the Green function gradG 0 , named as the (P) P bounded problem

(P) P bounded                              div (σ(x)) = 0 (D) σ(x) = -λ(x) • ε(x) (D) ε = gradT (D) T (x) = |E| gradG 0 (x) • P (∂D) (3.13)
A direct consequence from the resolution of Eqs.(3.12) and (3.13) is that

⟨ε⟩ E = A E • E, ⟨σ⟩ E = -B E • E (3.14) ⟨ε⟩ P = A p • P, ⟨σ⟩ P = -B p • P (3.15)
where ⟨ε⟩ E , ⟨ε⟩ P , ⟨σ⟩ E and ⟨σ⟩ P respectively denote the average thermal gradient and heat flux vectors over the inhomogeneity, and A E , A P , B E and B P are the second order concentration tensors.

         ⟨ε⟩ E = A E • E + A p • P ⟨σ⟩ E = -B E • E -B p • P (3.16)
Next, Substituting expression (3.16) in (3.10) renders:

P = D • E (3.17)
where

D = (i + B p -λ 0 • A p ) -1 : λ 0 • A E -B E (3.18)
with i being the second order identity tensor. At this stage, it is a simple matter to combine Eqs.(3.16) -(3.18) to finally establish the main result of this section:

         ⟨ε⟩ E = A E 0 • E, A E 0 = A E + A p • D ⟨σ⟩ E = -B E 0 • E, B E 0 = B E + B p • D (3.19)
where A E 0 and B E 0 denote respectively the average thermal gradient and heat flux concentration tensors of the bounded problem Eq.(3.11).

Conductivity and resistivity contribution tensors

Still in the context of linear conduction law, the extra thermal gradient (denoted by ∆ε) and heat flux vector (denoted by ∆σ) induced by the presence of inhomogeneity can be calculated in the dilute scheme as:

∆ε = -f H E 0 • Σ, ∆σ = -f N E 0 • E with f = | E | | Ω | (3.20)
where f denotes the volume fraction of the inhomogeneity and, H E 0 and N E 0 are both symmetric second rank tensors called the resistivity and conductivity contribution tensors, respectively1 . Moreover, the consistency laws ensuring that Σ and E are also the average heat flux and thermal gradient vectors within a Representative Volume Element (RVE):

Σ = (1 -f )⟨σ⟩ matrix + f ⟨σ⟩ E (3.21) E = (1 -f )⟨ε⟩ matrix + f ⟨ε⟩ E (3.22)
As matrix being homogeneous, by applying the Fourier Law, the average heat flux of the matrix can be written as:

⟨σ⟩ matrix = -λ 0 • ⟨ε⟩ matrix = 1 1 -f -λ 0 • E + f λ 0 • ⟨ε⟩ E (3.23)
Combining Eqs.(3.21) and (3.23), one has:

Σ = (1 -f )⟨σ⟩ matrix + f ⟨σ⟩ E = -λ 0 • E + f λ 0 • ⟨ε⟩ E + f ⟨σ⟩ E (3.24)
and then

Σ = -λ 0 • E + f ⟨σ⟩ E + λ 0 • ⟨ε⟩ E ∆σ (3.25)
Chapter 3

3.2

By taking into account the polarisation vector P inside the inhomogeneity E given in relation Eq.(3.10), consistency equation rewrites:

Σ = -λ 0 • E + f P ∆σ (3.26)
Comparison with relation Eq.(3.20) and Eq.(3.17) allows to obtain a new expression of the conductivity contribution tensor N E 0 which is approximated here by D uniquely:

N E 0 = -D = -(i + B p -λ 0 • A p ) -1 • λ 0 • A E -B E (3.27)
The contribution tensors H E 0 and N E 0 of a given inhomogeneity of any shape are interrelated as follows [START_REF] Sevostianov | On computation of the compliance and stiffness contribution tensors of non ellipsoidal inhomogeneities[END_REF]):

H E 0 = -r 0 • N E 0 • r 0 , N E 0 = -λ 0 • H E 0 • λ 0 (3.28)

Case of homogeneous inhomogeneity

In the case of a homogeneous material with thermal conductivity tensor λ(x) = λ E in the inhomogeneity x ∈ E, by applying the Fourier Law on average heat flux and thermal gradient ⟨σ⟩ E = -λ E • ⟨ε⟩ E , the partial concentration tensors have the following relations:

B E = λ E • A E , B P = λ E • A P (3.29)
It follows that Eq.(3.27) becomes:

N E 0 = -D = (λ E -λ 0 ) -1 + A P -1 • A E (3.30)
By using relations (3.20-3.25), one has:

∆σ = f ⟨σ⟩ E + λ 0 • ⟨ε⟩ E = -f N E 0 • E (3.31)
and consequently

⟨σ⟩ E + λ 0 • ⟨ε⟩ E = -N E 0 • E (3.32)
and then combining with the Fourier law, the average thermal gradient is derived in the form:

⟨ε⟩ E = (λ E -λ 0 ) -1 • N E 0 • E (3.33)
Finally, we get the following connection between the concentration and contribution tensors by comparing to relation (3.19):

A E 0 = (λ E -λ 0 ) -1 • N E 0 , B E 0 = λ E • A E 0 (3.34)
So that the concentration and contribution tensors can be calculated from each other.

Note that the particular consideration of the perfectly insulating inhomogeneity λ E → 0 leads to the fact that the heat flux concentration tensors B E and B P both vanish. In practical applications, pore inhomogeneities fully saturated by dry air can be approximated, in a first approach, by the limiting case of the perfectly insulating inhomogeneity, due to the low thermal conductivity of dry air compared to most of solids constituting geomaterials. Consequently, the conductivity contribution tensor N E 0 (see Eq.(3.30) and the resistivity contribution tensor H E 0 can be simplified as:

N E 0 = -r 0 + A P -1 • A E , H E 0 = r 0 • r 0 -A P -1 • A E • r 0 (3.35)
And average thermal gradient concentration tensor A E 0 may be deduced from relation Eq.(3.34)

A E 0 = -r 0 • N E 0 = r 0 • r 0 -A P -1 • A E (3.36)
In the case of an infinity conductivity inhomogeneity λ E → ∞, the thermal gradient concentration tensors A E and A P both vanish. Eq.(3.27) becomes:

N E 0 = -D = (i + B p ) -1 • B E (3.37)
and the dual resistivity contribution tensor becomes:

H E 0 = -r 0 • (i + B p ) -1 • B E • r 0 (3.38)
and the thermal gradient concentration tensor A E 0 vanishes.

Case of ellipsoidal inhomogeneity

For a general isolated ellipsoidal inclusion, the conductivity N E 0 and resistivity H E 0 contribution tensors have explicit forms in terms of Hill polarization tensors:

N E 0 = P E 0 + (λ E -λ 0 ) -1 -1 , H E 0 = Q E 0 + (r E -r 0 ) -1 -1 (3.39)
where Q E 0 and P E 0 are Hill polarization tensors. Their analytical expressions are available for inclusion in spheroidal shape (ellipsoid of revolution) embedded in a transversely isotropic matrix (see A.9).

The thermal gradient concentration tensor A E 0 is then derived:

A E 0 = i + P E 0 • (λ E -λ 0 ) -1 (3.40)
And analytical contribution tensors can be obtained by Eq.(3.34) which will be further compared to the numerical results to validate the methodology.

Numerical framework of the resistivity contribution tensor estimate

In this section, the focus is on the numerical procedure for the computation of the resistivity contribution tensors H E 0 . In order to simplify the notation, we drop the subscript 0 and the superscript E with the understanding that H i j represents the components of H E 0 . Moreover, for the sake of keeping this work focused and concise, we especially consider here the inhomogeneity in the case of individual and perfectly insulating pore (denoted also by E) implemented in the center of a cubically bounded RVE D with a transversely isotropic matrix2 D M = D\E. Different shapes of the pore, in particular the non-ellipsoidal concave ones, will be considered whose 3D geometries, as shown in Fig. 3.1, are realized by adopting a user-defined Matlab script. Note that directions of the symmetry between the matrix anisotropy and that of the pore are both aligned on the axis e 3 . Moreover, the RVE is meshed by use of the software Netgen (Schöberl, 1997) with quadratic 3D thermal elements (DC3D10) that are compatible with the Finite Elements computations via Abaqus/Standard software [START_REF] Smith | ABAQUS/Standard User's Manual[END_REF].

An example in the case of superspherical pore with p = 0.4 is shown in Fig. 3.2. As discussed in Section 3.2.2, such a (P) bounded problem (e.g. Eq.(3.11)) can be decoupled into (P) E bounded and (P) P bounded (e.g. Eqs.(3.12) and (3.13)). Hence, the simulation procedure will be simultaneously carried out in two groups and due to the transversely anisotropy, each of them requires two computations. More specifically,

• for (P) E bounded problem, temperature loading T E with boundary thermal gradient respectively along e 1 and e 3 :

T E = E i x i , E = E 0 e 1 or E = E 0 e 3 (3.41)
• for (P) p bounded problem, temperature loading T P with boundary polarization vector respectively along e 1 and e 3 :

T P = |E| ∂G ∂x i P i , P = P 0 e 1 or P = P 0 e 3 (3.42)

As the considered media being defined in a finite domain (i.e. porous media), the numerical average method is adopted for the estimation of the resistivity contribution tensor H E 0 . Following Eq.(3.35), the resistivity concentration tensors A E and A P will be firstly computed. The consistency laws ensure that the subjected macroscopic thermal gradients are also the average ones, one has:

⟨ε⟩ E D = E = (1 -f ) ⟨ε⟩ E D M + f ⟨ε⟩ E E ⟨ε⟩ P D = (1 -f ) ⟨ε⟩ P D M + f ⟨ε⟩ P E (3.43)
where ⟨ε⟩ E D and ⟨ε⟩ P D denote the subjected macroscopic thermal gradient respectively for the (P) E bounded and (P) P bounded 73 Chapter 3

3.2 problems, ⟨ε⟩ E E , ⟨ε⟩ P E , ⟨ε⟩ E D M and ⟨ε⟩ P D M are the corresponding average thermal gradient in the porous phase and those in the matrix one.

Note that in Eq.(3.43), the macroscopic thermal gradient ⟨ε⟩ E D of the (P) E bounded problem is known as predefined, while ⟨ε⟩ P D of the (P) P bounded problem should be calculated. Following the Gauss theorem, the latter one can be obtained as an integral over the external boundary ∂D :

⟨ε i ⟩ P D = 1 | D | M m=1 T P n i (m) S (m) (3.44)
where M is the total number of the surface elements, T P denotes the nodal temperature of integration points, S (m) is the area of the m-th one and n i defines the unit normal of the outer surfaces.

Since the pore is insulated, the heat flux in the porous phase vanishes (i.e. ⟨σ⟩ E E = ⟨σ⟩ P E = 0), it follows that:

⟨σ i ⟩ D M = ⟨σ i ⟩ D = 1 | D M | N n=1 (σ i ) (n) V n (3.45)
with N being the total number of the volume elements and V (n) giving the volume of the n-th one.

Following Fourier's Laws, ⟨ε i ⟩ E D M and ⟨ε i ⟩ P D M can then be calculated through:

⟨ε i ⟩ D M = -r i j • ⟨σ i ⟩ D M
Having in hand the above computed quantities, the average thermal gradient fields in the porous phase are expressed as:

⟨ε⟩ E E = E -(1 -f ) ⟨ε⟩ E D M f ⟨ε⟩ P E = ⟨ε⟩ P D -(1 -f ) ⟨ε⟩ P D M f (3.46)
The concentration and contribution tensors can then be computed by Eqs.(3.19) and (3.35). Due to the transversely isotropy, both of them have 3 non-zero components which are A 11 = A 22 , A 33 , H 11 = H 22 , H 33 in the present study.

Assessment and validation of the proposed numerical procedure

In this section, the proposed numerical procedure will be assessed and validated by comparing its predictions of the resistivity contribution tensor H E 0 with some analytical and numerical results in literature. To this end, a preliminary step aiming at defining an appropriate scale of the bounded RVE will be first carried out. It will be accomplished by studying the accuracy and the efficiency of the proposed numerical procedure in the case of a spherical pore planted in an isotropic matrix. Next, we will systematically consider the superspherical pores embedded in an isotropic matrix as well as the spheroidal ones surrounded by a transversely isotropic matrix to respectively justify its accuracy on the concavity of the pores and that on the anisotropy of the matrix material.

Spherical pore in isotropic matrix

Let us consider a spherical pore of radius a embedded in the center of a cubic RVE with side length 2L. By varying the so-called scale ratio a/L, the same mesh of each geometry with fixed a/L will be separately subjected to the corrected boundary conditions and the uncorrected ones (i.e. classical Hashin type boundary conditions). Due to the symmetry of the problem, only one independent component of resistivity contribution tensor H 11 is numerically computed that will be next compared to the corresponding analytical solution Eq.(A.81). The relative errors are displayed in Fig. 3.3. It can be observed that the computation convergence with the corrected boundary conditions occurs and stabilizes even when a/L ≃ 3 with an excellent accuracy that the relative errors are around the value of 10 -4 . Whereas it seems like that those obtained from the classical modeling just begin stabilize when the scale ratio a/L is much bigger (e.g. a/L ≥ 10) and present important relative errors with respect to the former ones.

This can be interpreted as, on one hand the corrections of the boundary conditions allow to accelerate the numerical convergence without degrading the computation accuracy, on the other hand when the scale ratio a/L is sufficiently big, the representative bounded model tends to be an infinite one such that the correction of the boundary conditions is hence less efficient or might be useless in the numerical modeling. By considering the computation results as illustrated in Fig. 3.3, we adopt in an ad-hoc manner the scale ratio a/L = 8 in the next part of this work. This is of course an approximation but will be shown as sufficiently accurate in the following numerical estimations. Chapter 3

3.2

Superspheroidal pore in isotropic matrix

Table .3.2 displays the numerical results of the resistivity contribution tensor for the superspheroidal pores surrounded by an isotropic matrix as well as their comparisons with the numerical results obtained by Sevostianov et al. (2016a). Different geometries and meshes are realized by varying the value of the concavity parameter p ∈ [0.2, 1]. We pay our first attention to the particular case with p = 1 (i.e. spherical pore). It can be observed that, by comparing with the computations of Sevostianov et al. (2016a), the FEM results obtained in this work present much smaller relative errors with respect to the analytical ones (i.e. H 11 = H 33 = 1.5). Moreover, the relative error of the FEM results with respect to those of Sevostianov et al. (2016a) is shown to be increased as we decrease the concavity parameter p toward the limit value of 0. This difference becomes much significant and cannot be ignored especially for the estimation of H 11 in the case of the concave pores (i.e. p < 0.5). In addition to the corrected boundary conditions that accelerate the computation convergence and preserve the accuracy of estimations, this might also be due to the fact that, as shown in Fig. 3.2, the meshes realized in this work are sufficiently refined in the transition zone between the matrix and porous phases that allows to reduce the effect of the geometric singularity as mush as possible 3 . In return, it reveals that the numerical estimations obtained in the present work from the proposed numerical procedure and the refined meshes can be considered as accurate and precise in the case of isotropic matrix. Relative error of FEM results with respect to the analytical one for p = 1.

Table 3.3 shows the numerical estimation of the resistivity contribution tensor for the superspherical pores 3 For more details on the mesh information, readers are referred to A.3 and the Table B.2 of Sevostianov et al. (2016a).

surrounded by an isotropic matrix as well as its comparison with respect to that obtained by [START_REF] Chen | Micromechanical Modelling of Heterogenous Porous Materials with Application to Oolitic Rocks[END_REF]. It is worthy to mention that, unlike the superspheroid pores that present the axisymmetry, the superspherical ones are symmetric in all of the three principal directions. Consequently, only H 11 is displayed in Table 3.3 since its three principal components are identical (i.e. H 11 = H 22 = H 33 ) and, according to the boundary conditions, the deviatoric ones are all null (i.e. H i j = 0 if i j). Once again, it can be observed that the FE computation realized in the present work is much more rigorous especially when the pore is concave (i.e. p < 0.5). In order to understand better the aforementioned conclusion, as the superspheroidal pore tends to be a penny crack crossed by a perpendicular needle along the symmetry axis x 3 when the concavity p → 0, we show in Table 3.4 the comparison between the FEM computations in a limiting case with p = 0.2, the associated results of Sevostianov et al. (2016a) and the analytical ones of an aligned penny crack embedded in the isotropic host matrix. The geometry of the latter one is approximated by considering a very small thickness e such that e/L = 0.001 (see Fig: 3.1d) and the corresponding analytical solution can be calculated from Eq.(A.82) with the anisotropy parameter ν = 1. It is convenient to note here that the obtained resistivity contribution tensor is normalized by the volume fraction V crack V sphere to avoid the volume effect.. According to the comparison, although a good agreement can be observed among all of the three numerical results and the analytical one, the FE computation obtained from the proposed numerical model presents a very small relative error than the one calculated by Sevostianov et al. (2016a).

Spherical pore in transversely isotropic matrix

This section deals with the assessment and validation of the numerical procedure by paying particular attention to the matrix anisotropy. The thermal problem is of lower order. This allows us to involve an additional parameter: matrix anisotropy of conductivity κ. In this light, we switch off the effect of the concavity parameter p by considering that the pore is in a spherical form and the matrix is transversely isotropic. Note once again that the proposed numerical procedure can also be applied to any type of anisotropy, which probably requires supplementary simulations and is not further pursued here. Table 3.5 shows the FE computation by varying the anisotropy parameter of the matrix , 10] as well as its comparison with repsect to the analytical results (see for instance Section 3.2.2.2 and A.9). A very good agreement can be found and the relative errors are shown to be minor. it justifies that the accuracy of the proposed numerical method is unaffected by the matrix anisotropy. In this section, we propose to carry out the study of the resistivity contribution tensor H E 0 in the case of transversely isotropic matrix comprising concave/convex pore to understand in more detail their combined effect. This will be carried out by separately varying the anisotropy parameter of the matrix, denoted by κ, and the concavity one p of the concave/convex pores (i.e. superspheroidal and superspherical pores) that will both be defined in relatively large intervals such as κ ∈ [0.1, 10] and p ∈ [0.2, 5]. Again, we restrict the study, particularly in the superspheroidal cases which verify the symmetry of revolution, to the assumption that the directions of the symmetry between the matrix anisotropy and that of the pore are aligned on the same direction e 3 (see for instance Figs.3.1a and 3.1b).

κ = λ 1 0 /λ 3 0 in the interval κ ∈ [0.1
In return, for the FE computations, two non-zero and independent components of H E 0 (i.e. H 11 and H 33 ) will be computed. This is achieved by the proposed numerical procedure that was described in Section 3.2.3 by fixing the scale ratio of the RVE a/L = 8 and only subjected to the corrected boundary conditions (with remark "with correction" in numerical results). For the sake of prediction accuracy, the numerical computation will be carried out based on the sufficiently refined meshes, for which the numbers of nodes and elements are detailed in A.3. Specifically, Section 3.2.5.1 investigates firstly the predictions in the case of superspherical pores. As expected, the combined effect of the matrix anisotropy and that of the pore concavity is quantitatively obtained. This effect on the response of the resistivity contribution tensor is also clearly illustrated in the p -κ -H i j spaces. In the remaining Section 3.2.5.2, similar study is realized for the superspheroidal pore embedded in the transversely isotropic matrix.

Superspherical pore in transversely isotropic matrix

Tables 3.6 and 3.7 respectively summarize the numerical estimations of H 11 and H 33 for the superspherical pores.

First of all, it can be found that, by fixing the concavity parameter p, their predictions decrease with the increase of the anisotropy parameter κ. A second very interesting result, observed in any case of a fixed value κ, is that the influence of the concavity parameter p is not monotone. More specifically, the predictions of H 11 and those of H 33 both increase with the decrease of the concavity when p is small, in which this dependence becomes more obvious especially when the pore is concave (i.e. p < 0.5). While this effect becomes opposite when the pore is severely convex (i.e. p → 5). Similar observations were done in the numerical study of [START_REF] Chen | Micromechanical Modelling of Heterogenous Porous Materials with Application to Oolitic Rocks[END_REF], albeit for an isotropic matrix. The turning points of monotonicity are marked with yellow box in the tables and not all of them fall on the spherical shape (i.e.p = 1), which indicate a certain correlation with the anisotropy degree of matrix κ. Moreover, we present in Tables 3.6 and 3.7 the comparison between the numerical estimations of spherical pores p = 1 with the corresponding theoretical results. An excellent agreement can be found with the relative error being around 10 -5 that could indirectly justify the accuracy of the proposed numerical model for convex shapes. R.E. 1 0.013% 0.005% 0.008% 0.008% 0.008% 0.008% 0.008% 0.008% 0.008% 1 Relative error of numerical H 33 with respect to analytical one for p = 1

For completeness, Figs.3.4 and 3.5 illustrate clearly the combined effect of the matrix anisotropy κ and the concavity parameter p. Each of their influences is also projected in the planes κ -H i j and p -H i j , respectively. As discussed before, these effects are found to be much significant when the pore is concave with an important anisotropy of the matrix4 . Specifically, in Fig. 3.4, we observe that the prediction of H 11 presents a more important evolution and a slower subsequent saturation in the plane κ -H 11 than that in the plane p -H 11 , which quantitatively and qualitatively Chapter 3

3.2

Concluding remarks

In the present work, we have numerically evaluated the resistivity contribution tensor of the concave pore inhomogeneity embedded in a transversely isotropic matrix. This has been realized by use of an original developed nu-merical homogenization method complying with the adapted boundary conditions based method recently formulated by [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF] (see also [START_REF] Du | Numerical computation of compliance contribution tensor of a concave pore embedded in a transversely isotropic matrix[END_REF]). The proposed numerical procedure was carried out for a bounded representative elementary volume and is shown to be efficient and accurate in the numerical modeling.

By paying particular attentions to the pore concavity and the matrix anisotropy, a major contribution of this work is found as the sufficiently exact computation results and analysis that illustrate their significant effect on the thermal conductivity properties. Specifically, the adapted boundary conditions based method was extended in the context of the matrix anisotropy thanks to the Green function and its gradients applied in the correction of the boundary conditions. The numerical homogenization method has been utilized in the proposed numerical procedure that has been firstly assessed and validated by comparing its predictions with the analytical and existing analytical and numerical results in particular cases. It is then used to the investigate the effect of the pore concavity on different transversely isotropic matrix from the quantitative estimates of the resistivity contribution tensor, which was found to be of critical importance especially in the case of concave pores.

Finally, the proposed numerical method is able to deal with any general anisotropy of the matrix material but such a study has not been attempted here for the sake of keeping the work focused and concise. In the perspective point of view, it is clear from the previous simulations that the effective properties such as those predicted from semi-analytical homogenization models could be further developed.

Effective thermal properties of a composite containing transversely isotropic matrix and irregular shapes of pores

The effective thermal properties of an anisotropic solid containing randomly or aligned oriented pores depend on the conductivity of the minerals, on porosity as well as pore shape. In the text to follow, the previous numerical results in section 3.2, where resistivity contribution tensors of concave pores in a transversely-isotropic material are obtained, are firstly approximated to analytical expressions. Hereafter, the effect of the concavity parameters and the anisotropic degree of matrix on the overall thermal properties of a porous material with transversely-isotropic solid phase are discussed. For this goal, various micromechanical schemes: Non-interaction approximation, Mori-Tanaka-Benveniste and Maxwell schemes are carried out to predict the effective thermal properties, using the approximated formulas of resistivity contribution tensors.

Resistivity contribution tensor of a concave pore

The prediction of the macroscopic mechanical response of porous materials is established based on semi-analytical approximations via mean-field homogenization schemes. One of the approaches to characterize contribution of microstructure (irregularly shaped pores, anisotropic matrix and son on) to the effective thermal properties of heterogeous materials is based on the evaluation of their conductivity and resistivity contribution tensors. The property contribution tensors of linear thermal problem are recalled in 3.1. Resistivity contribution tensor H of an individual pore of superspherical and axisymmetric superspheroidal shapes embedded in a transversely-isotropic material was computed numerically via finite element method (FEM) in a companion paper. For linear thermal problem, the resistivity contribution tensors H E 0 have two independent thermal coefficients: H E 11 and H E 33 for porous materials containing transeversly isotropic matrix. We present in next section a method that obtain the explicit approximate expressions of the two contribution coefficients by considering a variable concavity parameter and anisotropic degree of matrix.

The two reference shapes are mathematically described as:

• superspheroidal pore |x 1 | 2p + |x 2 | 2p + | x 3 ς | 2p = a 2p (3.47)
• axisymmetrical superspheroidal pore

x 2 1 + x 2 2 p + | x 3 ς | 2p = a 2p (3.48)
p, ς and a respectively denote non dimensional concavity parameter, aspect ratio and semi-lengths in plane 0x 1 x 2 .

Theses shapes are convex in the range p > 0.5 and concave for p < 0.5. Both shapes degenerate into a spheroid when

p = 1.
In what follows we will only consider ς = 1 and a = 1, the first shape is then a supersphere of unit semi lengths on x i axis, and the second shape obtained by a rotation about symmetry axis x 3 . Supersphere and axisymmetrical superspheroid coincide with sphere in the case p = 1 but strongly differ in the limiting case p → 0.

3.3.2. Approximation formula for resistivity contribution tensor of a superspherical or axisymmetrical superspheroidal pore embedded in a transversely isotropic host matrix

In this investigation, thermal analytical approximation formulas for the resistivity contribution tensor of porous materials, containing 3D pores of particular shapes and transversely isotropic matrix are presented. Two type of microstructures are considered: supersphere and axisymmetrical superspheroid. We restrict the study to the following assumptions:

• same directions of symmetry between matrix and pore inclusion (aligned case)

• study is focused on the concavity parameter p and anisotropic degree of the host matrix κ =

λ 0 1 λ 0 3 with λ 0 3 = 1W • K -1 • m -1
Using these expressions would considerably simplify the process of calculating the property contribution tensor to obtain the distributions of pore shape and the host matrix.

Volume and surface area of superspherical and axisymmetrical superspheroidal pores

Approximation formula may be obtained by using basic geometric information related to the considered reference shapes, supersphere and axisymmetrical superspheroid, defined in relations (3.47-3.48), with aspect ratio ς = 1 and a = 1. These informations are volume, total surface area and corresponding volumes write (Γ denotes Euler Gamma domains. Maximal relative errors of approximate results compared to FEM results of superspherical pore are given in

f
H so ii (κ, p) = A so (p)/(V so (p)) 2/3
A sphere / V sphere 2/3 f so ii (κ, p) H sphere ii

, no sum over i, i ∈ [1, 3] (3.54)

where sphere corresponds to the spherical case (particular case when p = 1).

Functions f so ii (κ, p) are the same formula as Eq.(3.53). The polynomial coefficients of f so 11 and f so 33 are presented in Table .3.12 with different piecewise ranges of κ and p. Comparisons between the results of approximate relation Eq. (3.54) and FE results (in Table.3.8 & 3.9) are presented in Fig. 3.9. The corresponding maximal relative errors of approximate results compared to FEM results for axisymmetrical superspheroidal pore are given in Table .3.11 which are lower than 1.5%.

So far, based on the FE results obtained in section 3.2 and analytical solution for sphere (when p = 1), we built the approximate analytical formulas to predict thermal contribution from all possible superspherical and axisymmetrical superspheroidal pores embedded in all possible transversely isotropic matrix in the parameters frame, with high accuracy. Figure 3.9: The two independent components H so 11 and HH so 33 of a superspheroidal pore embedded in TI matrix, as a function of concavity parameter p and anisotropic degree of host matrix κ. Comparison between FEM results (star points) and approximate relations (plain surface).

Evaluation of the effective thermal properties of materials with transversely isotropic matrix

The approximate expressions for the pore property contribution tensors H derived in the section 3.3.2 have been used as input in the different homogenization schemes such as Non Interaction Approximation, Mori-Tanaka-Benveniste and Maxwell schemes (respectively referred with superscripts NIA, MTB and MX), see [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF]; [START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF]; [START_REF] Kachanov | Micromechanics of materials, with applications[END_REF]), to obtain approximate analytical solution for the effective thermal properties. The equations (3.55) yielded by the approximate expressions (3.52) & (3.54) hold for a wide range of concavity parameter p ∈ [0.2, 1], anisotropic degree of matrix κ ∈ [0.1, 10]) and porosity φ.

r NIA = r 0 + φH E 0 , r MTB = r 0 + φ 1 -φ H E 0 , r MX = r 0 + 1 φ H E 0 -1 -Q Ω 0 -1 (3.55)
where φ denotes the porosity. Q Ω 0 denotes the Hill tensor of the effective inclusion of the Maxwell scheme, which is supposed of spheroidal shape (with aspect ratio γ e ) and aligned with the directions of the TI host matrix. Q Ω 0 is related to the Hill tensor P Ω 0 (see A.9) by the relation:

Q Ω 0 = r 0 • I -P Ω 0 • r 0 (3.56)
where r 0 is conductivity tensor of matrix.

In this section, the sensitivity analysis of effective thermal conductivity with respect to the shape of pore p and contrast in matrix properties κ is discussed. The porosity of clay matrix φ (denoted f I p in [START_REF] Giraud | Effective porothermoelastic properties of transversely isotropic rock-like composites[END_REF]) is approximately comprised in the range φ ≤ 0.30 for the sensitivity analysis.

Effective thermal conductivity of superspherical pores with MTB scheme

This section attempts to quantitatively investigate the performance of the solutions obtained from MTB approximation in the estimation of effective thermal properties. Effective thermal coefficients λ MTB 1 and λ MTB 3 obtained with MTB approximation of a random orientation distribution of superspherical pores in the transverse plane embedded in transversely isotropic matrix are presented in figures 3.10. Approximation formula Eq.3.52 deduced from FEM are compared to approximation of resistivity contribution tensor of an aligned oblate spheroidal pore of same volume (the semi axis length of supersphere is equal to the maximum axis radius of the oblate spheroid) Eq.(3.57):

H i j (p) ≈ H spheroid i j (γ(p)), 0.2 ≤ p ≤ 1 (3.

It may be observed that

• the contribution of these two shapes to thermal effective properties are pronouncedly different in the different concave range and contrast degree range of matrix. It is expected that the supersphere tends to three orthogonal needles (with zero volume and zero surface) when p tends to zero, whereas the axisymmetrical superspheroid tends to a circular crack with one central orthogonal needle (the latter has zero volume but non zero surface).

Note that, the most relevant shape compared to microstructure of porous materials is certainly supersphere.

• obtained effective porous material is transversely isotropic with same symmetry axis than matrix. The two shapes bring the opposite effects in different directions of thermal coefficients. In the transverse plane (for λ 1 ), superspherical pore has more obvious affect than axisymmetrical superspheroidal pore, while inverse for the conductivity along symmetry axis (λ 3 ).

Comparisons between NIA, MTB and Maxwell homogenization schemes

Effective thermal properties predicted by NIA, MTB and Maxwell are presented in Fig. 3.13 respectively for superspherical pores randomly oriented in transverse plane and aligned axisymmetrical superspheroidal. Predictions of effective conductivity are close at low porosity φ and lay gradually away at higher porosity. The Maxwell estimate with effective spheroid inclusion of aspect ratio γ e = 0.5 follows the MTB results and drops increasingly beneath the NIA at higher porosity. The shape of the effective inclusion of the Maxwell scheme is still an open issue when host matrix is anisotropic (see [START_REF] Sevostianov | On the shape of effective inclusion in the Maxwell homogenization scheme for anisotropic elastic composites[END_REF]; [START_REF] Giraud | Effective electrical conductivity of transversely isotropic rocks with arbitrarily oriented ellipsoidal inclusions[END_REF]; [START_REF] Sevostianov | Maxwell's methodology of estimating effective properties: Alive and well[END_REF]).

Introduction

In the past decades, the mechanics of composite porous solids was developed to predict the effective properties both in theoretical and computational approaches. In the context of theoretical modelling, some pioneering contributions (Hashin andShtrikman (1962a, 1963); [START_REF] Hill | Elastic properties of reinforced solids: Some theoretical principles[END_REF]) was realized, by adopting variational approach, for the prediction of the bounds for the effective elastic moduli of multi-phases materials with an isotropic matrix. Due to some mathematical difficulties, the method of pure analytical homogenization is limited especially for anisotropic medium and complex 3D microstructures. Then these problems generally require computational approaches.

Simplified "model composites" were studied by adopting appropriate numerical engineering methods, for which the microstructure is considered as periodic described via a Representative Volume Element (RVE). For instance, [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF] considered two different families of numerical methods to deal with the composites with periodic microstructure respectively based on the Finite Element Method and that of Fast Fourier Transforms. They gave their comparison results for spherical impenetrable particles placed randomly in a cubic unit cell with volume fraction up to 26.78%. [START_REF] Segurado | A numerical approximation to the elastic properties of sphere-reinforced composites[END_REF] then developed a modified random sequential adsorption algorithm (originally proposed by [START_REF] Widom | Random Sequential Addition of Hard Spheres to a Volume[END_REF]) to generate 3D cubic unit cells containing non-overlapping identical spheres, embedded in a continuous and isotropic elastic matrix, in which different materials with specific microstructures were studied, such as those containing rigid spherical inclusions and voids in an elastic matrix as well as typical composite made up of glass spheres in an epoxy resin. [START_REF] Kachanov | Handbook of Elasticity Solutions[END_REF] studied the contribution tensors of an inhomogeneity embedded in an isotropic solid that are of direct relevance to the effective elastic properties without the consideration of the Eshelby theory. In the context of Effective Media Theory (EMT), this contribution tensor approach was extended to the case of non ellipsoidal inhomogeneities for the estimation of the effective elastic properties and thermal conductivity. [START_REF] Bohm | Multi-Inclusion Unit Cell Studies of Reinforcement Stresses and Particle Failure in Discontinuously Reinforced Ductile Matrix Composites[END_REF] developed three-dimensional unit cell models that contain randomly positioned spherical inclusions and/or fibers for investigating the thermomechanical behavior of ductile matrix composites. [START_REF] Eroshkin | On micromechanical modeling of particulate composites with inclusions of various shapes[END_REF] proposed a methodology to analytically or numerically predict the effective moduli of the two-phases composites by taking into account the compliance contribution tensor H of single inclusion problem. [START_REF] David | Compressibility and shear compliance of spheroidal pores: Exact derivation via the Eshelby tensor, and asymptotic expressions in limiting cases[END_REF] derived the explicit expressions of the elastic compliance of a spheroidal pore in an isotropic solid via Eshelby's theory that lead some asymptotic expressions with respect to the pore compressibility as well as the shear compliance (see [START_REF] Wu | The effect of inclusion shape on the elastic moduli of a two-phase material[END_REF]) of "crack-like" and "needle-like" pores in a finite aspect ratio. [START_REF] Drach | Comparison of full field and single pore approaches to homogenization of linearly elastic materials with pores of regular and irregular shapes[END_REF]

E 11 1.0e -5 0 0 0 0 0 E 22 0 1.0e -5 0 0 0 0 E 33 0 0 1.0e -5 0 0 0 E 12 0 0 0 1.0e -5 0 0 E 13 0 0 0 0 1.0e -5 0 E 23 0 0 0 0 0 1.0e -5
script to calculate effective elastic properties through Eq.(4.7):

C e f f i jkl = σ i j m ε 0 kl m (4.7)
where σ i j m are the macroscopic stress field (i.e. average volume stress field) and ε 0 kl m denotes the macroscopic strains in the m th load case (i.e. applied strain components). Effective engineering constants can then be obtained from the effective compliance tensor.

Assessment of the numerical algorithm

Note here that matrix Poisson's ratio ν 0 = 0.3 was assumed in all isotropic cases unless otherwise specified.

Convergence in terms of realizations

The convergence analysis in terms of the number of realizations is conducted for materials containing isotropic matrix and randomly distributed spherical pores with porosity φ = 0.25 and aspect ratio a = b, c/a = 1.0, at fixed number of pores (n = 100) and discretization, as shown in 4.3. Then the convergence analysis in terms of the number of realizations is conducted for materials containing isotropic matrix and randomly distributed spheroidal pores with porosity φ = 0.25 and aspect ratio a = b, c/a = 2.0, at fixed number of pores (n = 100) and discretization, as illustated in 4.4.

Isotropic matrix with randomly spherical pores

In the context of the assessment of the proposed numerical homogenization procedure, the accuracy of some homogenization models is then studied for the composites made of an isotropic matrix embedded by isotropic spherical pores. The convergence of the RVE is studied as similar as that described in Section 4.2 with respect to the number of spherical pores of the microstructure (i.e. RVE), which is defined as the number of pores from which its increase does not clearly affects the macroscopic elastic character. Normalized effective Young's modulus , 5, 10, 20, 50, 100, 200, 300]) are estimated via direct FEA simulations in the case of the 4.3 porosity being equal to φ = 0.25. As shown in Fig. 4.5, the results were compared to the predictions of single pore solutions by referring MTB and Maxwell homogenization schemes for which H-tensors analytically computed in Section 2 (see also [START_REF] Du | Numerical computation of compliance contribution tensor of a concave pore embedded in a transversely isotropic matrix[END_REF][START_REF] Du | Effective elastic properties of transversely isotropic materials with concave pores[END_REF]) are taken into account. Note that some classical schemes, such as the Hashin Shtrikman Bounds (HSB) (Hashin andShtrikman (1962b, 1963) converges when the number of pores exceeds 100. The composite material containing randomly spherical pores remains quasi-isotropy when comparing the effective engineering constants in different principal axis, which is logical because that a completely random distribution of spherical pores theoretically leads to an isotropic effective response.

In practice, for each realization, C e f f is not strictly isotropic due to the finite number of pores during the numerical modeling. Consequently, the effective bulk and shear modulus for each realization are obtained by the isotropic projection C as follows:

C iso = 3 KJ + 2 GK with K = 1 3 C :: J = 1 9 C ii j j and G = 1 10 C :: K = 1 10 C i ji j -9 K (4.8)
Here, K and G are the isotropized effective bulk and shear moduli. The isotropic projectors are defined by

J = 1 3 i ⊗ i and K = I -J (4.9)
with J, K and I being the fourth-order hydrostatic, deviatoric projection tensors and identity tensor respectively.

In the light of convergence analysis, the number of pores 100 is chosen in the following calculations. Direct FEA simulations with porosity φ = 0.05, 0.10, 0.15, 0.20, 0.25 are applied to the microstructures with spherical pores in the case of monodispersion, as presented in Fig. 4.6. The direct FEA results of effective bulk and shear modulus are compared to those obtained from several analytical results [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF]; Kachanov and Sevostianov

Isotropic matrix with randomly spheroidal pores

We study in this section the convergence of the proposed numerical procedure in the case of porous media with an isotropic matrix and randomly distributed spheroidal pores with the porosity φ = 0.25 and aspect ratio a = b, c/a = 2.0. here that in the present study, we can observe that the totally random dispersion of the pores still brings the some transverse isotropy (higher stiffness in the normal direction E 3 > E 1 and quasi-same in the other two directions

E 1 = E 2 )
while an isotropic response should theoretically be homogenized. The analysis of deviation from isotropy is presented in A.12. This implies that this isotropic overall response is very difficult to obtain by the adopted RSA method even with a large number of pores due to the finite number of the pores (maximal value 300 in our study), whereas the difference in different axes is not very important.

Unlike the spherical case, the macroscopic effective elastic properties depends on the dispersion of the spheroidal pores. The orientation and the distribution of the pores should be therefore considered in the analytical and/or semi-analytical homogenization procedure especially for the complex microstructures with some manner of pores the anisotropy of matrix presents an important influence on the macroscopic effective response, which qualitatively follows the same one of the matrix material (higher stiffness in the normal direction E 3 than that of the E 1 = E 2 ones).

It also illustrates that, in the simulation point of view, a number of pores n ∈ [50, 100] is found to be sufficient for the related microstructure representation. In order to avoid repeating mesh generation, we chose n = 100 in the following part of this subsection. a good agreement with respect to the analytical/semi-analytical results specifically at relatively lower porosities. The effective inclusion is chosen as spherical shape for Maxwell scheme due to the random distribution of spherical pores in modeling, where its results mathematically coincide with those of MTB. Through comparative studies, it also shows that MTB and Maxwell homogenization models stand better than NIA as being more accurate in the porosity range of less than 0.25, but still lost some accuracy of Poisson ratio and shear moduli at higher volume fractions.

Transversely isotropic matrix with randomly spheroidal pores

This section deals with the numerical estimation of the macroscopic effective response of the porous material with a transversely isotropic matrix embedded by a finite number of randomly distributed and orientated spheroidal pores.

The transversely isotropic elasticity of the matrix is described by the parameters of Table 4.3. The aspect ratio is taken as γ = c/a = 2. Fig. 4.12 shows the effective Young and shear modulii as well as the Poisson ratios as dependence of the porosity as well as the comparison of their estimations with those in the case of spherical pores (i.e. γ = 1 in section 4.4.1). It can be first observed that the effective elastic moduli decrease with the increase of the porosity, except for the Poisson ratio ν 12 . The comparisons also qualitatively reveals that the different morphology of the pores (i.e. only for spherical form and/or spheroidal one considered here) may have no significant affect on the effective elastic properties of porous materials in the case of the randomly distributed and oriented pores which theoretically leading to the transversely isotropic response as the one of the matrix material. It is also interesting to note that following Fig. 4.12, this fact is more obvious at small porosities. fields. It can be also observed that the stress concentration generally occurs at the interface between matrix and pore and the interaction area of the nearby pores. The stress field near the interface perpendicular to loading path is usually higher than other domain. The maximum Von Mises stress under the load in the X/Y direction is higher than that in the Z direction, which is also consistent with the characteristics of the material.

Conclusion

In this section, the effective elastic properties of the porous material were studied by adopting the numerical homogenization method that was applied to a RVE containing a finite number of randomly distributed and oriented pores. The RVE is represented by a cubic unit cell which is numerically generated by use of a Random Sequential Adsorption (RSA) algorithm. The matrix was respectively considered as isotropic and transversely isotropic embedded ether by spherical pores or the spheroidal ones. The numerical computation was compared to the analytical predictions obtained by some available micromechanical models, such as Non Interaction Approximation (NIA),

Mori-Tanaka-Benveniste (MTB), Maxwell schemes as well as some available numerical results in the literature, except for the one of transversely isotropic matrix with spheroidal pores. It was found that the pore morphology and the matrix anisotropy significantly affect the effective properties of the porous material. The MTB and Maxwell schemes does not introduce a significant influence on the corresponding estimations in the case of spherical pores planted in either isotropic or transversely isotropic matrix, and the one of spheroidal pores in isotropic matrix. Finally, this numerical procedure has not been applied to the complex microstructure with the finite number of concave pores due to their so strongly geometrical concavity that leads to an irresolvable difficulty for the meshing of the RVE. Such a study is now underway and would be realized and reported elsewhere.

analytical homogenization for their assessment and validation. More specifically, at small pore concentrations, the Non-Interaction Approximation provides sufficiently good accuracy. However, when the porosity increases, other more "advanced" micromechanical schemes, such as MTB and Maxwell, become necessary, which reflects that the influence of interaction between pores cannot be neglected.

Note again that all of the numerical methods proposed in this thesis are able to deal with any general anisotropy of a multiphased material but such a study has not been attempted here for the sake of keeping the work focused and concise. Application to other material systems (e.g. composites, fibrous materials, etc.) within elastic domain will bring a huge expansion. Continuous research could focus on the numerical homogenization modelling with composite porous geomaterials containing irregular forms (such as supersphere, superspheroid, penny crack, etc.) of defects embedded in an anisotropic matrix. It could be worth mentioning that the thermal conductivity problem is somehow equivalent to the magnetic and electric problem. With the current development of electro-and magnetoactive polymeric materials via 3D printing such inclusion shapes might be interesting from the point of view of enhanced magnetization for instance. Direct extension to coupling problems in material science such as thermoelastic, magnetic-elastic and piezo-electric composites could also be envisaged in the long term in a straightforward manner. It is also relevant to mention in this connection that due to the meshing difficulties encountered in this thesis, some meshfree methods, such as the one of Fast Fourier Transform (FFT) could be adopted. Last but not least, the consideration of the interfaces between matrix and pores, different phenomena (e.g. initiation and propagation of cracks, plasticity, etc) may also constitute a challenging extension.

unit vector of symmetry axis of the material, let us introduce the second-order tensors

i N = n ⊗ n = n i n j e i ⊗ e j , i T = i -i N (A.5)
In the particular case of n = e 3 , (A.5) writes i N = e 3 ⊗ e 3 , i T = e 1 ⊗ e 1 + e 2 ⊗ e 2 (A.6)

One introduces fourth-order tensors

E 1 = i N ⊗ i N , E 2 = 1 2 i T ⊗ i T , E 3 = 1 √ 2 i N ⊗ i T , E 4 = 1 √ 2 i T ⊗ i N (A.7) E 5 = i T ⊗i T - 1 2 i T ⊗ i T , E 6 = i T ⊗i N + i N ⊗i T (A.8)
It may be shown that any transversely isotropic fourth-order tensor can be decomposed as

L = 6 i=1 l i E i = l i E i (A.9)
Considering symmetry axis equal to n = e 3 , Walpole matrix representation of tensor H = h i E i writes We detail herefafter calculation of compliance contribution H E 0 of a spheroidal pore E in TI basis, in terms of Observing now that δ ′ (λ z) = δ ′ (z)/λ 2 for all λ > 0, using the definition of the derivative in the sense of distributions [START_REF] Schwartz | Théorie des distributions[END_REF][START_REF] Gel'fand | Generalized Functions[END_REF]) and still adopting the spherical coordinates associated to the pole x 0 yield the following expression grad G(x) = -1 8π 2 ∥x∥ 2 .30) It may be noticed that A.30 corresponds to relation 5.58 given in [START_REF] Mura | Micromechanics of defects in solids[END_REF], p. 33. Fourier transform based derivation and change of variable on the unit sphere used in this paper are similar to the ones detailed in [START_REF] Mura | Micromechanics of defects in solids[END_REF].

H =                         h 2 +h 5 2 h 2 -h 5 2 h 4 √ 2 0 0 0 h 2 -h 5 2 h 2 +h 5 2 h 4 √ 2 0 0 0 h 3 √ 2 h 3 √ 2 h 1 0 0 0 0 0 0 h 6 0 0 0 0 0 0 h 6 0 0 0 0 0 0 h 5                         (A.
= -K -1 • ∂K ∂z • K -1 ) grad G(x) = 1 8π 2 ∥x∥ 2 2π φ=0 K -1 (u φ ) • u φ • C • e + e • C • u φ • K -1 (u φ ) ⊗ u φ -K -1 (u φ ) ⊗ e dφ (A
A.2.2. Exact 3D elastic Green function in the transversely isotropic case

We briefly recall exact Green function given in [START_REF] Pouya | Green's function solution and displacement potentials for transformed transversely isotropic materials[END_REF]. We only consider non degenerate case c-c 3 √ 2 -c 6 0 with n = e 3 (n denotes unit vector on symmetry axis of transversely isotropic material). See [START_REF] Pouya | Green's function solution and displacement potentials for transformed transversely isotropic materials[END_REF] for discussion on more general and non degenerate cases. We recall that Green function is not singular on the axis x 3 , except at the origin x 3 = 0. In the case ρ → 0 and

x 3 0, it writes G 0 i j x 3 e 3 = 1 2 π c 6 ∥x∥ g (ρ=0) 11 δ 1i δ 1 j + δ 2i δ 2 j + g (ρ=0) 33 δ 3i δ 3 j (A.44) We only recall in this section the solution of strain Hill polarization tensor of a spheroidal inclusion aligned in a transversely isotropic host matrix, in the non degenerate case (c -c 3 √ 2c 6 0) with n = e 3 . We use constants c, υ i defined in section (A.2.2) and we present Withers solution [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] as in [START_REF] Parnell | The Eshelby, Hill, Moment and Concentration Tensors for Ellipsoidal Inhomogeneities in the Newtonian Potential Problem and Linear Elastostatics[END_REF]. See also [START_REF] Laws | A note on penny-shaped cracks in transversely isotropic materials[END_REF]; [START_REF] Sevostianov | Effective elastic properties of matrix composites with transversely-isotropic phases[END_REF] for exact solutions derived from Fourier transfom integral. The function S(x) (with 3 ) characterizes the influence of the shape of the spheroid)

g (ρ=0) 11 = 1 2 A 1 υ 1 + A 2 υ 2 + 1 υ 2 3 , g (
S(x) = 1 1 -x 2 - x 1 -x 2          1 √ 1-x 2 arccos (x) , x < 1 1 √ x 2 -1 arcosh (x) ,
x > 1 Components p i of strain Hill tensor P E 0 in TI tensor basis write

p 1 = -2 2 i=1 L 2 i M i υ 5 i I 3 (υ i ), p 2 = 2 2 i=1 M i υ i I 1 (υ i ), p 3 = √ 2 2 i=1
L i M i υ 3 i I 3 (υ i ), p = p 3 (A.49)

p 5 = 2 i=1 M i υ i I 1 (υ i ) + D I 1 (υ 3 ) 2 , p 6 = 1 2 2 i=1
(1 + L i )M i υ 3 i (I 3 (υ i ) -2L i I 1 (υ i )) + Dυ I 3 (υ 3 ) 4 (A.50) In Figs.A.4 -A.9, we show some supplementary results of the compliance contribution tensor H for the ellipsoidal pores respectively with the aspect ratio γ = 1, 1/2 and 1/10. By simultaneously considering those of γ = 1/5 (see Section 2.3.5.1), we focus on the value of a/L at the starting point of the convergence between the related evolution respectively obtained from the proposed numerical procedure and the classical one. Finally, we fix a/L = 8 for the numerical computations in the cases of the concave pores. .9275 -0.2891 -0.2936 15.9062 47.1787 2.2130 5.3414 0.25 2.0047 -0.2903 -0.1911 6.3536 20.3159 2.2940 2.6196 0.3 2.1128 -0.2966 -0.1204 3.3695 11.7872 A.6. Background on property contribution tensors value of the logarithm is chosen such that the imaginary part has the smallest value and belongs to ] -π 2 , π 2 ] with the same branch cut as the square root.

I k (η) = 1 z=-1 z k z 2 + η 2 (1 -z 2 )
dz and J k (η 1 , η 2 ) = J 0 (η 1 , η 2 ) J 2 (η 1 , η 2 ) J 4 (η 1 , η 2 )

if η 1 = η 2 = 1 2 2 3 2 5 if η 1 = η 2 1 arcosh η 1 +η 1 √ η 2 1 -1 η 3 1 √ η 2 1 -1 η 1 √ η 2 1 -1-arcosh η 1 η 1 (η 2 1 -1) 3 2 (2+η 2 1 )
√ η 2 1 -1-3 η 1 arcosh η 1 (η 2 1 -1)

5 2
if η 1 η 2 (η 2 1 -1) I 0 (η 1 )-(η 2 2 -1) I 0 (η 2 ) η 2 2 -η 2 1 (η 2 1 -1) I 2 (η 1 )-(η 2 2 -1) I 2 (η 2 ) η 2 2 -η 2 1 (η 2 1 -1) I 4 (η 1 )-(η 2 2 -1) I 4 (η 2 ) η 2 2 -η 2 1 Table A.9: Calculation of integrals J k (η)

Related compliance contribution of the spheroidal pore is then deduced from relations (A.12-A.16).

The exact solution of compliance contribution tensor of the penny shaped crack, H = lim ω→0 ω (Q 0 ) -1 , is detailed in [START_REF] Barthélémy | Simplified approach to the derivation of the relationship between hill polarization tensors of transformed problems and applications[END_REF]. We only recall final results. One uses

                                     p 1 1 = π 2 (γ 1 + γ 2 ) C 1111 -(γ 2 1 + γ 1 γ 2 + γ 2 2 ) C 2323 C 2323 C 3333 p 1 2 = π 4 (γ 1 + γ 2 ) 1 γ 1 γ 2 C 3333 + 1 C 2323 p 1 3 = - √ 2 π 4 (γ 1 + γ 2 ) C 1133 + C 2323 C 2323 C 3333 p 1 5 = π 8 1 γ 1 γ 2 (γ 1 + γ 2 ) C 3333 + 1 γ 1 + γ 2 + 1 γ 3 1 C 2323 p 1 6 = π 8 C 1111 -2 γ 1 γ 2 C 1133 -γ 1 γ 2 (γ 2 1 + γ 1 γ 2 + γ 2 2 ) C 3333 γ 1 γ 2 (γ 1 + γ 2 ) C 2323 C 3333 - γ 3 C 2323 (A.65a) (A.65b) (A.65c) (A.65d) (A.65e)
γ 1 , γ 2 , γ 3 have been previously defined. Compliance contribution tensor of a penny shaped crack aligned in a transversely isotropic host matrix writes .66) where p 1 1 , p 1 2 , p 1 3 and p 1 6 are given in (A.65a), (A.65b), (A.65c) and (A.65e). One deduces the components denoted h c 1 , h c 6 and used in approximate relations of compliance contribution tensor of an aligned axisymmetric superspheroidal pore: The expanded formula for gradG i 0 component writes: .74) where λ 1 0 , λ 2 0 and λ 3 0 are the three conductivity parameters of anisotropic matrix.

H = lim ω→0 ω (Q 0 ) -1 = - 1 C 2 3333 p 1 1 + 2 C 2 1133 p 1 2 + 2 √ 2 C 3333 C 1133 p 1 3 E 1 - 1 4 C 2 2323 p 1 6 E 6 (A
h c 1 = - 1 C 2 3333 p 1 1 + 2 C 2
gradG i 0 = 1 4 π λ 1 0 λ 2 0 λ 3 0 x i λ i 0 1 x 2 1 λ 1 0 + x 2 2 λ 2 0 + x 2 3 λ 3 0 3/2 (A
A.11. Complementary assessment of the numerical results for full field modelling

The convergence analysis in terms of the number of realizations is conducted for transversely isotropic materials containing isotropic matrix and randomly distributed spheroidal pores with porosity φ = 0.25 and aspect ratio a = b, c/a = 2.0, at fixed number of pores (100) and discretization, as illustated in A.12. A.12. Deviation from isotropy

The deviation from isotropy δ iso of the effective stiffness tensor C is evaluated with the normalized common Euclidean distance: .83) where Frobenius norm of the tensor C: ∥C∥ F = tr C : C T Table A.10 presents the values of the deviation for different concerned shapes of pores embedded in isotropic matrix. The spherical pores lead to the quasi-isotropic effective response with δ iso < 0.32% for monodisperse distributions. As one increases the aspect ratios, the deviation from isotropy increases. Even for a high number of pore n = 300 in our work, the transversely isotropic shape of the pore affects the evaluation of the elastic properties.

δ iso = C -C iso F ∥ C∥ F (A
The deviation is δ iso < 1.9% for prolate pores with aspect ratios γ = 2.0. This implies that even if theoretically one can consider very large number of orientations, in practice, with our present RSA method, a significant effect of the pore shapes remains and lead to an anisotropic response but acceptable.

  Coefficients of piecewise functions f so 11 and f so 33 of axisymmetrical superspheroidal pore . . . . . . . . 4.1 Elastic parameters for the transversely isotropic matrix . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Prescribed overall strain field corresponding to the six load cases . . . . . . . . . . . . . . . . . . . . 4.3 Elastic parameters for the transversely isotropic matrix . . . . . . . . . . . . . . . . . . . . . . . . . A.1 G 0 i jk : numerical results at (x 1 = -1, x 2 = 0.8, x 3 = 1.5) . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 G 0 i j and G 0 i jk : numerical results at (x 1 = 0, x 2 = 0, x 3 = 1.) . . . . . . . . . . . . . . . . . . . . . . . . A.3 Reference transversely isotropic elastic parameters: C 0 i components in TI tensor basis and components of Hill tensor related to an oblate spheroid of aspect ratio γ = 0.5 . . . . . . . . . . . . . . . . . . . . A.4 Number of nodes and elements in the meshes of the cubic models comprising different superspheroid pores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.5 Number of nodes and elements in the meshes of the cubic models comprising different superspherical pores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.6 Evolution of components A i jkl for the superspheroidal pore embedded in a transversely isotropic matrix on the Logarithm value of the concavity parameter log(p) such that p ∈ [0.2, 5] . . . . . . . . A.7 Numerical estimation of A i jkl for the superspherical pore embedded in a transversely isotropic corrected model with different values of concavity p ∈ [0.2, 5]. . . . . . . . . . . . . . . . . . . . . . A.8 Calculation of integrals I k (η) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.9 Calculation of integrals J k (η) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.10 Deviation from isotropy δ iso monodisperse spherical & spheroidal pores embedded in isotropic matrix. 2.8 Geometries of different representative models: (a) the spherical model with ellipsoidal pore; (b) the cubic model comprising the superspherical pore; (c):the cubic model embedded with superspheroidal pore. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  H i jkl for the superspheroidal pore embedded in a transversely isotropic matrix on the Logarithm value of the concavity parameter log(p) such that p ∈ [0.2, 5] . . . . . . . . 2.15 Evolution of components H i jkl for the superspherical pore embedded in a transversely isotropic matrix on the Logarithm value of the concavity parameter log(p) such that p ∈ [0.2, 5] . . . . . . . . . . . . 2.16 Superspherical pore (relation (2.61) with ς = 1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17 2D representation in diametral plane of a 3D axisymmetrical superspheroidal pore, with ς = 1 and symmetry axis x 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  1212 component of the T I compliance contribution tensor of an axisymmetric superspheroidal pore embedded in TI matrix, as a function of concavity parameter p. H 1212 is used to check accuracy of transverse isotropy by comparing to (H 1111 -H 1122 )/2.. . . . . . . . . . . . . . . . . . . . . . . . . . 2.22 Left figure: comparison between H 1212 component and (H 1111 -H 1122 )/2 to illustrate tetragonal symmetry of H tensor. Right figure: relative distance between H tensor and its TI projection as a function of concavity parameter. This relative distance is equal to zero in the case of the spherical inclusion p = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.23 Effective transverse E MTB 1 (top, left) and normal E MTB 3 (top, right) Young's elastic moduli, effective axial shear modulus G MTB 31 as a function of porosity φ, MTB approximation, aligned axisymmetrical superspheroidal pores. App 1 (plain lines): approximation formula (relations 2.70-2.71-2.72) for axisymmetrical spheroidal pores, App 2 (dashed lines) : approximation oblate spheroid with same volume (relations 2.78-2.83). Note that the two approximations coincide only for normal Young's modulus E MTB 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.24 Effective transverse E MTB 1 (top, left) and normal E MTB 3 (top, right) Young's elastic moduli , effective shear coefficient G MTB 31 as a fonction of porosity φ, MTB approximation, superspherical (plain lines) and aligned axisymmetrical superspheroidal (dashed lines) pores. . . . . . . . . . . . . . . . . . . . . 2.25 Effective transverse E ef 1 (top, left) and normal E ef 3 (top, right) Young's moduli, effective shear coefficient G ef 31 (bottom) as a fonction of porosity φ, for aligned axisymmetric superspheroidal pores randomly oriented in the isotropic transverse plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.26 Effective transverse E ef 1 (top, left) and normal E ef 3 (top, right) Young's moduli, effective shear coefficient G ef 31 (bottom) as a fonction of porosity φ, for aligned superspherical pores randomly oriented in the isotropic transverse plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Geometries of different representative models: (a) superspherical pore p = 0.4; (b) superspheroidal pore p = 0.4; (c) spherical pore p = 1.0; (d) aligned penny shaped crack with thickness 0.002mm. . . 3.2 Mesh of the superspherical pore with the concavity p = 0.4 in the center of cuboid matrix . . . . . . 3.3 Relative errors of H 11 with respect to the analytical results for the spherical pore embedded in an isotropic matrix for different scale ratios a/L ∈ [3, 10]. . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Numerical estimation of H 11 for the superspherical pore embedded in transversely isotropic corrected model with different concavity parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10] . 3.5 Numerical estimation of H 33 for the superspherical pore embedded in transversely isotropic corrected model with different concavity parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10] . 3.6 Numerical estimation of H 11 for the superspheroidal pore embedded in transversely isotropic corrected model with different concavity parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10] . 3.7 Numerical estimation of H 33 for the superspheroidal pore embedded in transversely isotropic corrected model with different concavity parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10] .

  of porosity φ, MTB approximation, superspherical pore, with constant anisotropic degree of host matrix κ = 0.2 or concavity parameter p = 0.35. Comparison between approximation formula Eq.(3.52) (plain lines) and oblate spheroid with same volume Eq.(3.57) (dashed lines). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . of porosity φ, MTB approximation, superspheroidal pore, with constant anisotropic degree of host matrix κ = 0.2 or concavity parameter p = 0.35. Comparison between approximation formula Eq.(3.54) (plain lines) and oblate spheroid with same volume Eq.(3.57) and Eq.(3.59) (dashed lines). . . . . . . . . . . . . . . . . . of porosity φ, NIA approximation, with constant anisotropic degree of host matrix κ = 0.2 or concavity parameter p = 0.35. Comparison between superspherical pores (plain lines) and axisymmetric superspheroidal pores (dashed lines). . . . . . . . 3.13 Effective conductivity λ se i (supersphere) and λ so i (superspheroid) as a function of porosity φ, with different approximations, with constant anisotropic degree of host matrix κ = 0.2 and concavity parameter p = 0.35. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 General iteration of microstructure generation based on RSA . . . . . . . . . . . . . . . . . . . . . . 4.2 RVE meshed with 3D elements in NETGEN, number of spherical pores = 50, porosity = 0.2 . . . . . 4.3 Effective elastic constants for five different realizations of monodisperse microstructures with 100 spherical pores (γ = 1) and isotropic matrix, porosity φ = 0.25 . . . . . . . . . . . . . . . . . . . . . 4.4 Effective elastic constants for five different realizations of monodisperse microstructures with 100 spheroidal pores (γ = 2) and isotropic matrix, porosity φ = 0.25 . . . . . . . . . . . . . . . . . . . . 4.5 Effective elastic properties of a material containing various number of spherical pores and isotropic matrix, porosity φ = 0.25 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Effective bulk and shear modulus of a material containing spherical pores and isotropic matrix, with different porosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Effective elastic properties of a material containing various number of spheroid pores and isotropic matrix, aspect ratio a = b, c/a = 2.0, porosity φ = 0.25 . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Effective bulk and shear modulus of a material containing spheroidal pores and isotropic matrix, with different porosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Effective elastic properties of a material containing various number of spherical pores and transversely isotropic matrix, porosity φ = 0.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Effective elastic properties of a material containing spherical pores and transversely isotropic matrix, with different porosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Von Mises Stress field corresponding to different loadcases of a material containing spherical pores and transversely isotropic matrix, porosity φ = 0.25 . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Effective elastic properties of a material containing spheroidal pores γ = 2 and transversely isotropic matrix, with different porosity, compared with spherical case . . . . . . . . . . . . . . . . . . . . . . 4.13 Von Mises Stress field corresponding to different loadcases of a material containing spheroidal pores (γ = 2) and transversely isotropic matrix, porosity φ = 0.25 . . . . . . . . . . . . . . . . . . . . . . . A.1 Compliance contribution tensor for superspherical pore embedded in isotropic matrix with p = 0.4, as dependence of number of nodes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Numerical estimations of A i jkl for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/5 and different scale ratio a/L ∈ [4, 18]. . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.4 Numerical estimations of H i jkl for the spherical pore γ = 1 embedded in a transversely isotropic matrix with different scale ratio a/L ∈ [4, 18]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  a/L ∈ [4, 18]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.6 Numerical estimations of H i jkl for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/2 and different scale ratio a/L ∈ [4, 18]. . . . . . . . . . . . . . . . . . .

  AbbreviationsBSE BackScattered Electron FEM/FEA Finite Element Method/Approach FFT Fast Fourier Transformed GA Gaussian Approximation HSB Hashin Shtrikman Bound MTB Mori-Tanaka-Benveniste NIA Non-Interaction Approximation PBC Periodic Boundary Conditions PCW Ponte Castaneda and Willis r.h.s Right Hand Side RSA Random Sequential Addition RVE/REV, VER Representative Volume Element, Volume Elémentaire Représentatif SEM Scanning Electron Microscopy TI Transversly Isotropic Symbols p Average polarisation tensor D Finite domain E Inhomogeneity domain Ω infinite domain E, Σ Remote strain tensor, Remote stress tensor G, G Green function A Strain concentration tensor B Stess concentration tensor

Figure 1 :

 1 Figure1: Les multiples échelles (" macro-méso-micro") d'homogénéisation[START_REF] Temizer | A numerical method for homogenization in non-linear elasticity[END_REF]

  Fig.2 : laboratoire ANDRA de Meuse-Haute Marne), ou roches de couvertures de cavités de stockage. La compréhension fine du comportement mécanique , en relation avec microstructure complexe est de manière générale cruciale dans le contexte du stockage souterrain.

Figure 2 :

 2 Figure 2: Stockage de déchets radioactifs (http://www.cigeo.com)

Figure 3 :Figure 4 :

 34 Figure 3: Illustration des différentes échelles de roche argileuse Yven et al. (2007)

Figure 1

 1 Figure 1.2: Laboratory ANDRA in Meuse-Haute Marne (France) -storage of the radioactive waste (http://www.cigeo.com)

Figure 1 . 3 :

 13 Figure 1.3: Illustration of different scales of clay rock Yven et al. (2007)

  natural and man-made materials. The observations by backscattered electron images (BSE) of harzburgite in Wark et al. (2003) (Fig. 2.1) and by scanning electron microscopy (SEM) of ironstone in Grgic (2001) (Fig. 2.2) and of oolitic limestone in Grgic (2011) (Fig. 2.3) illustrate the typical concave pores for microstructures of geomaterials.

Figure 2 . 1 :

 21 Figure 2.1: The outlines of melt pools exposed in BSE images of harzburgite showing typical concave pores[START_REF] Wark | Reassessment of pore shapes in microstructurally equilibrated rocks, with implications for permeability of the upper mantle[END_REF]

Figure 2

 2 Figure 2.2: Low magnification SEM images of ironstone showing the concave sideritic interoolithic cement. (a) Crack surface. (b) Polished thin section[START_REF] Grgic | Modélisation du comportement à court et à long terme des roches de la formation ferrifère lorraine[END_REF]

  Figure 2.5: Diagram of infinite inclusion problem under remote strain/stress field

  2.6 as an example).

Figure 2

 2 Figure 2.6: 3D mesh of the cuboid RVE including the superspherical pore with p=0.4, left is general view of half cube cutting in plan Z, right is local map of highlighted region

  )) only few analytical results are available for non-ellipsoidal 3-D shapes (see discussion in the book[START_REF] Kachanov | Micromechanics of materials, with applications[END_REF]).[START_REF] Argatov | Rigid toroidal inhomogeneity in an elastic medium[END_REF] calculated stiffness contribution tensor of an absolutely rigid thin toroidal inhomogeneity,[START_REF] Kachanov | Rice's internal variables formalism and its implications for the elastic and conductive properties of cracked materials, and for the attempts to relate strength to stiffness[END_REF] obtained analytical solutions for compliance contribution tensors of a crack growing from a pore and a cracks with partial contact between the faces, Krasnitskii et al. (2019) evaluated elastic fields associated with a rigid torus. Several results have been obtained combining numerical and analytical techniques.Trofimov et al. (2017b) used finite element calculations to analyse the effect of shape of several representative convex polyhedra on the overall elastic properties of particle-reinforced composites.Trofimov et al. (2017a);Trofimov and Sevostianov (2017) quantified effect of waviness of a helical fiber and its elastic properties on the property contribution tensors of such a fiber. In the context of inhomogeneities of concave shape, analytical approximation of compliance contribution tensor of a superspherical pore was first obtained by Sevostianov and Giraud (2012) using numerical results of[START_REF] Sevostianov | On computation of the compliance and stiffness contribution tensors of non ellipsoidal inhomogeneities[END_REF]. Their result was corrected by[START_REF] Chen | Evaluation of the effective elastic and conductive properties of a material containing concave pores[END_REF] using higher accuracy numerical estimates.Sevostianov et al. (2016a);[START_REF] Chen | Combined effect of pores concavity and aspect ratio on the elastic properties of a porous material[END_REF] considered other types of concave pores.[START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF] analysed elastic fields associated with an inhomogeneity of superspherical pore.

Figure 2

 2 Figure 2.7: Diagram of infinite & bounded domains

  Figure 2.8: Geometries of different representative models: (a) the spherical model with ellipsoidal pore; (b) the cubic model comprising the superspherical pore; (c):the cubic model embedded with superspheroidal pore.

FEM

  simulations that were recently obtained in Chen et al. (2017); Sevostianov et al. (2016a) (see also Trofimov et al.

Figure 2 . 9 :

 29 Figure 2.9: Geometry and mesh of 1/8 spherical inclusion in the center of spherical matrix

Figure 2 .

 2 Figure 2.10: Numerical estimations of H i jkl for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/5 and different scale ratio a/L ∈ [4, 18].

  Figure2.11: Relative errors of numerical computations of H i jkl with respect to the analytical results[START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/5 and different scale ratios a/L ∈[4, 18].

  Figure 2.12: Geometry and mesh for a cubic model comprising a superspherical pore with p = 0.4.

Figure 2

 2 Figure 2.15: Evolution of components H i jkl for the superspherical pore embedded in a transversely isotropic matrix on the Logarithm value of the concavity parameter log(p) such that p ∈ [0.2, 5]

  et al.(2015, 2018);Sevostianov et al. (2016a). The authors supplemented finite element modeling with analytical approximations for compliance contribution tensors of pores of such shapes. These results were used to calculate overall elastic properties of materials with multiple concave pores: oolitic rock[START_REF] Kalo | Effective elastic moduli of a heterogeneous oolitic rock containing 3-d irregularly shaped pores[END_REF] and 3-D printed S i 3 N 4 ceramics[START_REF] Lurie | Mechanical behavior of porous S i 3 N 4 ceramics manufactured with 3d printing technology[END_REF]. All the mentioned results have been obtained for materials with isotropic matrix.The number of explicit results on elastic properties of heterogeneous materials with anisotropic matrix is substantially smaller. Piezoelectric properties of transversely isotropic materials containing circular fibers aligned with the axes of symmetry of the matrix have been calculated using various homogenization techniques by[START_REF] Sevostianov | On the modeling and design of piezocomposites with prescribed properties[END_REF].[START_REF] Sevostianov | Effective elastic properties of matrix composites with transversely-isotropic phases[END_REF] calculated compliance contribution tensor for a spheroidal inhomogeneity of arbitrary aspect ratio embedded in a transversely-isotropic material.[START_REF] Levin | Elastic properties of inhomogeneous transversely isotropic rocks[END_REF] calculated effective elastodynamics properties of transversely isotropic rocks containing aligned spherical and strongly oblate spheroidal pores. Effective porothermoelastic properties of transversely isotropic rocks such as mudstones, argillites, shales have been studied in the frame of Effective Media Theory (EMT), by[START_REF] Giraud | Effective poroelastic properties of transversely isotropic rocks-like composites with arbitrarily oriented ellipsoidal inclusions[END_REF].[START_REF] Cosenza | Effect of the local clay distribution on the effective elastic properties of shales[END_REF] calculated overall properties of transversely-isotropic clay containing spherical inhomogeneities.[START_REF] Vasylevskyi | On micromechanical modeling of orthotropic solids with parallel cracks[END_REF] calculated overall properties of a transversely-isotropic material containing parallel circular cracks.[START_REF] Seyedkavoosi | Randomly oriented cracks in a transversely isotropic material[END_REF]

  Refer to appendices A.1 and A.6 for background on tensors and property contribution tensors. Recently Du et al. (2020) calculated components of the compliance contribution tensors of superspheroidal and axisymmetric superspheroidal pores embedded in a transversely-isotropic material numerically.

Figure 2

 2 Figure 2.16: Superspherical pore (relation (2.61) with ς = 1)

Figure 2

 2 Figure2.17: 2D representation in diametral plane of a 3D axisymmetrical superspheroidal pore, with ς = 1 and symmetry axis x 3

  Figure 2.23: Effective transverse E MTB 1 (top, left) and normal E MTB 3 (top, right) Young's elastic moduli, effective axial shear modulus G MTB 31

  Figure 2.24: Effective transverse E MTB 1 (top, left) and normal E MTB 3 (top, right) Young's elastic moduli , effective shear coefficient G MTB 31

Figure 2 Figure 2

 22 Figure 2.25: Effective transverse E ef 1 (top, left) and normal E ef 3 (top, right) Young's moduli, effective shear coefficient G ef 31 (bottom) as a fonction of porosity φ, for aligned axisymmetric superspheroidal pores randomly oriented in the isotropic transverse plane.

  By taking into account the linearity of Eqs.(3.11) -(3.13) and combining Eqs.(3.14) and (3.15), one has:

  Figure 3.1: Geometries of different representative models: (a) superspherical pore p = 0.4; (b) superspheroidal pore p = 0.4; (c) spherical pore p = 1.0; (d) aligned penny shaped crack with thickness 0.002mm.

Figure 3

 3 Figure 3.2: Mesh of the superspherical pore with the concavity p = 0.4 in the center of cuboid matrix

Figure 3 . 3 :

 33 Figure 3.3: Relative errors of H 11 with respect to the analytical results for the spherical pore embedded in an isotropic matrix for different scale ratios a/L ∈ [3, 10].

Figure 3

 3 Figure 3.6: Numerical estimation of H 11 for the superspheroidal pore embedded in transversely isotropic corrected model with different concavity parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10]

Figure 3 . 7 :

 37 Figure 3.7: Numerical estimation of H 33 for the superspheroidal pore embedded in transversely isotropic corrected model with different concavity parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10]

Figure 3 . 8 :

 38 Figure 3.8: The two independent components H se 11 and H se 33 of a superspherical pore embedded in TI matrix, as a function of concavity parameter p and anisotropic degree of host matrix κ. Comparison between FEM results (star points) and approximate relations (plain surface).

  Figure 3.10: Effective conductivity λ MTB 1 and λ MTB 3 as a function of porosity φ, MTB approximation, superspherical pore, with constant anisotropic degree of host matrix κ = 0.2 or concavity parameter p = 0.35. Comparison between approximation formula Eq.(3.52) (plain lines) and oblate spheroid with same volume Eq.(3.57) (dashed lines).

Figure 3

 3 Figure3.13: Effective conductivity λ se i (supersphere) and λ so i (superspheroid) as a function of porosity φ, with different approximations, with constant anisotropic degree of host matrix κ = 0.2 and concavity parameter p = 0.35.

  compared two approaches to predict the effective elastic properties of solids with regular and irregular shaped pores, by utilizing the numerical homogenization applied to periodic RVEs containing multiple arrangements of pores and the Mori-Tanaka and Maxwell micromechanical models based on the compliance contribution tensor (H-tensor). Contributions of convex polyhedral particles to the overall elastic properties of composite materials are studied in Trofimov et al. (2017b) by adopting the theoretical and numerical homogenization methods. Recently, Trofimov et al. (2018) analyzed the overall elastic properties for the composite materials containing multiple superspherical and superspherical rigid inhomogeneities based on the H tensor. The

Figure 4 . 3 :

 43 Figure 4.3: Effective elastic constants for five different realizations of monodisperse microstructures with 100 spherical pores (γ = 1) and isotropic matrix, porosity φ = 0.25

Figure 4 . 5 :

 45 Figure 4.5: Effective elastic properties of a material containing various number of spherical pores and isotropic matrix, porosity φ = 0.25

(Figure 4 . 6 :

 46 Figure 4.6: Effective bulk and shear modulus of a material containing spherical pores and isotropic matrix, with different porosity

Figure 4 . 7 :

 47 Figure 4.7: Effective elastic properties of a material containing various number of spheroid pores and isotropic matrix, aspect ratio a = b, c/a = 2.0, porosity φ = 0.25

Figure 4

 4 Figure 4.8: Effective bulk and shear modulus of a material containing spheroidal pores and isotropic matrix, with different porosity

Fig. 4 .Figure 4

 44 Fig.4.10 shows the numerical estimation for the effective elastic properties of the concerned porous material as dependence of the porosity, which are also compared with the semi-analytical results obtained from some available

Fig. 4 .

 4 Fig.4.11 shows the distribution of Von Mises stress of one realization in the case of the porosity φ = 0.25 and different prescribed macroscopic strains loads. It shows that the localization of the stress generally occurs in the interaction area of the neighbor pores. The maximal Von Mises stress in the case of the the X/Y loads is higher than that of the Z one, which is also consistent with the transversely isotropic characteristics of the matrix material.

  Figure 4.11: Von Mises Stress field corresponding to different loadcases of a material containing spherical pores and transversely isotropic matrix, porosity φ = 0.25

Fig. 4 .

 4 Figure 4.13: Von Mises Stress field corresponding to different loadcases of a material containing spheroidal pores (γ = 2) and transversely isotropic matrix, porosity φ = 0.25

  (A.25), (A.26) and the derivative of the inverse ( ∂K -1 ∂z

  Hill polarization tensor of a spheroidal inclusion aligned in a transversely isotropic host matrix

  (A.46) Functions I 1 , I 3 (a misprint in relation 5.124 Parnell (2016) has been fixed) write:I 3 (υ i ) = 4π υ i S(υ i γ), I 1 (υ i )

Figure A. 2 :

 2 Figure A.2: Numerical estimations of A i jkl for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/5 and different scale ratio a/L ∈ [4, 18].

  Figure A.3: Relative errors of numerical computations of A i jkl with respect to the analytical results Withers (1989) for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/5 and different scale ratios a/L ∈ [4, 18].

Figure A. 8 :

 8 Figure A.8: Numerical estimations of H i jkl for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/10 and different scale ratio a/L ∈ [4, 18].

  Figure A.10: Evolution of components A i jkl for the superspheroidal pore embedded in a transversely isotropic matrix on the concavity parameter p such that p ∈ [0.2, 5]

  Figure A.11: Evolution of components A i jkl for the superspherical pore embedded in a transversely isotropic matrix on the Logarithm value of the concavity parameter log(p) such that p ∈ [0.2, 5]

  Calculation of integrals I k (η)

  ) = -0.197702p 2 + 0.220097p + 0.939377 f se-b 1122 (p) = -0.24264p 2 + 0.0415192p + 1.0399f se-b 1133 (p) = -0.357454p 2 + 0.52452p + 0.827103 f se-b 3333 (p) = -0.326089p 2 + 0.439731p + 0.861657 f se-b 1212 (p) = -2.04372p 4 + 7.78585p 3 -11.3964p 2 + 7.68465p -0.839061 f se-b 2323 (p) = -2.65601p 4 + 9.6441p 3 -13.4224p 2 + 8.59804p -0.982841 (A.69)Finite element results of an axisymmetrical superspheroidal pore embedded in a TI matrix are given in table 3 of paper[START_REF] Du | Numerical computation of compliance contribution tensor of a concave pore embedded in a transversely isotropic matrix[END_REF]. Approximate relations for the concavity parameter range 0.2 < p < 1 write f so-a 1111 (p) = 54.285p 4 -85.8182p 3 + 44.7893p 2 -6.31294p + 0.294191 f so-a 1122 (p) = 50.5731p 4 -84.3085p 3 + 46.7634p 2 -7.38197p + 0.378348 f so-a 1133 (p) = 39.0659p 4 -42.9843p 3 + 7.43095p 2 + 5.44332p -0.646506 (A.70) f so-b 1111 (p) = -3.74036p 4 + 13.3244p 3 -18.0592p 2 + 11.0818p -1.45743 f so-b 1122 (p) = -4.07865p 4 + 14.7289p 3 -20.227p 2 + 12.5889p -1.82352 f sp-b 1133 (p) = -1.92702p 4 + 6.39843p 3 -7.8095p 2 + 3.95979p + 0.293327 (A.71) A.8. Three dimensional Green function for Poisson's equation in infinite anisotropic medium As shown previously, the corrected boundary condition Eq.(3.11) relies on the gradient of Green function. It is worthy to recall that the three dimensional Green function G 0 (x) for Poisson's equation in infinite anisotropic medium could be analytically calculated by:

Figure

  Figure A.12: Effective elastic constants for five different realizations of monodisperse microstructures with 100 spheroidal pores (γ = 2) and transversely isotropic matrix, porosity φ = 0.25

  Figure A.13: Effective elastic properties of a material containing various number of spherical pores and isotropic matrix, porosity φ = 0.2

Table of Contents

 of Compliance contribution tensor components of the superspherical pores in isotropic matrix . . . . . . 2.2 Compliance contribution tensor components of the superspherical pores in transversely isotropic matrix 2.3 Elastic parameters for the transversely isotropic matrix . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Numerical estimations for the independent components of H tensor in the case of isotropic host matrix embedded with a superspherical pore and the comparison with the results obtained in Trofimov et al. C 0 i jkl components of C 0 tensor and related c i components in transversely isotropic tensor basis E i . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical estimations of H i j for the superspheroidal pore embedded in isotropic host matrix with respect to Sevostianov et al. (2016a) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Numerical estimations of H 11 for the superspherical pore embedded in isotropic host matrix with respect to Chen (2016) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical estimations of H i j for the crack embedded in isotropic host matrix with respect to the analytical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Numerical estimations of H i j for the spherical pore embedded in transversely isotropic host matrix (κ ∈ [0.1, 10]) with respect to the analytical results. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Numerical estimation of H 11 for the superspherical pore embedded in transversely isotropic corrected model with different concavity parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10] . 3.7 Numerical estimation of H 33 for the superspherical pore embedded in transversely isotropic corrected model with different concavity parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10] . 3.8 Numerical estimation of H 11 for the superspheroidal pore embedded in transversely isotropic corrected model with different concavity parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10] . 3.9 Numerical estimation of H 33 for the superspheroidal pore embedded in transversely isotropic corrected model with different parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10] . . . . . . 3.10 Coefficients of piecewise functions f se 11 and f se 33 of superspherical pore . . . . . . . . . . . . . . . . . 3.11 Maximal relative errors of approximate results compared to FEM results, ∥(H Approx

	3.4
	List of Tables
	Résumé Abstract Table of Contents List of Tables List of Figures List of Symbols and Abbreviations Résumé de la thèse 2.1 2.9 Maximal relative errors of approximate relations compared to FEM results, (H Approx i jkl	-H FEM i jkl )/H FEM i jkl	iii v vii xi xiii xviii xx ∞
	i jkl	-H FEM i jkl )/H FEM i jkl	∞
	xi

(2018). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Numerical estimation of H i jkl for the superspheroidal pore embedded in a transversely isotropic corrected model with different values of concavity p ∈ [0.2, 5]. . . . . . . . . . . . . . . . . . . . . . 2.6 Numerical estimation of H i jkl for the superspherical pore embedded in a transversely isotropic corrected model with different values of concavity p ∈ [0.2, 5]. . . . . . . . . . . . . . . . . . . . . . 2.7 Reference transversely isotropic elastic parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Reference transversely isotropic elastic parameters: for the superspherical pore with p ∈ [0.2, 1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Maximal relative errors of approximate relations compared to FEM results, (H Approx for the axisymetrical superspheroidal pore with p ∈ [0.2, 1] . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 linear elasticity and linear thermal conduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2

  [START_REF] Kalo | Effective elastic moduli of a heterogeneous oolitic rock containing 3-d irregularly shaped pores[END_REF]) et céramiques S i3N4 imprimées à 3 -D[START_REF] Lurie | Mechanical behavior of porous S i 3 N 4 ceramics manufactured with 3d printing technology[END_REF]). Leurs résultats ont également été utilisés par[START_REF] Sevostianov | On the compliance contribution tensor for a concave superspherical pore[END_REF] pour dériver des expressions analytiques approximatives pour les composantes du tenseur de contribution de conformité. Cette dérivation a été récemment corrigée et étendue au tenseur de contribution de résistivité[START_REF] Chen | Evaluation of the effective elastic and conductive properties of a material containing concave pores[END_REF];Sevostianov et al. (2016a).

Le tenseur de contribution à la résistivité a été introduit par

[START_REF] Kachanov | Explicit cross-property correlations for porous materials with anisotropic microstructures[END_REF]

dans le contexte de la connexion des propriétés croisées entre les propriétés élastiques et conductrices de matériaux hétérogènes.

[START_REF] Kushch | Dipole moments, property contribution tensors and effective conductivity of anisotropic particulate composites[END_REF]

a développé des relations explicites entre le tenseur de contribution de résistivité et les moments dipolaires.

Trofimov et al. (2017b)

a comparé les prédictions des propriétés élastiques globales des composites renforcés avec des particules de différentes formes polyédriques par des approches d'analyse directe par éléments finis et des schémas micromécaniques. Les résultats des deux approches sont en bon accord pour des fractions volumiques jusqu'à 30% pour toutes les combinaisons de matériaux étudiées. En résumé, les travails mentionnés ci-dessus consistent à étudier/chercher les propriétés effectives du matériau concerné:

  Comme nous l'avons précédemment indiqué, les solutions d'Eshelby sur le problème d'inhomogénéité ellipsoïdale dans une matrice isotrope infinie[START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF]) ont engendré un grand nombre de développements et d'applications pour évaluer les propriétés mécaniques. Le premier problème d'Eshelby est Problème d'inclusion (transformation)[START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF]), qui considère une région d'inclusion dans un solide élastique linéaire homogène infini sollicité par une déformation permanente uniforme due à la déformation non élastique comme dilatation thermique , transformation de phase, etc. La déformation prescrite se réfère à des contraintes propres sans la contrainte de la matrice environnante (déformations de transformation sans contrainte). Il convient de souligner que l'inclusion et

la matrice ont les mêmes constantes élastiques. Le problème d'inclusion d'Eshelby est intéressant pour résoudre les champs élastiques à la fois dans l'inclusion et dans la matrice. Le second problème d'Eshelby est le Problème d'Inhomogénéité (Eshelby (1961)), contrairement au Problème d'Inlcusion, qui a des modules d'élasticité différents de ceux du matrice élastique linéaire. Les champs élastiques du milieu hétérogène soumis à un chargement à distance et à une contrainte propre sont d'un grand intérêt. Les solutions détaillées des problèmes d'Eshelby peuvent être trouvées dans le livre de

  ). Dans le schéma de Maxwell, le champ lointain induit par la présence d'inhomogénéités est assimilé au champ lointain produit par un domaine fictif d'une certaine forme avec des propriétés effectives inconnues. 'agissant les simulations numériques réalisées dans cette thèse, nous nous intéresserons aux deux schémas en modélisation micromécanique respectivement basés sur la solution de tenseurs de contribution en prenant en compte inclusion/pore isolé et des simulations directes appliquées à des microstructures considérées (i.e. VERs correspondants). Plus spécifiquement, au cours de calculs de tenseurs de contribution, nous adopterons une

	Résultats et Conlusions chapitre 3) avec des pores isolants. Une étape importante dans la simulation de tels matériaux anisotropes consiste supersphéroïde, fissure, etc.) qui sont également noyés dans une matrice anisotrope. Un certain nombre d'effets
	à une étude de convergence de l'algorithme numérique. La méthode avec des conditions aux limites adaptées sera supplémentaires pourraient être pris en compte, tels que les interfaces entre la matrice et les pores ainsi que les
	Les propriétés élastiques et thermiques des milieux composites poreux sont régies par les propriétés mécaniques adopté mais nécessite des modifications par rapport au problème élastique linéaire précédemment étudié. De manière instabilités et l'endommagement matérielles (i.e. plasticité, fragilité, etc.) qui pourraient se produire.
	transversalement isotrope réagissant à diverses conditions de chargement : champs de déformations ou de température. de concavité et de l'anisotropie du matériau matriciel, leurs effets combinés sur le tenseur de contribution de propriété Algorithme numérique des constituants, la morphologie des pores, la porosité et l'interaction entre les pores. Comme indiqué ci-dessus, la méthode de champs complet (i.e. l'approche d'homogénéisation numérique) sera adoptée dans le contexte d'homogénéisation semi-analytique. Ceci vise à estimer numériquement les propriétés effectives de matériaux hétérogènes concernés dans cette thèse qui comportent des inhomogénéités non ellipsoïdales concaves et une matrice intéressante, nous allons observer que pour la réponse thermique linéaire, ces simulations convergent plus rapidement que la réponse élastique linéaire. Ceci pourrait être partiellement lié au fait que la solution thermique (e.g. le tenseur de contribution de résistivité, etc) est du second ordre tandis que la solution élastique (e.g. le tenseur de contribution de compliance, etc) est du quatrième ordre. Depuis une analyse couplée après un grand nombre d'étude du paramètre Chapter 1
	Nous allons spécifiquement étudier l'effet de la forme, de la porosité, de la distribution poreuse et de l'orientation des pores sur la réponse thermomécanique à des déformations et températures finies à différentes échelles. Nous allons des constituants concernés seront obtenus. En basant sur des nombreuses données obtenues par les simulations numériques, nous allons adopter des approximations afin de mathématiquement exprimer ces effets aux tenseurs Introduction
	également considérer des microstructures anisotropes périodiques basées sur des algorithmes RSA puis effectuer une de contribution associés. Contrairement au ceux du problème élastique étudié, ils devraient inclure une considération
	étude comparative avec les résultats obtenus par l'approche d'homogénéisation semi-analytique. d'un paramètre supplémentaire. Il s'agit du degré anisotrope de la matrice. Ces expressions approximées seront aussi
	Nous nous concentrons dans un premier temps sur la question de la contribution des inhomogénéités aux utilisées comme entrée dans les différents schémas micromécaniques tels que l'approximation de non interaction, les
	propriétés élastiques effectives. Afin de réaliser une telle tâche, dans le chapitre 2, en adoptant l'approche schémas de Mori Tanaka-Benveniste et de Maxwell, pour analyser les propriétés thermiques effectives. Les deux 1.1. General Introduction
	d'homogénéisation numérique et la méthode des conditions aux limites adaptées, les tenseurs de concentration et variables, l'une servant à la morphologie poreuse, l'autre permettant de faire varier les matériaux de la phase solide
	de contribution seront numériquement calculés avec nos attentions particulièrement portées aux formes concave matricielle, montrent une importance plus significative en particulière si la concavité est plus faible. D'après une
	des inhomogénéités (i.e. formes supersphériques et supersphéroïdales). Ces conditions aux limites adaptées serons analyse paramétrique, on peut conclure que le degré d'anisotropie du matériau matriciel κ et le paramètre de concavité
	rigoureusement reformulées, dans le cas du problème élastique linéaire, via la méthode intégrale basée sur la p présentent un effet combiné significatif lors de l'estimation des coefficients thermiques effectifs. A ce stade, on note
	transformée de Fourier notamment en résolvant le problème de singularité sur l'axe de l'isotropie transverse. La une fois encore que le matériau homogénéisé étudié a une propriété transversalement isotrope ayant le même axe de
	méthode numérique proposée s'est avérée efficace et précise après plusieurs évaluations et validations appropriées symétrie que la matrice.
	en comparant ses prédictions, dans certains cas particuliers, avec les résultats analytiques et certains numériques Le dernier chapitre 4 de cette étude a été consacré à une discussion comparative de deux approches de modélisation
	disponibles. Les prédictions résultantes peuvent converger pour un domaine matriciel relativement petit et le pour prédire les propriétés élastiques effectives : la méthode champs moyen (i.e. méthode d'homogénéisation semi-
	processus prend moins de temps en maintenant une précision suffisamment précise. Les tenseurs de déformations analytique) et la méthode champs complet (i.e. l'homogénéisation numérique directement applique au VER).
	et de contraintes ainsi que ceux de contribution de compliance & résistivité seront calculés par cet algorithme Pour cela, les microstructures sont considérées comme périodiques ayant des pores sphériques et/ou sphéroïdaux
	numérique. Les résultats obtenus quantitativement illustrent l'effet de la concavité des pores et de l'anisotropie du aléatoirement orientés et positionnés dans la phase matricielle. Comme indiqué ci-dessus, ceci sera réalisé par
	matériau matriciel. Nous allons constater que les formes concaves produisent une influence sensiblement plus grande l'algorithme RSA. Les prédictions obtenues par ces deux méthodes seront comparées. Ceci nous permettra d'avoir
	que les formes convexes. Ces effets serons mathématiquement exprimés en adoptant des approximations appropriées certaines conclusions provisoire, telles que, les deux méthode peuvent gérer les spécificités de microstructure
	ensuite consacré à créer le maillage avec qui pourraient être considérées comme une entrée des différents schémas micromécaniques afin de quantitativement périodique. Plus précisément, dans le cas de petites concentrations de pores, l'approche de la NIA fournit précision
	des éléments quadratiques de type C3D10 pour les problèmes élastiques et ceux de type DC3D10 pour les problèmes estimer les propriétés effectives étudiées, tels que les schémas d'homogénéisation NIA, MTB et Maxwell. Une suffisamment pertinente. À mesure que la porosité augmente, les schémas micromécaniques plus avancés tels que
	thermiques). Les simulations ont été effectuées via des calculs parallèles en utilisant le centre de calculs de l'Université grande nombre de résultats seront montrés dans ce travail en portant nos attentions particulières aux application des ceux de MTB et de Maxwell sont de plus en plus nécessaires, ce qui signifie que l'influence de l'interaction entre les
	de Lorraine -Mésocentre EXPLOR (un Ensemble de Calcul Scientifique Pour la Lorraine). L'homogénéisation méthodes précédemment décrites dans le cas de géomatériaux poreux (i.e. roches poreuses, etc). Une observation pores n'est pas négligeable.
	numérique et les procédures des traitement des données sont réalisées par des scripts en Python et Matlab. intéressante est le fait que la réponse élastique effective des matériaux anisotropes dépend fortement de la forme des Notons à nouveau que la méthode numérique proposée dans cette thèse est capable de traiter tous les problème
	pores, de la porosité et de l'interaction des pores. L'effet de la forme du pore et celui de l'anisotropie de la matrice ne avec une anisotropie générale d'un matériau matriciel, mais une telle étude n'a pas été effectuée pour une raison de
	peuvent pas être considérés indépendamment. concentration et de concision du travail. Au niveau des perspectives, une recherche continue pourrait se concentrer
	Par la suite, de nombreux résultats importants seront obtenus dans le cadre de la conduction thermique (voir sur les simulations directes appliquées aux milieux composites contenant des défauts irréguliers (tels que supersphère,

S

méthode récemment développé basée sur des conditions aux limites adaptées

[START_REF] Sevostianov | Elastic fields generated by inhomogeneities: Far-field asymptotics, its shape dependence and relation to the effective elastic properties[END_REF]

afin d'accélérer la convergences des calculs numériques en réduisant le biais dû au modèle borné. La correction des conditions aux limites est donnée en fonction du tenseur de Green et de son gradient comme dépendant de l'élasticité anisotrope du matériau de la matrice. Ceci sont rigoureusement calculés basée sur la transformée de Fourier notamment pour régulariser les singularités sur la axe de l'isotropie transverse. Un grand nombre de simulations ont été réalisées. Les grandeurs clés ici sont le tenseur de contribution de propriété d'une inhomogénéité qui donne le champ élastique supplémentaire produit par l'introduction de l'inhomogénéité due au champ linéaire (i.e. champ élastique ou thermique) appliqué à distance. Différents schémas d'homogénéisation linéaire susmentionnés ont été proposés comme solutions semi-analytiques en adoptant des approximations appropriées. Tous les résultats numériques décrits dans ce travail via les tenseurs de contribution sont directement liés aux propriétés effectives, mais sont obtenus en dehors du cadre de la théorie d'Eshelby.

Concernant les simulations directes (i.e. méthode champ complet), différents outils numériques (e.g. scripts en Matlab, Python, Fortran, Netgen, etc), comme indiqué ci-dessus, ont été développés afin de pouvoir estimer les propriétés effectives liées aux différentes problématiques (e.g. problèmes mécaniques ou thermiques). Tout d'abord, différentes microstructures avec des inclusions/pores sous diverse formes (concave ou non) ont été générées en utilisant des scripts programmés en Matlab et Fortran. Le logiciel Netgen est

Table 2 .

 2 1: Compliance contribution tensor components of the superspherical pores in isotropic matrix

	p	H FEM 1111	H FEM 1122	H FEM 1212	H GA 1111	H GA 1122	H GA 1212	H R.E. 1 1111	H R.E. 1122	H R.E. 1212
	0.3 3.710 -0.869 1.691 3.730 -0.897 1.701 0.53% 3.08% 0.60%
	0.4 2.647 -0.634 1.356 2.639 -0.656 1.371 0.29% 3.29% 1.13%
	0.5 2.320 -0.562 1.270 2.343 -0.565 1.274 0.99% 0.55% 0.31%
	0.7 2.085 -0.505 1.231 2.093 -0.506 1.252 0.40% 0.11% 1.71%
	1	2.003 -0.477 1.240 2.005 -0.477 1.241 0.10% 0.04% 0.06%

1 Relative Error It can be observed in the table 2.1 a good agreement between current results of numerical calculations and the predictions of G.A. in

[START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF]

. The author would like to indicate that, some other comparisons between the present FEM results and those available in literature

[START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF]

;

Chen et al. (

table 2

 2

	.2.

Table 2

 2

.2: Compliance contribution tensor components of the superspherical pores in transversely isotropic matrix

Table 2 .

 2 4: Numerical estimations for the independent components of H tensor in the case of isotropic host matrix embedded with a superspherical pore and the comparison with the results obtained in[START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF].

Table 2

 2 Since the transverse isotropy of the matrix material and that of the superspheroidal pore are both around the e 3 , we can theoretically conclude that the overall response such as the compliance contribution tensor must satisfy the equality (H 1111 -H 1122 )/2 = H 1212 . This property is well verified from the computation results that indirectly demonstrate their accuracy. Fig.2.14 illustrates the H i jkl components as a function of the concavity parameter p. It can be observed that the evolutions, in particular those of H 3333 and H 1313 , are more important when p < 0.5. In other words, the concave form of the porous heterogeneity has a significant effect on the compliance contribution tensor H. Whereas it becomes qualitatively stabled when p > 0.5We investigate the numerical estimation of H i jkl for the superspherical pores surrounded by a cubically bounded matrix whose geometry and mesh is illustrated in Fig.2.12b. Detailed results for different values of concavity p ∈ [0.2, 5] are summarized in Table.2.6. First of all, it should be emphasized that the supersphere, as described in Eq.(2.49), is obviously not transverse isotropic. As a consequence, the aforementioned equality in the case of superspheroidal pores is no longer satisfied. Figure 2.14: Evolution of components H i jkl for the superspheroidal pore embedded in a transversely isotropic matrix on the Logarithm value of the concavity parameter log(p) such that p ∈ [0.2, 5]

	(i.e. convex forms).

Table.2.5 the numerical estimations of each independent non-zero component of H i jkl for the superspheroidal pores with different values of concavity parameter p. .5: Numerical estimation of H i jkl for the superspheroidal pore embedded in a transversely isotropic corrected model with different values of concavity p ∈ [0.2, 5].

2.3.6.2. Superspherical pores

Table 2.6: Numerical estimation of H i jkl for the superspherical pore embedded in a transversely isotropic corrected model with different values of concavity p ∈ [0.2, 5].

Table 2

 2

			.7: Reference transversely isotropic elastic parameters
	E 0 1 (GPa) E 0 3 (GPa)	ν 0 12	ν 0 31	G 0 31 (GPa)
	20.44	11.306		0.1027 0.1798	1.5851
	Table 2.8: Reference transversely isotropic elastic parameters: C 0 i jkl components of C 0 tensor and related c i components in transversely isotropic
	tensor basis E i
	C 0 1111 (GPa) C 0 3333 (GPa)	C 0 1122 (GPa)	C 0 1133 (GPa) C 0 2323 (GPa)
	22.3639	12.9994		3.8275	4.7092	1.5851
	c 1 (GPa)	c 3 (GPa)	c 3 = c 4 (GPa)	c 5 (GPa)	c 6 (GPa)
	12.9994	26.1914		6.65983	18.5363	3.1702
	Numerical results of Du et al. (

Table 2 .

 2 H 1212 , H 2323 . Comparisons between approximate relations (2.69) and finite element results are presented in figure 2.19. Maximal relative errors of approximate relations (2.69) compared to FEM results are given in table 2.9, they are lower than 4.%. 9: Maximal relative errors of approximate relations compared to FEM results, (H

	H E 1111	H E 1122	H E 1133	H E 3333	H E 1212	H E 1313
	0.0349 0.03677 0.0350 0.0356 0.0061 0.0060
					Approx i jkl	-H FEM i jkl )/H FEM i jkl	∞	for the superspherical pore
	with p ∈ [0.2, 1]

Figure

2

.19: The 6 independent components H i jkl of the tetragonal compliance contribution tensor of a superspherical pore embedded in TI matrix, as a function of concavity parameter p. Comparison between FEM results (dashed lines) and approximate relations (plain lines).

  H

	Chapter 2 Chapter 2													2.5 2.5
	Constant ω is determined by fit of finite element results (see section 2.3 and Du et al. (2020) for details on finite H 1111 H 3333 •10 -2 H 1212
	approx FEM element modelings). One obtains 0.12 7.5											approx FEM
	0.12			A = 0.221795, B = 0.591786, ω = 22 3			approx 2	(2.76)
				7
	Analytical solution for the compliance contribution tensor of a spheroidal pore aligned with the directions of a TI
	matrix is recalled in appendix. It may be noticed that approximation (2.71) for H 3333 component numerically coincides 0.11 2
	with the corresponding component of a spheroidal pore with same volume than superspheroidal pore : 6.5
	0.11								1
				H 3333 (p) ≈ H	spheroid FEM 3333 FEM, (H 1111 -H 1122 )/2 (γ(p)), 0.2 ≤ p ≤ 1			(2.77)
	0.1			6
	with	0.2	0.4	0.6 0.2	0.8 0.4	1	0.6	0.2	0.8	0.4		1	0.6	0.8	1
				p					p						p
	γ(p) = The analytical solution for the aligned spheroidal pore H V so (p) V sphere = Γ 1+2p 2p •10 -2 H 1122 -1.8 Γ 3 2p Γ 1 p •10 -2 i jkl -3 (it includes the particular case of the sphere, H 1133 (2.78) spheroid -2 FEM FEM approx approx
	H i jkl sphere -2.2
	-2.4							H 1133
										-4
	-2.6
	-2.8
										-5
	-3
	0.2	0.4	0.6	0.8		1			0.2			0.4		0.6	0.8	1
				p											p
							H 2323
	Maximal relative errors of approximate relations (2.70-2.72-2.77) compared to FEM results are given in table
												FEM
	2.10, they are lower than 5.%.	0.8								approx
				H E 1111	H E 1122	H E 1133	H E 3333	H E 1313
				0.00133 0.00174 0.00297 0.04890 0.03425 0.6
	Table 2.10: Maximal relative errors of approximate relations compared to FEM results, (H superspheroidal pore with p ∈ [0.2, 1]	sphere 2323 Approx i jkl	-H FEM i jkl )/H FEM i jkl	∞	for the axisymetrical (2.75)
				0.4
				0.2
				0.2	0.4		0.6		0.8			1
							52 53	p

with γ(p = 1) = 1) is deduced from the exact Hill tensor recalled in appendix A.6. Comparisons between approximate relations (2.70-2.71-2.72-2.77) and finite element results are presented in figure 2.20. Compliance contribution tensor of aligned axisymmetrical superspheroidal and spheroidal pores of same volume have quasi the same normal component H 3333 . In other words, only the volume characterizes this component, not affected by concavity (superspheroid) of convexity (spheroid). It must be emphasized that it is not the case for all the other components H i jkl , including the shear component H 2323 , for which the concavity parameter p is of major importance (volume is not sufficient to characterize compliance contribution tensor). Shear component H 1212 in the plane of transverse isotropy has been used to check accuracy of the symmetry of revolution by comparing to H 1111 -H 1122 )/2.

(see figure

2

.21), it may be noticed that both coincide as expected.

Figure

2

.20: The 5 independent components H i jkl of the T I compliance contribution tensor of an axisymmetric superspheroidal pore embedded in TI matrix, as a function of concavity parameter p. Comparison between FEM results (dashed lines) and approximate relations (plain lines). Note that H 1212 is used to check accuracy of transverse isotropy by comparing to (H 1111 -H 1122 )/2.. Figure 2.21: H 1212 component of the T I compliance contribution tensor of an axisymmetric superspheroidal pore embedded in TI matrix, as a function of concavity parameter p. H 1212 is used to check accuracy of transverse isotropy by comparing to (H 1111 -H 1122 )/2..

  Transverse isotropic projection of compliance contribution tensor Π T I H E 0 will be used instead of the compliance contribution of the superspherical pore H E 0 . See appendix for detail, only components H 1111 , H 1122 and H 1212 are modified, other components are equal, Left figure: comparison between H 1212 component and (H 1111 -H 1122 )/2 to illustrate tetragonal symmetry of H tensor. Right figure: relative distance between H tensor and its TI projection as a function of concavity parameter. This relative distance is equal to zero in the case of the spherical inclusion p = 1.

	H T I 1133 = H 1133 , H T I 3333 = H 3333 , H T I 2323 = H 2323 (H T I i jkl denotes Π T I H E 0 i jkl ) H T I 1111 = 3 H 1111 + H 1122 + 2 H 1212 4 , H T I 1122 = H 1111 + 3 H 1122 -2 H 1212 4 H T I 1212 = H 1111 -H 1122 + 2 H 1212 4 Figure (2.22) illustrates the numerical difference between the tetragonal tensor H E 0 and its TI projection Π T I H E (2.81) (2.82) 0 . 0.2 0.4 0.6 0.8 1 0.1 0.15 0.2 0.25 0.3 0.35 p H 1212 supersphere: H tetragonal H 1212 (H 1111 -H 1122 )/2 0.2 0.4 0.6 0.8 1 0 0.02 0.04 0.06 0.08 0.1 p ∥H -Π T I (H)∥/∥H∥ 2.5 2.5.4.1. Aligned axisymmetric superspheroidal pores Effective elastic coefficients E MTB 1 , E MTB 3 , G MTB 31 obtained with MTB approximation are presented in figure (2.23). Approximation formula (2.70-2.71-2.72) deduced from FEM are compared to approximation of compliance contribution tensor of an oblate spheroidal pore of same volume (semi axis length of axisymmetrical superspheroid is equal to the greater semi axis length of the oblate spheroid) ): H i jkl (p) ≈ H spheroid i jkl (γ(p)), 0.2 ≤ p ≤ 1 (2.83) with aspect ratio of oblate spheroid γ(p) defined by relation (2.78) It may be observed that: • normal Young's modulus E MTB 3 may be estimated by a very simple approximation using compliance Figure 2.22: Chapter 2 contribution of an oblate spheroidal pore of same volume than the axisymmetrical superspheroidal pore, in the

concavity range p < 0.5. Numerically it is mainly related to normal component of compliance contribution tensor H 3333 . Comparison between axisymmetrical superspheroidal and oblate spheroidal pores shows that both approximations lead to the same effective coefficient E MTB 3 . In other words, concavity or convexity has no significant effect on this coefficient, when comparing same pores of same volume. It must be emphasized that this result is very specific and cannot be generalised. It only holds for this particular shape and the normal Young's modulus.. • oppositely, transverse Young's modulus E MTB 1 and shear coefficient G MTB 31 are strongly related to the concavity parameter p. Comparison of estimates based on concave and convex pores of same volume (respectively axisymmetrical superspheroid and oblate spheroid), in the range p < 0.5, see figure (2.23), shows significant

differences. This result is expected as an approximation based on an oblate spheroid of same volume is not precise for all components H 1111 , H 1122 , H 1133 , H 2323 , particularly in the concave range 0.2 < p < 0.5. It confirms, for a 3D shape embedded in an anisotropic matrix that the concavity parameter is of major importance when estimating effective elastic properties.

Table 3 .

 3

1: linear elasticity and linear thermal conduction

Table 3 .

 3 2: Numerical estimations of H i j for the superspheroidal pore embedded in isotropic host matrix with respect toSevostianov et al. (2016a)

	p	H FEM 11	H FEM 33	H (2016) 11	1	H (2016) 33	1 R.E.(H 11 ) 2 R.E.(H 33 ) 2
	0.2	1.4027 14.8636 2.0247 15.1972 44.345%	2.244%
	0.25	1.3486	6.4986	1.7156		6.6319	27.212%	2.052%
	0.3	1.3655	3.8970	1.6400		3.9392	20.102%	1.083%
	0.35	1.3900	2.8087	1.5486		2.8234	11.412%	0.522%
	0.4	1.4108	2.2794	1.5290		2.2896	8.377%	0.450%
	0.5	1.4410	1.8303	1.5063		1.8322	4.528%	0.106%
	0.7	1.4750	1.5828	1.5091		1.5880	2.310%	0.330%
	1	1.4998	1.4998	1.5012		1.4963	0.095%	0.238%
	p = 1 theory	1.5	1.5	1.5		1.5
	R.E. 3	0.012% 0.012% 0.083% 0.250%
	1 Numerical results obtained in Sevostianov et al. (2016a)
	2 Relative error of FEM results with respect to those of Sevostianov et al. (2016a)
	3

Table 3 .

 3 3: Numerical estimations of H 11 for the superspherical pore embedded in isotropic host matrix with respect to[START_REF] Chen | Micromechanical Modelling of Heterogenous Porous Materials with Application to Oolitic Rocks[END_REF]

	p	H FEM 11	H Chen 1 11	H R.E. 2 11
	0.2	3.9379	2.8465 27.715%
	0.25	2.6193	2.1166 19.191%
	0.3	2.0892	1.7603 15.743%
	0.35	1.8385	1.6689	9.226%
	0.4	1.7065	1.6421	3.776%
	0.5	1.5856	1.5719	0.860%
	1	1.4999	1.5111	0.749%
	3	1.5531	1.5588	0.369%
	1 theory 1.5000	1.5000
	R.E. 3 0.008% 0.740%
	1 Chen (2016)
	2 Relative error of numerical H 11 with
	respect to Chen (2016)
	3 Relative error of numerical H 11 with
	respect to analytical one for p = 1

Table 3 .

 3 4: Numerical estimations of H i j for the crack embedded in isotropic host matrix with respect to the analytical resultsModelsTheory FE crack Fig:3.1d FE superspheroid p = 0.2 Ref. superspheroid 1

	H 11	0	0.0015	0.0598	0.0863
	H 33	0.6366	0.6403	0.6336	0.6478
	Relative Error H 33	NA	0.582%	0.483%	1.751%
	1 Normalised results of Superspheroidal case with p = 0.2 in Sevostianov et al. (2016a)

Table 3 .

 3 5: Numerical estimations of H i j for the spherical pore embedded in transversely isotropic host matrix (κ ∈ [0.1, 10]) with respect to the analytical results.

	κ	H FEM 11	H FEM 33	H theory 11	H theory 33	H R.E. 1 11	H R.E. 33
	0.1 12.1317 2.8453 12.1322 2.8450 0.004% 0.013%
	0.2 6.4179 2.2625 6.4184 2.2626 0.007% 0.005%
	0.5 2.7985 1.7518 2.7988 1.7519 0.008% 0.008%
	0.8 1.8333 1.5710 1.8335 1.5712 0.008% 0.008%
	1	1.4999 1.4999 1.5000 1.5000 0.008% 0.008%
	2	0.8022 1.3269 0.8023 1.3271 0.008% 0.008%
	5	0.3467 1.1813 0.3468 1.1813 0.008% 0.008%
	8	0.2238 1.1321 0.2239 1.1322 0.008% 0.008%
	10	0.1815 1.1133 0.1815 1.1134 0.008% 0.008%
	1 Relative error of numerical H 11 with respect to analytical one
				78

Table 3 .

 3 6: Numerical estimation of H 11 for the superspherical pore embedded in transversely isotropic corrected model with different concavity parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10]Table 3.7: Numerical estimation of H 33 for the superspherical pore embedded in transversely isotropic corrected model with different concavity parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10]

	H 11

1 Relative error of numerical H 11 with respect to analytical one for p = 1

Table 3 .

 3 8: Numerical estimation of H 11 for the superspheroidal pore embedded in transversely isotropic corrected model with different concavity parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10] Table3.9: Numerical estimation of H 33 for the superspheroidal pore embedded in transversely isotropic corrected model with different parameters p ∈ [0.2, 5] and anisotropic degrees of matrix κ ∈ [0.1, 10] Relative error of numerical H 33 with respect to analytical one for p = 1 Figs. 3.6 and 3.7 show the evolution of the H i j components on the concavity parameter p and the one of the matrix anisotropy κ, As same as the superspherical cases, it can also be finally concluded that both of them significantly affect the resistivity contribution tensor H E 0 .

	H 11

  se ii (κ, p) = L 00 + L 10 κ + L 01 p + L 20 κ 2 + L 11 κp + L 02 p 2 +L 30 κ 3 + L 21 κ 2 p + L 12 κp 2 + L 03 p 3 +L 40 κ 4 + L 31 κ 3 p + L 22 κ 2 p 2 + L 13 κp 3 + L 04 p 4

	(3.53)
	The polynomial coefficients L mn (m, n ∈ [0, 4]) of f se 11 and f se 33 up to degree 4 are presented in Table.3.10 with different
	piecewise ranges of κ and p.

Table 3

 3

			.10: Coefficients of piecewise functions f se 11 and f se 33 of superspherical pore
	coef		f se 11				f se 33
	κ ∈	[0.1, 1.0]	[1.0, 10]	[0.1, 1.0] [1.0, 10]	[0.1, 1.0]	[1.0, 10]	[0.1, 1.0] [1.0, 10]
	p ∈	[0.2, 0.5]	[0.2, 0.5] [0.5, 1.0] [0.5, 1.0] [0.2, 0.5]	[0.2, 0.5] [0.5, 1.0] [0.5, 1.0]
	L 00	1.07427	1.45602	0.24384	0.37308	2.39465	1.85791	0.65934	0.16292
	L 10	1.35187	0.12840	0.08821	0.00591	-1.96658	-0.30303 -0.75780 -0.03036
	L 01	-8.46522	-9.53368	1.22295	0.99061 -12.17791 -11.15540 0.81355	1.87068
	L 20	-1.34853	-0.01234 -0.01261 -0.00018	1.93006	0.04303	0.33645	0.00144
	L 11	-4.07483	-0.40572 -0.07860 -0.00414	3.96823	0.85177	1.05374	0.04605
	L 02	36.03311	36.13293 -0.46815 -0.36658 39.31755	39.88421 -0.73064 -1.51874
	L 30	0.50130	0.00049	0	0	-0.81052	-0.00325	0	0
	L 21	3.48496	0.02659	0	0	-2.54793	-0.06165 -0.34665 -0.00151
	L 12	3.10502	0.43796	0	0	-3.40873	-0.95001 -0.28565 -0.01501
	L 03	-57.47065 -55.48038	0	0	-57.16757 -61.00594 0.25638	0.48429
	L 40	0	0	0	0	0	0	0	0
	L 31	-1.01925	-0.00083	0	0	1.04211	0.00189	0	0
	L 22	-1.66488	-0.01069	0	0	0.01614	0.02762	0	0
	L 13	-0.13651	-0.17788	0	0	2.37955	0.38600	0	0
	L 04	33.03312	31.62765	0	0	31.21625	35.19119	0	0
	Comparisons between approximate relation Eq.(3.52) and finite element results are presented in figure 3.8.

Table . 3

 . .11. The small errors less than 0.96% are negligible. Satisfactory agreement between proposed approximate formula and numerical results are observed for different superspherical shapes and thermal conductivities.

	Chapter 3			3.3
	H se 11	H se 33	H so 11	H so 33
	0.951% 0.871% 0.500% 1.489%
		92

Table 3 .

 3 11: Maximal relative errors of approximate results compared to FEM results, ∥(H Similar approximation formula Eq.(3.54) is proposed for components H 11 (H 22 = H 11 ) and H 33 of axisymmetrical superspheroidal pore but with another set of geometric parameters.

	Approx ii	-H FEM ii	)/H FEM ii	∥

Table 3

 3

			.12: Coefficients of piecewise functions f so 11 and f so 33 of axisymmetrical superspheroidal pore
	coef		f so 11				f so 33
	κ ∈	[0.1, 1.0]	[1.0, 10]	[0.1, 1.0] [1.0, 10]	[0.1, 1.0]	[1.0, 10]	[0.1, 1.0]	[1.0, 10]
	p ∈	[0.2, 0.5]	[0.2, 0.5] [0.5, 1.0] [0.5, 1.0]	[0.2, 0.5]	[0.2, 0.5]	[0.5, 1.0] [0.5, 1.0]
	L 00	0.26483	0.31506	-0.25095 -0.33687	19.08281	11.76301	4.33631	1.29115
	L 10	-0.00238	-0.01251 -0.05148 0.00924	-12.97547	-1.45203	-4.57413 -0.31561
	L 01	-4.87821	-5.68135	3.96451	4.28190 -141.94940 -84.24440 -10.88093 0.41270
	L 20	0.01106	0.00069	0.02711 -0.00055	12.82125	0.17599	3.30432	0.02697
	L 11	-0.47228	0.06239	0.11204 -0.00771	49.14392	5.02357	10.33640	0.71128
	L 02	36.00467	39.01536 -4.30156 -4.67035 463.93148	268.14315	14.16878 -2.63909
	L 30	0	0	0	0	-8.36487	-0.01212	-1.08391 -0.00095
	L 21	0.11231	-0.00302 -0.02921 0.00054	-22.09062	-0.30483	-4.87185 -0.03998
	L 12	1.47851	-0.05392 -0.05868 -0.00141 -80.78645	-6.96489	-7.98576 -0.54316
	L 03	-70.65161 -74.94268 1.58820	1.72562 -704.18395 -394.21125 -8.68698	3.02476
	L 40	0	0	0	0	2.32501	0.00034	0	0
	L 31	0	0	0	0	5.51144	0.00850	1.14473	0.00102
	L 22	-0.21038	0.00223	0	0	12.74062	0.16059	1.47231	0.01196
	L 13	-1.15256	-0.00330	0	0	51.70313	3.68239	2.26321	0.15184
	L 04	45.74778	47.86995	0	0	407.85909	221.70354	2.05902	-1.09374

Table 4 .

 4 2: Prescribed overall strain field corresponding to the six load cases

			load case
	Components of prescribed strain
	1	2	3	4	5	6

  (E 1 , E 2 , E 3 ), Poisson's ratio (ν 12 , ν 13 , ν 23 ) and shear moduli (G 12 , G 13 , G 23 ) of materials containing different number of spherical pores (n = [1

  ), Ponte Castaneda and Willis (PCW)[START_REF] Ponte Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF]) one, MTB and Maxwell ones make the same predictions in the considered porous media with isotropic matrix and spherical pores (e.g. E 1 ≃ E 2 ≃ E 3 ).

	0.7						0.7
	0.68						0.68
	0.66						0.66
	0.64						0.64
	0.62						0.62
	0.6						0.6
	0.58						0.58
	0.56						0.56
	0.54						0.54
	0.52						0.52
	0.5						0.5
	0	50	100	150	200	250	300	0	50	100	150	200	250	300
				1
				0.98
				0.96
				0.94
				0.92
				0.9
				0.88
				0.86
				0.84
				0.82
				0.8
				0	50	100	150	200	250	300

Table . 4

 . .3, the transverse isotropy of the matrix is described by the following parameters: Firstly, the convergence analysis of transversely isotropic matrix in terms of the number of pores is conducted for the microstructures containing randomly distributed spherical pores with the porosity φ = 0.20. As shown in Figs 4.9, Figure4.9: Effective elastic properties of a material containing various number of spherical pores and transversely isotropic matrix, porosity φ = 0.2

	4.4.1. Transversely isotropic matrix with randomly spherical pores
	0.8						0.8
	0.75						0.75
	0.7						0.7
	0.65						0.65
	0.6						0.6
	0.55						0.55
	0.5						0.5
	0	50	100	150	200	250	300	0	50	100		150	200	250	300
				1.8
				1.7
				1.6
				1.5
				1.4
				1.3
				1.2
				1.1
				1
				0.9
				0.8
				0	50	100	150		200	250	300
					E 1 (GPa) E 3 (GPa)	ν 12	ν 31		G 13 (GPa)
			Properties of material	20.44	11.31		0.1027 0.1798	1.585
				Table 4.3: Elastic parameters for the transversely isotropic matrix

  10)Relations with usual H i jkl components writeh 1 = H 3333 , h 2 = H 1111 + H 1122 , h 3 = √ 2H 3311 , h 4 = √ 2H 1133 , h 5 = H 1111 -H 1122 , h 6 = 2H 2323(A.11)

Table A .

 A 1: G 0 i jk : numerical results at (x 1 = -1, x 2 = 0.8, x 3 = 1.5) 12k413.858 -737.743 1271.82 Considering the data given in table (2.3) and position vector (x 1 = -1, x 2 = 0.8, x 3 = 1.5) one obtains 10 7 G 0 =

		k = 1	k = 2	k = 3
	10 7 g 11k 871.471 -2504.54 -24616.9
	10 7 g 22k 2910.25 -520.84 -24044.6
	10 7 g 33k 68988.8 -55191. -14305.5
	10 7 g 23k 5133.61 2897.27 -1858.36
	10 7 g 31k 587.147 5133.61	2322.95
	10 7 g 			
	        	39800.47 -903.6818 -7004.161 -903.6818 39393.815 5603.329 -7004.161 5603.329S 134599.77	        	(A.43)

Table A .

 A 2: G 0 i j and G 0 i jk : numerical results at (x 1 = 0, x 2 = 0, x 3 = 1.)

		k = 1	k = 2	k = 3
	10 7 g kk (no summation) 60722.101 60722.101 502034.39
	10 7 g 11k	0.	0.	-60722.101
	10 7 g 22k	0.	0.	-60722.101
	10 7 g 33k	0.	0.	-502034.39
	10 7 g 23k	0.	70648.848	0.
	10 7 g 31k	70648.848	0.	0.
	10 7 g 12k	0.	0.	0.
	lim x→1 S(x) = 1

Table A .

 A 3: Reference transversely isotropic elastic parameters: C 0 i components in TI tensor basis and components of Hill tensor related to an oblate spheroid of aspect ratio γ = 0.5 GPa -1 ) p 3333 (GPa -1 ) p 1122 (GPa -1 ) p 1133 (GPa -1 ) p 2323 (GPa -1 )Numerical results are given in table (A.3).A.3. Information of meshes for the FEM computations in the case of cubic model containing concave poreWe show in Fig.A.1 the analysis of mesh convergence in terms of the number of nodes for one example, the compliance contribution tensor for superspherical pore embedded in isotropic matrix with p = 0.4. It can be observed that the numerical results converge when the number exceeds 1.5 million.

	1.486												10.8
													10.7
	1.484	c 0 1111 (GPa)		C 0 3333 (GPa)	C 0 1122 (GPa) 10.6		C 0 1133 (GPa)	C 0 2323 (GPa)
	1.482			22.3639		12.9994		3.8275	10.5		4.7092	1.5851
				c 1 (GPa)		c 2 (GPa)	c 3 = c 4 (GPa) 10.4		c 5 (GPa)	c 6 (GPa)
	1.48
				12.9994		26.1914		6.65983	10.3		18.5363	3.1702
	4 1.478	6	p 1 (GPa -1 ) 8	10	p 2 (GPa -1 ) 12 14	p 3 = p 4 (GPa -1 ) 16 18 4 10.2	p 5 (GPa -1 ) 6 8	10 p 6 (GPa -1 ) 12	14	16	18
			0.0656356		0.0199043	-0.0100943		0.0292588	0.106492
	-0.115	(a) A 1111 p 1111 (0.0245816	0.0656356	-0.00467723 3.15		-0.00713773	0.0532462 (b) A 3333
	-0.12												3.1
	and then -0.125												3.05
	-0.13					P 1111 =	p 2 + p 5 2	, P 1122 =	p 2 -p 5 2 3	, P 3333 = p 1	(A.51)
	-0.135												2.95
	4	6	8	10 (c) A 1133	P 1133 = 12	p 4 √ 2 14	, P 3311 = 16 18	p 3 √ 2	, P 2323 = 4 6	p 6 2	8	10 (d) A 3311	12	14	16	18 (A.52)
	-0.08												1.575
	-0.082
													1.57
	-0.084
	-0.086												1.565
	-0.088
	-0.09												1.56
	4	6	8		10	12		14	16	18			4	6	8	10	12	14	16	18
					(e) A 1122											(f) A 1212
							2.085
							2.08
							2.075
							2.07
							2.065
							2.06
							4	6	8	10	12		14	16	18

(g) A 1313

•

  Spherical pore: γ = 1 Numerical estimations of H i jkl for the spherical pore γ = 1 embedded in a transversely isotropic matrix with different scale ratio a/L ∈ [4, 18]. Relative errors of numerical computations of H i jkl with respect to the analytical results[START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] for the spherical pore γ = 1 embedded in a transversely isotropic matrix with different scale ratios a/L ∈ [4, 18]. Numerical estimations of H i jkl for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/2 and different scale ratio a/L ∈ [4, 18]. Relative errors of numerical computations of H i jkl with respect to the analytical results[START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/2 and different scale ratios a/L ∈ [4, 18].

		0.12 10 0 0.11 10 0 0.0634		0.23 10 0 0.45 10 0 1.65
		0.0633
		0.118 0.1		10 -1 0.4
	Relative error Relative error	0.114 0.116 10 -4 10 -2 0.0629 10 -4 0.08 0.09 10 -2 0.063 0.0631 0.0632	Relative error Relative error	0.22 0.225 10 -3 10 -2 10 -4 0.3 0.35 10 -2 1.55 1.6
	4 0.112 4 -0.03 -0.029 -0.028 -0.027 -0.026 -0.025 4 0.069 0.07 0.071 0.072 0.073 0.074 0.075 Figure A.4: 140 6 8 10 12 14 16 18 (a) H 1111 6 8 10 12 14 16 18 (c) H 1122 6 8 10 12 14 16 18 (e) H 1212 143 4 6 8 10 12 14 16 18 0.07 0.0628 10 -6 (a) H 1212 4 6 8 10 12 14 16 18 4 6 8 10 12 14 16 18 Figure A.5: 2 (a) H 1111 5 10 15 -0.025 -0.02 -0.015 -0.01 (c) H 1122 5 10 15 0.04 0.045 0.05 0.055 0.06 0.065 (e) H 1212 4 6 8 10 12 14 16 18 (a) Figure A.7: • Oblate pore: γ = 1/10 (a) H 1111 4 6 8 10 12 14 16 18 -6.96 -6.94 -6.92 -6.9 -6.88 -6.86 -6.84 10 -3 (c) H 1122 4 6 8 10 12 14 16 18 0.0348 0.0349 0.035 0.0351 0.0352 Figure A.6: 2 (e) H 1212	4 0.215 -0.024 4 -0.0235 -0.023 -0.0225 4 0.216 0.218 0.22 0.222 0.224 10 -4 1.5 0.25 4 2 4 -0.035 -0.03 -0.025 -0.02 2 0.2 0.21 0.22 0.23 0.24 0.25 2 -0.037 -0.036 -0.035 -0.034 -0.033 4 4 6 4 0.488 0.49 0.492 0.494 0.496 0.498	6 6 4 4	6 6 6 6 5 8 6 6 6	8 8	10 (b) H 3333 10 (d) H 1133 10 (f) H 1313 10 (b) H 1313 8 10 10 (b) H 3333 10 12 12 12 12 12 12 (d) H 1133 8 10 12 (f) H 1313 8 10 12 (b) (b) H 3333 8 10 12 (d) H 1133 8 8 8 10 12 (f) H 1313	14 14 14 14 14 14 14 14 14 14	16 16 16 16 16 16 15 16 16 16 16	18 18 18 18 18 18 18 18 18 18

Table A .

 A 6: Evolution of components A i jkl for the superspheroidal pore embedded in a transversely isotropic matrix on the Logarithm value of the concavity parameter log(p) such that p ∈ [0.2, 5]

				2.4092 1.8562
	0.35 2.2013 -0.3019 -0.0646 2.0944	8.0402 2.5032 1.5821
	0.4 2.2674 -0.3064 -0.0218 1.4539	6.0956 2.5737 1.4733
	0.45 2.3159 -0.3101 0.0106	1.0948	4.9690 2.6259 1.4255
	0.5 2.3494 -0.3117 0.0355	0.8777	4.2632 2.6610 1.4019
	0.6 2.3910 -0.3139 0.0670	0.6557	3.5146 2.7046 1.3865
	0.7 2.4121 -0.3140 0.0856	0.5505	3.1416 2.7258 1.3870
	0.8 2.4228 -0.3134 0.0974	0.4942	2.9321 2.7358 1.3933
	0.9 2.4281 -0.3126 0.1054	0.4616	2.8035 2.7403 1.4018
	1	2.4312 -0.3122 0.1112	0.4417	2.7200 2.7431 1.4112
	1.5 2.4300 -0.3090 0.1265	0.4123	2.5554 2.7386 1.4534
	2	2.4281 -0.3081 0.1340	0.4152	2.5211 2.7356 1.4858
	2.5 2.4267 -0.3074 0.1388	0.4234	2.5159 2.7338 1.5098
	3	2.4272 -0.3075 0.1425	0.4321	2.5203 2.7342 1.5285
	4	2.4290 -0.3075 0.1476	0.4470	2.5342 2.7361 1.5551
	5	2.4314 -0.3078 0.1512	0.4584	2.5478 2.7385 1.5731

Table A .

 A 7: Numerical estimation of A i jkl for the superspherical pore embedded in a transversely isotropic corrected model with different values of concavity p ∈ [0.2, 5].

	p	A FEM 1111	A FEM 1122	A FEM 1133	A FEM 3311	A FEM 3333	A FEM 1212	A FEM 1313
	0.2 11.0342 -1.5470 0.4701 1.6971 11.5571 7.5756 3.1521
	0.25 6.9234 -1.0207 0.2974 1.1060 7.3017 4.7760 2.1095
	0.3	5.1440 -0.7928 0.2216 0.8517 5.4701 3.6464 1.7214
	0.35 4.1841 -0.6706 0.1799 0.7189 4.4972 3.1302 1.5559
	0.4	3.6067 -0.5955 0.1537 0.6410 3.9260 2.8884 1.4771
	0.45 3.2362 -0.5423 0.1364 0.5910 3.5668 2.7761 1.4360
	0.5	2.9741 -0.4975 0.1243 0.5499 3.3036 2.7236 1.4132
	0.6	2.7249 -0.4420 0.1130 0.5078 3.0561 2.6957 1.3944
	0.7	2.5883 -0.3976 0.1086 0.4781 2.9065 2.6973 1.3912
	0.8	2.5102 -0.3629 0.1078 0.4595 2.8163 2.7094 1.3951
	0.9	2.4623 -0.3351 0.1089 0.4482 2.7586 2.7256 1.4024
	1	2.4314 -0.3123 0.1111 0.4417 2.7202 2.7435 1.4113
	1.5	2.3776 -0.2394 0.1262 0.4389 2.6479 2.8318 1.4578
	2	2.3733 -0.1983 0.1402 0.4500 2.6375 2.9019 1.4953
	2.5	2.3816 -0.1723 0.1513 0.4623 2.6441 2.9613 1.5252
	3	2.3911 -0.1535 0.1602 0.4734 2.6531 3.0068 1.5483
	4	2.4095 -0.1283 0.1735 0.4912 2.6718 3.0745 1.5817
	5	2.4230 -0.1114 0.1829 0.5042 2.6860 3.1190 1.6040
	• Superspherical pores

  Finite element results of a superspherical pore embedded in a TI matrix, are given in table 4 of paper[START_REF] Du | Numerical computation of compliance contribution tensor of a concave pore embedded in a transversely isotropic matrix[END_REF]. Approximate relations for the concavity parameter range 0.2 < p < 1 writef se-a 1111 (p) = 2.72394p 2 -2.6248p + 1.63141 f se-a 1122 (p) = 1.23332p 2 -1.20253p + 1.29293 f se-a 1133 (p) = 3.75077p 2 -3.50795p + 1.81628 f se-a 3333 (p) = 2.77572p 2 -2.46747p + 1.5398 f se-a 1212 (p) = 17.1508p 4 -43.7713p 3 + 38.3878p 2 -13.1282p + 2.36576 f se-a 2323 (p) = 50.3296p 4 -90.0984p 3 + 59.3072p 2 -15.8086p + 2.1938

	1133 p 1 2 + 2	√	2 C 3333 C 1133 p 1 3	, h c 6 = -	1 2323 p 1 4 C 2 6	(A.67)
	A.7. Coefficients for approximation formula of
	compliance contribution tensors

Table A .

 A 10: Deviation from isotropy δ iso monodisperse spherical & spheroidal pores embedded in isotropic matrix.

	Porosity Random spherical pores	Random spheroidal pores
		+ ISO matrix	+ ISO matrix
	0.05	0.07%	0.54%
	0.1	0.16%	1.24%
	0.15	0.32%	1.75%
	0.2	0.17%	1.60%
	0.25	0.13%	1.82%

This comparison is made except for the gradient G 0 i j,k on the symmetry axis due to the singularity of the reference solution derived from[START_REF] Pouya | Green's function solution and displacement potentials for transformed transversely isotropic materials[END_REF].

Note also that the corrected boundary conditions can also be applied to any type of anisotropy. It probably requires supplementary simulations and is not further pursued here.

The property contribution tensors then turn into conductivity N and resistivity H contribution tensors in contrast with the elasticity, the reader is referred for instance to[START_REF] Du | Numerical computation of compliance contribution tensor of a concave pore embedded in a transversely isotropic matrix[END_REF] or or section

2.3

It is important to emphasize that the proposed approach can be applied to any form of the bounded media and any type of matrix anisotropy.

It is worthy to emphasize that the so-called important anisotropy is defined in the present work with a small value of κ and in the direction of e 3 , which is aligned with the symmetry of the concave pore

It is worthy to mentioned here that due to the meshing difficulties in the case of concave pores, the superspherical and superspheroidal form will not be concerned at present. This could be overcome by using parametric equations or the FFT based numerical homogenization but such a study is not our pursue here.

Acknowledgements

Acknowledgements JFB thanks Cerema and Université Gustave Eiffel for allowing the creation of the joint research team entitled " Équipe de Recherche Commune sur les Matériaux pour une Construction Durable (ERC MCD)", within which these research works were carried out. IS acknowledges financial support from National Aeronautics and Space Administration (NASA) Cooperative Agreement NNX15AL51H.

reveals that the effect of the matrix anisotropy κ seems to be more important than that of the pore concavity p on the component H 11 . However, unlike the previous observation, Fig. 3.5 illustrates that the matrix anisotropy plays an important role in the estimation of the component H 33 , which is, as expected, logical due to the collinearity between the geometrical symmetry and the one of the matrix anisotropy.

Superspheroidal pore in transversely isotropic matrix

In this section, we investigate the combined effect of the concavity of the superspheroidal pore and the matrix anisotropy on the resistivity contribution tensor H E 0 , whose non-zero and independent components H 11 and H 33 are computed and summarized in Tables 3.8 and 3.9, respectively. Similar to the previous superspherical case discussed in Section 3.2.5.1, the estimations of both H 11 and H 33 components decrease with the increase of the anisotropy parameter κ, and present a non-monotone evolution due to the effect of the concavity parameter p. The turning points are also marked with yellow box in the tables and indicate a certain correlation with the anisotropy degree of matrix κ.

Even this non-monotonicity is slight, it is worthy to note here that, unlike the superspherical pores discussed before, it occurs when the pore is "extremely concave" (e.g. p ≃ 0.25) for the component H 11 and when the pore is "severely convex" for H 33 (e.g. p ≥ 1). A possible interpretation is that, unlike the superspherical cases, the superspheroidal pores present a symmetry of revolution that is colinear with the direction of the matrix anisotropy, which induces an augmentation of the macroscopic anisotropic response of the RVE in the direction of e 3 . function, see [START_REF] Chen | Evaluation of the effective elastic and conductive properties of a material containing concave pores[END_REF]; Sevostianov et al. (2016a); [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF]; [START_REF] Du | Effective elastic properties of transversely isotropic materials with concave pores[END_REF] for details)

where superscripts se and so respectively refer to supersphere and superspheroid, p denotes the concavity parameter.

Except for some particular values of concavity parameter (p = 1 4 , 1 2 , 1), the total surface area needs to be calculated by numerical integration. As in [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF]; [START_REF] Du | Effective elastic properties of transversely isotropic materials with concave pores[END_REF], the surface area A so (p) of the axisymmetrical superspheroid is given by the single integral accounting symmetry of revolution:

The projection area S proj (p) of both 3D shapes onto planes x i x 3 (i = 1, 2, x 3 denotes symmetry axis of the axisymmetrical superspheroid) writes (Beta denotes Euler Beta function):

(3.51)

Approximation formula for superspherical pore

In this section, we consider a possibility to involve the geometric parameters presented in section 3.3.2.1 to yield the analytical formula of the resistivity contribution tensor H for superspherical inhomogeneities of different shapes and matrix of different contrast degree in thermal conductivities. We restrict this study to the range 0.2 ≤ p ≤ 1.0 and 0.1 ≤ κ ≤ 10. After analysing different approximate results interrelated the geometric parameters, we conclude the most accurate representations for superspherical case, which is summarized as follows:

where sphere denotes spherical case (particular case when p = 1). The non-dimensional factors are based on those proposed by [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF]. The analytical representation of limiting case of sphere (p = 1) are used for components H 11 and H 33 (with x 3 symmetry axis of TI matrix ). The analytical solution for the aligned spheroidal pore which includes the particular case of the sphere, is deduced from the exact Hill tensor recalled in appendix A.9. Functions f se ii (κ, p) are given in Eq.(3.53) by fitting the numerical results obtained in section 3.2.5.1. They are piecewise functions in polynomial forms related to both variables κ and p considering both concave and convex

• Both the anisotropic degree of host matrix κ and concavity parameter p have significant importance when estimating effective thermal coefficients. It's not possible to separate effect of anisotropy from the effect of concavity. Changes in variables κ and p have a greater impact on effective λ MTB 3 than λ MTB 1 which mean that the effect in axial conductivity aligned with the symmetry axis of matrix is larger than the transversal ones.

• As porosity φ increases, effective conductivity gradually decreases. This is completely logical because the thermal conductivity of insulating inhomogeneities tends to be zero. The increase in volume proportion of inhomogeneities will increase its impact.

Effective thermal properties of aligned axisymmetric superspheroidal pores with MTB scheme

Effective thermal coefficients λ MTB 1 and λ MTB 3 obtained with MTB approximation of aligned axisymmetric superspheroidal pores are presented in Fig. 3.11. Approximation formula Eq.3.54 deduced from FEM are compared to approximation of resistivity contribution tensor of an oblate spheroidal pore of same volume (the semi axis length of superspheroid is equal to the maximum axis radius of the oblate spheroid) Eq.(3.57) and Eq.(3.59).

Concluding remarks

The main novelty of present work is the prediction of effective thermal properties of heterogeneous materials containing 3D irregular pore embedded in an anisotropic matrix by consideration of a wide range of concavity parameter and matrix anisotropy. The thermal analytical approximation formulas for the resistivity contribution tensor of the reference shapes, supersphere and axisymmetrical superspheroid are provided based on the FEM data, which are then used as input in the different micromechanical schemes such as NIA, MTB and Maxwell schemes, to analyse the effective thermal properties. The current approach takes into account both the anisotropy of microstructural shape and the contrast conductivity of host matrix. The advantage of the proposed method is that it allows us to consider composites with any number of concavity parameter, contrast degree and porosity in the range of 0.2 ≤ p ≤ 1.0, 0.1 ≤ κ ≤ 10 and ϕ ≤ 0.3 with the explicit expressions. On the basis of parametric analysis, one can conclude that both the anisotropic degree of host matrix κ and concavity parameter p have significant importance when estimating the effective thermal coefficients. The different shapes and anisotropic degrees of matrix have different impact on the effective conductivity in transverse and perpendicular plane. The obtained effective porous material is transversely isotropic with same symmetry axis than matrix.

Chapter 4

Numerical homogenization applied to periodic RVEs

This chapter deals with the numerical estimations of the overall elastic properties of the porous geomaterials by adopting the Finite Element Method (FEM) based numerical homogenization directly applied to its microstructure, which is represented by a periodic Representative Volume Element (RVE) being able to characterize all of its necessary information. After an introduction in Section 4.1, the Random Sequential Adsorption (RSA) based approach is briefly recalled in Section 4.2 for the numerical generation of the periodic RVEs with a finite number of randomly positioned and oriented pores. It is then followed in Section 4.3 with some numerical results as well as their comparison to some analytical ones available in the literature and the semi-analytical predictions obtained from the previously proposed micromechanical modeling (see section 2.4). Section 4.4 presents some novel numerical results of composites containing random arrangement of spherical or spheroidal pores embedded in transversely isotropic medium. We finally present some concluding remarks in Section 4.5.

obtained results were compared with the direct FEM simulations performed on periodic RVEs containing parallel oriented pores as well as rigid particles.

However, all of the above-mentioned numerical results were obtained in the case of isotropic matrix. To the author's knowledge, pure analytical results are only available in some particular cases, such as those of ellipsoidal inhomogeneities embedded in isotropic matrix, ellipsoidal inhomogeneities aligned in transversely isotropic matrix with the same direction of symmetry, etc. This work aims to predict the effective elastic properties with different shapes of pores randomly embedded in an anisotropic solid by adopting the numerical homogenization method (i.e.

full field simulations). The obtained results are devoted to the assessment of the semi-analytical models obtained in the previous chapters (see Chapter 2).

Microstructure description and algorithm of direct FEA simulations

The microstructures studied in this work are numerically obtained based on the random sequential addition (RSA) algorithm as described in [START_REF] Torquato | Random sequential addition of hard spheres in high Euclidean dimensions[END_REF]; [START_REF] Anoukou | Random distribution of polydisperse ellipsoidal inclusions and homogenization estimates for porous elastic materials[END_REF]. As shown in Fig. 4.1, the general idea of RSA algorithm is that, during the process of the microstructure generation, the current candidate pore will be accepted if it does not overlap any previously defined pores. Otherwise, it will be rejected and a new candidate pore will be randomly generated until attaining the desired porosity. The volumes of the phases are

and with volume fraction φ = V (r) /V.

4.2

The aspect ratios of different families of ellipsoidal pores are γ r 1 = c r /a r and γ r

In our work, we consider pores of mono-disperse but same shape which means the parameters are equal to the reference one:

Hence, the lengths of the semi-axes of the reference pore are obtained by the following equations:

Table .4.1 gives some examples of the lengths of the semi-axes of the reference pores for different shapes: Following [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF], periodic boundary conditions are employed in the analysis. The displacement field u in such a unit cell (i.e. RVE) can be expressed by Eq.4.2

where x denote the microscopic position of a point in the unit cell. E gives the overall strain field which would be the actual applied homogeneous strain and u * is the periodic displacement field which accounts for the fluctuation due to the presence of heterogeneities. All components of periodic u * takes the identical values at the point on the boundary of the unit cell. Since the local strain field ε(u(x)) derives from u(x), it can be split as shown in Eq.4.3.

Note that the periodicity of u * decides the naught of the average of ε * on the unit cell and therefore the average strain field is equal to the overall applied strain field:

with ⟨.⟩ representing the average of a field in a domain Ω of volume V, such that:

Moreover, the Finite Element Method is adopted to resolve the microscopic stress and strain fields that will be used in the numerical homogenization process. Note here that the microstructure (RVE) will be subjected to uniaxial tensile and shearing loadings for respectively studying the corresponding macroscopic elastic responses of the concerned porous geomaterials. Ten-nodes quadratic tetrahedral elements (i.e. C3D10 in ABAQUS notation) are used in this thesis for the meshing of the RVE that are realized via NETGEN (Schöberl (1997)) as shown in Fig. 4.2.

where L is the length of the cubic RVE, k = (1, 2, 3) denotes the degree of freedom and E k j represents the applied displacement. These PBCs could be automatically generated by NETGEN with several necessary declaration in the corresponding input file of the RVE geometry that are compatible with the parameter setting of the ABAQUS/Standard software (see [START_REF] Wu | Applying Periodic Boundary Conditions in Finite Element Analysis[END_REF]; [START_REF] Omairey | Development of an ABAQUS plugin tool for periodic RVE homogenisation[END_REF]). As shown in Table .4.2, six sets of kinematically admissible fields (i.e. the displacement and strain fields) are applied to the RVE for the computations in the cases of uniaxial tension and three shear loads along different directions. Note that the prescribed macroscopic strains are set to be enough small (e.g. 1.0e -5 ) that the volumes could be considered constant.

Once the numerical computations are performed, the result files are processed by using a user-defined Matlab Chapter 4

4.3 dispersion. Hence, in the following part of this section, we will adopt the H tensor to estimate effective modulus of materials containing randomly oriented spheroidal pores as a comparison. In the framework of the NIA, such a study can be realized by referring the following relationships [START_REF] Wu | The effect of inclusion shape on the elastic moduli of a two-phase material[END_REF]):

( summation over i, j = 1, 2, 3) (4.10)

where T denotes the Wu's strain concentration tensor [START_REF] Wu | The effect of inclusion shape on the elastic moduli of a two-phase material[END_REF]). The relationship among T, H and the compliance tensor of the matrix material S 0 is written as H = T : S 0 [START_REF] Sevostianov | Relations between compliances of inhomogeneities having the same shape but different elastic constants[END_REF]). K and G are respectively the pore compressibility and the pore shear compliance (two invariants related to the shape rather than the orientation).

Note here that the explicit expressions for the parameters K and G in the cases of spherical and randomly oriented spheroidal pores are given in [START_REF] David | Compressibility and shear compliance of spheroidal pores: Exact derivation via the Eshelby tensor, and asymptotic expressions in limiting cases[END_REF].

Following the Mori-Tanaka scheme, one has:

Note also that [START_REF] Sevostianov | Generalization of Maxwell homogenization scheme for elastic material containing inhomogeneities of diverse shape[END_REF] re-writes Maxwell's schemes with property contribution tensors H and N. [START_REF] Sevostianov | On the shape of effective inclusion in the Maxwell homogenization scheme for anisotropic elastic composites[END_REF] gives some advice on the choice of the aspect ratios of the effective inclusion. One can conclude that in the limiting cases of totally randomly oriented inhomogeneities(resp. perfectly parallel identical inhomogeneities), the effective inclusion in Maxwell can be taken as a sphere (resp.a representative inhomogeneity).

Hence,in the case of a microstructure (i.e. RVE) with randomly oriented spheroidal pores, the effective inclusion could be considered spherical.

The number of pores 100 is chosen in following numerical computations according to the convergence analysis as shown in Fig. 4.7. Direct FEA simulation is applied to the microstructures with different porosities (e.g. φ = 0.05, 0.10, 0.15, 0.20, 0.25) and monodispersed spheroidal pores (with aspect ratio defined as 2). As illustrated in Fig. 4.8, the estimations of the effective bulk and the shear modulus are compared to those obtained from the micromechanical homogenization methods. The comparison is also made by taking into account the direct FEA results reported in [START_REF] Drach | Comparison of full field and single pore approaches to homogenization of linearly elastic materials with pores of regular and irregular shapes[END_REF].

Chapter 5

Conclusions and perspectives

This work was devoted to studying the effective elastic and thermal properties of porous composite geomaterial by taking into account the anisotropy simultaneously introduced by its complex microstructure and the matrix material.

In order to keep the work focused and concise, porous media has been specifically studied by paying our particular attention to the effect of the pore morphology and that of the transversely isotropic character of the matrix phase.

Superspherical and axisymmetric superspheroidal concave pores were taken to be two candidates representing the microstructural complexity.

We first focused on the issue of the extra contribution by inhomogeneities to the effective elastic properties.

In order to achieve such a task, in Chapter 2, the concentration and contribution tensors have been numerically estimated that are of direct relevance with the effective elastic properties of a transversely isotropic solid embedded by a non-ellipsoidal inhomogeneity, as mentioned above, by considering the superspherical and superspheroidal forms, respectively. These two concave shapes are in direct relation with the tomographic observation of the concerned geomaterials microstructures. In contrast to the infinite media that is usually considered in the classical and theoretical framework, the numerical computations were carried out by adopting the Finite Element Method based homogenization in a so-called bounded domain thanks to the adapted boundary conditions based method, initially proposed by [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF] in the case of isotropic matrix. By considering the corresponding Green functions and their gradients, such kind of boundary conditions have been rigorously reformulated in this thesis via the Fourier transform based integral method in particular by solving the singularity problem on the axis of the transverse isotropy.

This allows to accelerate the computation convergence without sacrificing its accuracy by incorporating the matrix anisotropy and the correction of the bias induced by the bounded character of the mesh domain. By adopting this method, the strain & stress concentration tensors and the compliance & stiffness contribution tensors have then been computed using a much smaller matrix domain. The results have quantitatively illustrated the effect of pore concavity and matrix anisotropy and showed that the concave shapes produce substantially larger influence than the convex ones. Moreover, the resulting predictions can be found to converge for a relatively small matrix domain and 123 Chapter 5

5.0 the process has been shown to be less time consuming by holding a sufficiently accurate precision. This was also proved after several appropriate assessments and validation by comparing its predictions, in some particular cases, with analytical results and some available numerical ones. Subsequently, based on the obtained numerical results, the contribution tensors were approximated in terms of the pore concavity parameter p in the range 0.2 ≤ p ≤ 1 for the two reference shapes: supersphere and axisymmetrical superspheroid. The accuracy is found to be less than 5% for all of their components. These semi-analytical formulas were next adopted as the input of the different micromechanical schemes for the predictions of the corresponding effective elastic properties, such as Non Interaction Approximation (NIA), Mori-Tanaka-Benveniste (MTB) and Maxwell ones. It was qualitatively and quantitatively found that the effective elastic response of anisotropic materials is highly dependent on the pore shape, porosity and the interaction of pores, whose combined effect was also shown to be as significant.

Secondly, effective thermal properties of porous media with anisotropic matrix and a concave insulating pore were studied in Chapter 3. As similar as in the elastic problem, the aforementioned method of adapted boundary conditions was reformulated and applied in order to first estimate the resistivity contribution tensor of the concerned geomaterial.

A first interesting observation is that the computations converge much quicker than those of the previously studied elastic problem, which should be due to the fact that the resistivity contribution tensor is in second order while the compliance contribution tensor is in fourth order. Coming to the consideration of the complex microstructure and the matrix anisotropy, a large number of computations have been realized leading to the conclusion that both of them significantly affect the thermal conductivity of the solid, which was found to be of critical importance when the pore is concave and the matrix is relatively much anisotropic. Next, based on these numerical estimations, some appropriated approximations were applied to the development of semi-analytical functions of the resistivity contribution tensor, by simultaneously considering the above-mentioned two parameters, the concavity of pore and the matrix anisotropy.

They were then adopted in the semi-analytical homogenization modelling, for which different micromechanical schemes such as NIA, MTB and Maxwell ones have been taken into account to study the corresponding effective thermal properties. It revealed that both the pore concavity and the matrix anisotropy have significant importance on the effective thermal coefficients especially at their lower values and, perform different influences on the transverse and perpendicular planes. The obtained effective homogenized material is transversely isotropic with the same symmetry axis as that of the matrix material.

Additionally, Chapter 4 was devoted to the application of the full field simulations based on the direct numerical homogenization on the periodic microstructure. Note that only spherical and spheroidal pores are taken into account due to the meshing difficulties of the RVE in the case of finite number of concave pores. The periodic RVE was generated by using the algorithm based on the Random Sequential Adsorption (RSA) method [START_REF] Torquato | Random Heterogenous Materials[END_REF][START_REF] Lopez-Pamies | The nonlinear elastic response of suspensions of rigid inclusions in rubber: Ii -a simple explicit approximation for finite-concentration suspensions[END_REF][START_REF] Cheng | A homogenization model for porous ductile solids under cyclic loads comprising a matrix with isotropic and linear kinematic hardening[END_REF][START_REF] Zerhouni | Quantifying the effect of two-point correlations on the effective elasticity of specific classes of random porous materials with and without connectivity[END_REF]) that allows to control some physical parameters during the numerical modelling, such as the porosity, the orientation (random or not) of the pores, the number of the pores, etc. By referring to the concerned geomaterials (e.g. clay rocks), the porosity is approximately defined as less or equal to 30%. The resulting predictions were then compared with those obtained by the method of the semi-

Appendix

A.1. Background on tensors

Notations : Barred letters A, C, D, Q refer to fourth order tensors, bold letters ε, σ, i refer to second order tensors, underlined letters x, x refer to first order tensors. Einstein's summation convention over repeated indices is used unless otherwise indicated. ⊗, : and :: respectively represent tensor product, (dot product), contracted products on two and four indices. i, I, J and K = I -J respectively represent the second-rank identity tensor, the fourth-rank symmetric identity tensor, and fourth-rank spherical and deviatoric isotropic projectors (δ i j denotes Kronecker delta symbol, δ i j = 1 if i = j, δ i j = 0 otherwise).

It may be interesting to introduce standard notation and the corresponding simplified algebra for fourthorder transversely isotropic tensor (see [START_REF] Walpole | Fourth-rank tensors in the thirty-two crystal classes: multiplication tables[END_REF]). See also post of Sébastien Brisard on github, http://sbrisard.github.io/, intitled Decomposition of transverse isotropic, fourth-rank tensors. By denoting n the components c i and p i of matrix stiffness tensor C 0 and strain Hill polarisation tensor P E 0

and components h i of compliance contribution tensor are deduced from calculation rules in transversely isotropic basis

See similar derivations in [START_REF] Sevostianov | Effective elastic properties of matrix composites with transversely-isotropic phases[END_REF] which considers different TI basis (see formula 2.56).

A fourth order tensor T may be transversely isotropised by projection onto a transversely isotropic basis to obtain

In case of a symmetry axis is x 3 , we apply the transformation Q α representing a rotation of angle φ about x 3 axis to fourth order tensor T and non zero components of matrix Q α write

Transverse isotropic projection Π T I (T) (A.17) corresponds to the following average over orientations, in the transverse The general expression of the Green tensor derived from a reasoning based on the Fourier transform [START_REF] Mura | Micromechanics of defects in solids[END_REF] or plane-wave expansion [START_REF] Willis | Micromechanics determination of electroelastic properties of piezoelectric materials containing voids[END_REF] writes

where δ is the scalar Dirac distribution and K(ξ) is the acoustic tensor defined from the fourth order elastic tensor C by

The expression (A.24) can be further simplified by using the spherical coordinates associated to the pole x 0 in which the integration variable ξ writes with the variable change z = cos θ ξ = cos θ e + sin θ u φ = z e + 1z 2 u φ (A.26)

where e = x/ ∥x∥ and u φ is the unit vector parametrized by φ describing the unit circle centered on 0 and orthogonal to x.

Recalling that δ(λ z) = δ(z)/λ for all λ > 0, it finally comes that Bonnet ( 2009)

. Gradient of the Green tensor

Taking the gradient of (A.24) provides the following third order tensor grad G

In the particular case investigated

One defines 6 functions R α and R * α (α = 1, 2, 3)

In the non degenerate case c

As indicated in [START_REF] Pouya | Green's function solution and displacement potentials for transformed transversely isotropic materials[END_REF], the singularity on x 3 axis (ρ = 0) may be removed by using the following expression for w 3

Green function's for the infinite transversely isotropic medium writes

g 1 δ i j + g 2 x i x j + g 3 x i δ 3 j + x j δ 3i + g 4 δ 3i δ 3 j (A.40)

with

is straightforward by using formal calculation tools. We provide in Tables A.4 and A.5 the mesh information during the FEM computation for cubically bounded representative elementary volume containing respectively the superspheroidal and the superspherical pore. The displayed numbers of nodes and those of 3D quadratic elements shows that each mesh is well refined for the corresponding computations to obtain a precision of computation as accurate as possible. for the ellipsoidal pore with the aspect ratio γ = 1/5. As can be found in the computation results of H, similar evolution trends in function of the scale ratio a/L, can be observed. It should be note here that, unlike the compliance contribution tensor H, the strain concentration tensor A does not present the major symmetry between the components A 1133 and A 3311 . As a consequence, it has 7 non-zero independent components.

Additionally, we show in Fig. A.3 the relative errors of their numerical estimations with respect to the analytcial solutions proposed by [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF] as well as those obtained from the classical numerical computations without the correction of the boundary conditions. A very good precision can be found for the proposed numerical procedure in the case of smaller scale ratios a/L ∈ [4, 8], whereas the correction of the boundary conditions is shown to be useless when a/L is sufficiently important.

Property contribution tensors are used in micromechanics to describe the contribution of a single inhomogeneity to the property of interest [START_REF] Kachanov | Micromechanics of materials, with applications[END_REF]). Compliance contribution tensors have been first introduced in the context of pores and cracks by [START_REF] Hori | Overall moduli of solids with microcracks: Load induced anisotropy[END_REF] as the extra average strain produced by a pore. The average strain, over representative volume | Ω |, can be represented as a sum

where S 0 where is the compliance tensor of the matrix and Σ is uniform remotely applied stress. The material is assumed to be a linear elastic; hence, the extra strain ∆ε due to presence of an inhomogeneity E is a linear function of the applied stress:

where | E | is the pore volume and H E 0 is fourth-rank compliance contribution tensor of the pore. The H E 0 tensor is determined by the shape and size of the inhomogeneity, as well as properties of the matrix and of the inhomogeneity material. This tensor is also affected by elastic interactions. In the non-interaction approximation, it is taken by treating the inhomogeneities as isolated ones. In the case of multiple inhomogeneities, the extra strain produced by m-th inhomogeneity is ∆ε (m) = f (m) H E(m) 0 : Σ so that the extra compliance due to all the inhomogeneities is given by

Formula (A.55) highlights the fundamental importance of the compliance contribution tensors: these tensors have to be summed up and averaged in the context of the effective elastic properties. The sum

properly reflects compliance contributions of individual inhomogeneities and constitutes the general microstructural parameters in whose terms the effective compliance should be expressed. Components of this tensor were calculated for 2 -D pores of various shape and 3 -D ellipsoidal pores in isotropic material by [START_REF] Kachanov | Effective Moduli of Solids With Cavities of Various Shapes[END_REF]. For

the general case of elastic inhomogeneities, these tensors were introduced and calculated (for ellipsoidal shapes) by Sevostianov and Kachanov (Sevostianov and Kachanov (2002a)). Components of the compliance contribution tensor We recall the compact solution of strain Hill polarization tensor of a spheroidal inclusion aligned in a transversely isotropic host matrix recently presented in [START_REF] Barthélémy | Simplified approach to the derivation of the relationship between hill polarization tensors of transformed problems and applications[END_REF] (reader may refer to this reference for the detailed derivation and python script of the complete solution). The corresponding compliance contribution tensor may be deduced from relations presented in section (A.1). In what follows, the aspect ratio of the spheroidal inclusion is denoted ω, and the symmetry axis n of the host matrix and the spheroidal inclusion is taken as n = e 3 . The particular case of the spherical pore is then deduced by setting ω = 1. The solution writes in Walpole TI tensor basis

γ 1 and γ 2 denote the square roots of the roots of the polynomial

The square root of a complex argument is defined with a positive real part and it is consistent with the cmath library of Python. The inverse hyperbolic cosine of a complex argument, denoted by arcosh , has one branch cut, extending left from 1 along the real axis to -∞, continuous from above. Note that arcosh Z = ln (Z + √ Z 2 -1) where the principal

In particular isotropic case which means λ i 0 = λ 0 (i = 1, 3), the equations Eqs.(A.72) -(A.73) reduce to

(A.75)

A.9. Hill polarization tensor and resistivity contribution tensor of a spheroidal inclusion aligned in a transversely isotropic host matrix

One considers a transversely isotropic matrix of conductivity tensor λ 0 (n denotes unit vector on the symmetry axis, in this paper n = e 3 ):

The Hill polarisation tensor P E 0 of an spheroidal inclusion aligned in the directions of a transversely isotropic matrix (i.e. spheroid and matrix have the same symmetry axis) writes:

where γ is the aspect ratio of radius and with shape function g(ξ) (see [START_REF] Barthélémy | Effective permeability of media with a dense network of long and micro fractures[END_REF]; [START_REF] Giraud | Effective electrical conductivity of transversely isotropic rocks with arbitrarily oriented ellipsoidal inclusions[END_REF])

By inserting (A.77) into (3.2) and with the relationship (3.2), resistivity contribution tensor H E 0 of an insulating λ E = 0 aligned spheroidal pore writes:

i N (A.79) and the particular case of the spherical pore γ = 1 embedded in transversely isotropic matrix is derived

i N (A.80) and the particular case of the spherical pore γ = 1 embedded in an isotropic matrix ν = 1 is recovered

For the limiting case of an aligned penny shaped crack embedded in a TI matrix, the resistivity contribution tensor writes

A.10. RSA algorithm [START_REF] Anoukou | Random distribution of polydisperse ellipsoidal inclusions and homogenization estimates for porous elastic materials[END_REF] propose an extension of the RSA algorithm to obtain random periodic (or not) distributions of mono-and polydisperse non-overlapping ellipsoidal inclusions of arbitrary shape and orientation. Here we briefly recall the algorithm and please refer to [START_REF] Anoukou | Random distribution of polydisperse ellipsoidal inclusions and homogenization estimates for porous elastic materials[END_REF] for detail. The inputs of the algorithm are: the dimensions of the cuboidal cell L 1 , L 2 and L 3 , the volume fraction f , the number of phases N p , the number N re f of reference inclusions, the microstructural parameters ω (r) 1 , ω (r) 2 , χ (r) , f r , and two offsets distance ξ 1 and ξ 2 used in the calculation of the minimum distance parameters s 1 and s 2 , respectively. The outputs of the algorithm are: the position vector of the center of the ellipsoid v r i , the semi-axes lengths a (r) i , b (r) i , c (r) i and the Euler angles ϕ (r) i , θ (r) i , ψ (r) i of the ellipsoidal inclusion i belonging to the phase r. Specifically, for each phase r = 2, ..., N p , the RSA algorithm can be decomposed into four steps:

• Step 1: Compute the semi-axis lengths of inclusion (i).

•

Step 2: In the sequential addition, generate a random center position vector v (r) i for inclusion i in phase r. Compute the minimum distance ∆ 1 between a new inclusion i and any previously accepted inclusion j = 1, ..., i-1 including its 26 periodic images, and compare this distance to the minimum value s 1 .

• Step3: Use the algorithm for finding the minimum distance between ellipsoid and plane to determine ∆ 2 and ∆2 . Compare the distance ∆ 2 with the minimum value s 2 . If ∆ 2 ≤ s 2 , then generate a new center position of inclusion i, otherwise accept the inclusion i.

• Step4: Ensure periodicity of the cuboidal cell by considering periodic images of the inclusion.