Prof Jean Sulem

Prof Mahdia Hattab

Besuelle to attend my final defense but also for letting my defense be an enjoyable moment,

List of Figures

XV

List of Tables

L'analyse des images fournies par ces techniques a montré que ces roches sont formées de trois constituants principaux :

• Les oolithes ayant une forme plus ou moins sphérique.

• Les pores inter-oolithiques ayant des formes irrégulières.

• Le ciment inter-oolithique. Pour le problème de conductivité thermique, le produit vectoriel correspond au produit du vecteur de température et du vecteur normal.

• Méthode analytique basée sur la fonction explicite de l'ellipsoïde Pour vérifier l'approximation des pores par des ellipsoïdes, on a évalué ensuite les tenseurs de contribution de souplesse et de résistivité thermique pour les ellipsoïdes avec la même procédure numérique d'une part, et en utilisant une solution analytique en calculant le vecteur normal à partir de la fonction implicite de l'ellipsoïde 1 2 3 ( , , ) 0 f z z z  . Le vecteur normal unitaire à la surface de l'inhomogénéité est alors proportionnel au gradient de la fonction scalaire. L'indice «0» se réfère à la phase matricielle et l'indice «1» se réfère à la partie d'inclusion. Ainsi, l'importance du premier problème Eshelby est due à la détermination du tenseur de contribution de la souplesse en fonction des tenseurs Hill et Eshelby, c'est-à-dire pour l'inclusion ellipsoïdale, la solution des tenseurs de souplesse peut être obtenue analytiquement comme suit :

ℍ = -𝕊 0 ∶ ℕ ∶ 𝕊 0 (0.6)
où N est un tenseur constant de quatrième ordre dépendant de la forme et donné par: Bruggeman, 1935 ;Hill, 1965). Elle permet la transition de l'échelle microscopique à l'échelle mésoscopique.

ℕ -1 = (ℂ 1 -ℂ 0 ) -1 + ℙ (0.

XXVIII

La deuxième étape permet la transition de l'échelle mésoscopique à l'échelle macroscopique par le schéma d'homogénéisation de Maxwell. Le schéma d'homogénéisation de Maxwell (Maxwell, 1873) a été introduit pour la première fois dans le contexte de la conductivité électrique effective d'un matériau contenant plusieurs inhomogénéités sphériques. Maxwell a proposé une solution pour ce problème en comparant les résultats du champ lointain de la perturbation du champ électrique appliqué de l'extérieur calculé de deux façons différentes:

-Le champ lointain résultant est la somme de tous les champs lointains générés par toutes les inhomogénéités individuelles dans une région homogénéisée Ω.

-Le champ résultant à des points éloignés est égal au champ lointain généré par la grande région homogénéisée Ω avec des propriétés effectives inconnues.

General Introduction

In general, the short and long-term behavior of rocks is modeled either by macro-mechanical (i.e., phenomenological) approaches or by micromechanical approaches. Micromechanical approaches have a particular and practical advantage since a relatively small number of assumptions are needed. In this case, the complexity of the simulated macroscopic behavior does not derive from the complexity of the introduced formalism, but rather from a few microstructural assumptions and statistical considerations on the constitutive elements (grains, crystals, pores, cracks).

The macroscopic thermomechanical behavior of the materials is largely affected by the microstructure. The micro-macroscopic relationship is described via micromechanical models developed within the framework of the homogenization theory of heterogeneous media (Maxwell homogenization method for example). The methods of homogenization (scale change) are used to define effective properties, elastic or conductive, on the "macroscopic" scale where the heterogeneous material can be considered as an equivalent homogeneous material. The aim of this work is then to study the influence of the microstructure of heterogeneous porous rocks on the behavior at the macroscopic scale.

The micro-macroscopic models used are enriched by microstructural data obtained from microscopic observations. Hence, the originality of this work is due to the unique combination of microstructural observations (SEM, X-ray 3D nano-tomography), multi-scale mechanical experiments (nano-indentation tests) and micro / macroscopic modeling. In addition, this work is innovative because it involves the study of natural materials (rocks) with realistic (heterogeneous) microstructures containing inhomogeneities (pores, cracks, inclusions) of different shapes, orientations and properties. Particular attention will be paid to oolitic porous rocks such as Lavoux limestone and iron ore, composed of an assembly of porous grains more or less spherical (oolites), pores and inter-granular crystals. These rocks were widely studied in GeoRessources laboratory as part of the feasibility study for the geological storage of carbon dioxide CO2 and the study of collapses of underground iron mines in Lorraine (France). Therefore, a detailed description of the mineralogical, petro-physical and hydromechanical properties of these rocks exists and it can be used as a reference in this work.

XXXIV

In Chapter I, we characterize the microstructure of the Lavoux limestone and the iron ore using three-dimensional imaging techniques such as scanning electron microscopy (SEM) and 3D computed nanotomography. X-Ray computed nanotomography allows distinguishing different components of considered rocks: porous network (size, spatial distribution and volume fraction), oolites and inter-oolite crystals. In addition, nano-indentation tests were performed to determine mechanical properties such as elastic parameters, which are necessary in homogenization models.

In Chapter II, we present statistical data describing several geometrical parameters (volume, radius, surface, sphericity, orientation…) of the components of considered rocks. We used a simplified model within the framework of Maxwell homogenization scheme where the porous oolites are approximated by spheres, and the pores of irregular shapes are approximated by ellipsoids. This approximation is performed using the Principal Component Analysis method (PCA), which provides the geometric properties such as the length of the semi-axes and the orientation of the resulting ellipsoids.

In Chapter III, we verify the approximation of irregularly shaped pores by ellipsoids by evaluating property contribution tensors of these pores. Thus, compliance contribution tensors for 3D irregular pores and their ellipsoidal approximations are calculated using the finite element method (FEM). These tensors were compared and a relative error is estimated to evaluate the accuracy of the approximation. The FEM numerical method is verified by comparing the numerical solution of compliance contribution tensors of the ellipsoids to the known analytical solution of these same shapes based on Eshelby's theory.

In Chapter IV, the numerical method used in the elastic problem is extended to thermal conductivity problem, where the approximation of irregularly shaped pores is also verified by evaluating thermal resistivity contribution tensors. Calculated compliance and resistivity contribution tensors were used to compute effective elastic and thermal properties of a material containing irregularly shaped pores by a two-step Maxwell homogenization scheme. Finally, computed properties of a material containing irregularly shaped pores are compared to those of the same material containing ellipsoidal approximations to evaluate once again the accuracy of the approximation of irregularly shaped pores by tri-axial ellipsoids.

To conclude, some essential results that have been done in this work will be reminded and a brief perspective of the future work will be mentioned

Chapter I

Microscopic observations and micromechanical properties of oolitic porous rocks SEM has been used to study the microstructure of a sample of Lavoux limestone. Several characteristic information was obtained:

(i) Topography or the surface features of the material.

(ii) Morphology i.e. the shape and the size of the particles composing the material.

(iii) Element composition of the material and the corresponding relative amount.

We used a JSM-7600F SEM (Figure 1.1) at GeoRessources Laboratory (University of Lorraine, Nancy, France), to visualize the sample of Lavoux limestone. This SEM combines two special technologies: an electronic column with a detection through the lens (semi-in-lens detectors) and a Schottky in-lens field effect gun for up to 400nA probe current. Thus, this SEM can achieve an ultra-high resolution with a wide range of beam current (1pA to more than 400nA). Furthermore, the JSM-7600F SEM has a real magnification of 10 6 with a resolution of 1 nm.

Equipped with a set of secondary electron and backscattered electron detectors, and an energy dispersive spectrometer, this SEM has the ability to present detailed data on the morphological characteristics and the chemical composition. Coupled with the SEM, we used an Energy Dispersive X-ray Spectroscopy (EDS) detector to perform qualitative and quantitative chemical analysis of the surfaces.

Mineralogical and morphological analysis of the Lavoux limestone using scanning electron microscopy (SEM) showed that it is made of oolites made up of concentric layers of microcalcite and bound by calcite cement, with a size that varies between 100 and 1000 micrometers (Figure 1.2). Three main types of calcite crystals can be found: large syntaxial crystals of several hundreds of micrometers, small equant calcites with an average size of 10 micrometers precipitated around the oolites in vadose conditions and microcalcites with an average size of 3 micrometers in the core of oolites.

The chemical characterization performed on Lavoux limestone sample by EDS showed that the material is mainly composed of calcite and contains a very small fraction of clays and dolomite.

Indeed, the nominal composition in atomic percentage is: O (75.3%), Mg (0.16%), AL (0.5%), Si (0.29%), Ca (23.42%), Fe (0.32%). Therefore, the Lavoux limestone is composed of more than 98% carbonates (calcite) and a small proportion of clays.

I.1.2 Description of Iron ore

We studied also the microstructure of another heterogeneous oolitic rock called iron ore from the Lorraine iron basin. The iron ore basin is over 100 km long and 50 km wide, from Luxembourg in the north to Nancy (France) in the south. The iron ore layers (iron content: 30% to 35%), 3-7m thick, are separated by intercalary rocks (argillites, siltites and calcarenite). This material has been widely studied by Grgic (2001;2005;2013) who described this material as an assemblage of grains of goethite: ovoid or spherical grains of 0.5 to 2 mm. This assembly of grains is organized in two different ways, depending on the size of the oolite. At their maximum size, oolites are organized in the form of an intersecting structure of variable size. The result is a strong heterogeneity and anisotropy of the ore of variable dimensions in space (Pineau, 1978). In general, oolites are cemented by an assembly of calcite called inter-oolitic cement, siderite and phyllosilicates in varying proportions. The pores of a few nanometers are filled with a largely crystallized spatonic carbonate cement, siderite limestone that generally provide good cohesion to the assembly (Grgic, 2001). The mineralogical composition of the interoolitic cement can vary from area to area, therefore iron ore can be considered as a heterogeneous rock from this point of view (Bubenicek, 1970).

Moreover, we observed the microstructure of the iron ore under SEM (Figure 1.3). SEM observations showed that the oolites of this material are made of concentric layers with an average size of 2 mm. In addition, a calcitic cement ensures the cohesion of the oolites and less inter-oolitic pores than the Lavoux limestone were observed.

An overview of the iron ore under the MEB coupled with the EDS showed that the cement is rich in ferrous phyllosilicates and carbonates. According to the chemical composition of the carbonates, it is found that the calcium carbonate is very pure (99% calcite). In contrast, iron carbonate is a mixture of siderite (predominantly 75%), calcite and magnesite. The hardness of a material is defined as the surface resistance presented by a material subjected to a surface penetration under the action of a given load. The choice of the indenter is made so that its deformation is neglected compared to the deformation of the tested material. After the contact of these two bodies, the indented material is affected by superficial modifications that depend not only on the intrinsic characteristics of the material but also on mode and conditions of application of the indenter. As a result, the hardness is in fact a complex property that cannot be directly related to other mechanical properties of the material.

In general, classical indentation tests are perfectly adapted to macroscopic measurement of hardness, by evaluating optically the residual imprint left by the indenter. However, an impression of a few microns in width is more difficult to measure optically. In addition, it is necessary to measure the hardness not post-test, but directly during the test, in order to determine the properties of the material under load (Tabor, 1970). Thus, it is for these two reasons that new machines were created 40 years ago. With these devices, the vertical force applied to the indenter is displaced continuously during the test. The displacement is then used to calculate the surface of the impression under load. This technique is called instrumented indentation. We used nanoindentation tests to determine the micromechanical properties of the three families of calcite of Lavoux limestone. This technique is developed in the 1970s and 1980s (Doerner and Nix, 1986;Loubet et al., 1984Loubet et al., , 1986;;Bulychev et al., 1975).

The field of application of the nanoindentation method is wide and covers a large range of materials. The method was applied to high temperature silicon by Suzuki and Ohmura (1996).

Then, it was used to measure the hardness of vitreous polymers at the nanometric scale by Hochstetter et al. (1999). In addition, the physical properties (modulus of elasticity and hardness) of molten silicon, calcium glass, aluminum crystals, tungsten and quartz were also calculated from the charge-discharge curves of the test of nanoindentation by Oliver and Pharr (1992).

Furthermore, several applications of this technique on geomaterials can be cited. The nanomechanical behavior of C-S-H (calciumsilicate hydrates) in cementitious materials was evaluated by Costantinides and Ulm (2007) thanks to a statistical analysis of hundreds of nanoindentation tests. In addition, nanoindentation tests have been used to present a systematic approach to treat interfaces between localized (fine grained) and peripheral domains at the atomic scale of crystalline solid (Karpov et al., 2006). Similarly, this technique has recently been used to determine the properties (elastic modulus, hardness) of different constituents of partially saturated argillite (Auvray et al., 2013(Auvray et al., , 2015(Auvray et al., , 2017)): a series of measurements was carried out by varying the relative humidity with a nanoindentation system placed in an air-conditioned chamber. This

I.2.2 Function and equipment:

The most known nanoindentation system that uses the technique of surface referencing is CSM nanoindentation system. Indeed, by controlling the relative position of the reference in real time with respect to the penetration depth of the indenter, the technique leads to several advantages: great accuracy on depth measurements, measurement in a short time and protection of examined surface against external influences e.g. acoustic disturbances.

After the nanoindentation test, it is possible to plot a force-penetration curve used to calculate the hardness of the material. The hardness is the ratio between the applied force F on the indenter and the projected contact area Ac. It is given by the following equation (Fischer-Cripps, 2002):

i c F H A  (1.1)
Since the hardness depends on the geometry of the indenter and the parameters of the test, then it is not an intrinsic parameter of the material. We used a Berkovich indenter with a triangular-shaped diamond pyramid (Figure 1.4).

) imp F H d   (1.2)
where F is the load in mN, imp d is the depth of the impression in nm and H is the hardness in GPa.

I.2.3 Mechanical properties measured by nanoindentation

I.2.3.1 Determination of elastic modulus

The first mechanical property that may be obtained from a typical nanoindentation test is the hardness. The second mechanical property obtained from this test is the elasticity modulus. To know the expression of this modulus, we have to study and discuss the Load-Displacement curve plotted after the nanoindentation test. During an indentation test, the acquisition system records the applied load as a function of the penetration depth of the indenter. These two parameters are continuously measured during the loading and the unloading phases. The result is a loaddisplacement curve that must be used later to deduce the properties of the tested material (elasticity, plasticity, viscoelasticity). An example of this type of curves is shown in Figure 1.5. The displacement hm of the indenter is measured normally with a capacitive sensor. Thus, we obtain the raw data of load and displacement. To obtain the value of penetration h, it is necessary to identify the point of contact, which we call in h . This curve has 2 distinct parts: The first one is AB which represents the loading phase and corresponds to the penetration of the indenter. Thus, the indenter reaches its maximal load Fmax at the penetration hmax that represents the maximal penetration of the indenter. This maximal state is represented on the curve by the point B. Loading curve can be represented mathematically using the following equation:

F Ah   (1.3)
where h is the penetration and  is the parameter of the power law describing the curve. The second part of the curve BC represents the unloading phase and corresponds to the withdrawal of the indenter. The intersection of the tangent to the unloading curve, and the abscissa axis is the plastic displacement hl. The slope of the tangent to the discharge curve represents the contact stiffness t S between the indenter and the sample. This calculation is possible thanks to the consideration that the contact between the indenter and the sample is purely elastic during the unloading (Bulychev et al., 1975). Hence, it is possible then to determine the reduced elasticity modulus of the contact between the indenter and the sample 𝐸 𝑟 from the unloading curve.

Considering that the beginning of the unloading curve is dominated by the elasticity of the material, one can use the equations of Sneddon calculated using Hankel transform, to evaluate the reduced elasticity modulus using the following equation:

2 t r c S E A   (1.4)
Sneddon calculated the expressions of the load and the displacements for conical, spherical and cylindrical indenters (Sneddon, 1965). Pharr showed that this expression can be used for any axisymmetrical indenter. For other type of indenters, the formula (1.4) is multiplied by a correction coefficient  , e.g. for Berkovich indenter  = 1.034 and for Vickers indenter,  = 1.012 (King, 1987;Oliver and Pharr, 2004):

2 t r c S E A    (1.5)
where t S is the contact stiffness.

Furthermore, one can deduce the elastic modulus of the material from the reduced elastic modulus of the contact between the indenter and the material since this modulus represents a combination of the contribution of the elasticity of both the material and the indenter. Indeed, the elastic modulus of the material can be obtained using the following equation:

1 𝐸 𝑟 = 1-𝜈 2 𝐸 + 1-𝜈 𝑖 2 𝐸 𝑖 (1.6)
where Er is the reduced elastic modulus, E and  are respectively the elasticity modulus and the Poisson's coefficient of the material, and i E and i  are respectively the elasticity modulus and the Poisson's coefficient of the indenter.

To summarize, 2 mechanical properties may be measured from the Load-Displacement curve: The hardness and the elastic modulus of the material. Thus, in order to calculate these two parameters, we need to define several parameters at the beginning:

• The load F.

• The elastic properties of the indenters ( i E and i



).

• The contact stiffness t S .

• The projected contact area c A The first three parameters are supposed to be known or may be easily measured. However, the evaluation of the contact area is not immediate. We will discuss in the following section the process of their calculation.

I.2.3.2 Determination of the contact area

This area is the parameter needed to calculate the elastic modulus of the material. This calculation is very delicate because usually we don't have access to this area but to the depth of the penetration left by the indenter. We will present the Oliver and Pharr approach which is the most used approach to describe the relation between the penetration of the indenter and the contact area: Oliver and Pharr described the depth of the penetration of the indenter h as the sum of the contact depth hc between the indenter and the sample and another depth hs obtained due to the elastic deflection considered outside the contact (Figure 1.5) (Oliver and Pharr, 1992;2004):

ℎ = ℎ 𝑐 + ℎ 𝑠 (1.7)
The depth h is measured during the test, so to determine the depth hs, Oliver and Pharr consider that the behavior of the material around the indenter is elastic. Therefore, one can refer to the equations of Sneddon to define hs for a conical indenter:

ℎ 𝑠 = 𝜋-2 𝜋 (ℎ -ℎ 𝑟 ) (1.8)
Moreover, from the expression of the force-penetration curve for a conical indenter indenting an elastic medium, we can relate h and r h to the load and the contact stiffness as follows:

ℎ -ℎ 𝑟 = 2 𝐹 𝑆 𝑡 (1.9) or ℎ 𝑐 = ℎ -𝜀 𝐹 𝑆 𝑡 (1.10)
where ε = 2 (π-2) / π = 0.727 is a constant obtained from the elastic theory and depending on the indenter geometry. To link the contact height hc to the projected contact area, Oliver and Pharr propose to use an area function to take into account the geometric imperfection at the end of the indenter (Oliver and Pharr, 1992;2004):

𝐴 𝑐 = 𝐴 0 ℎ 𝑐 2 + 𝐴 1 ℎ 𝑐 + 𝐴 2 ℎ 𝑐 1/2 … … + 𝐴 8 ℎ 𝑐 1/128 (1.11)
where 𝐴 0 represents the simple geometric relationship between c A and ℎ 𝑐 , 𝐴 1…8 are the coefficients to take into account the imperfection of the end of the indenter. In order to determine these coefficients, a nano-indentation test should be performed on a homogeneous material of which the elasticity modulus (silica in general) is exactly known. In other words, small indentation depths allow the determination of phase properties, while larger indentation depths lead to the response of the homogenized medium (Figure 1.7). The indentation process consists of applying an increasing normal on the surface of the sample using a diamond

Berkovich indenter with an elastic modulus of 1141 GPa and Poisson's ratio of 0.07 (Oliver and Pharr 1992). In addition, the indenter area function is determined as a function of the contact depth c h and is given by: 𝐴 𝑐 = 24.5ℎ 𝑐 2

(1.12)

Since our goal is to know the micromechanical properties of different families of the calcite of the Lavoux limestone, we performed three sets of tests:

• Set I of 400 indentations that cover the entire sample using grid indentation technique.

• Set II of 36 indentations focusing on the different layers of oolites.

• Set III of 6 indentations focusing on superficial surface of oolites, cement micro-calcite and cement macro-calcite.

Note that the first set was used to obtain a general overview on the elastic modulus of the sample, while the 2 others set were used to understand the behavior of each components.

One describes the general measurement protocol for the three series as follows:

(i) The number of indentation is fixed for each set.

(ii) The load is applied respecting the prescribed velocity of indentation, maximal load (loading stage).

(iii) The load is hold for certain time were the displacement and force are recorded continuously (iv) The indenter is removed (unloading stage).

(v) This procedure is performed in a repetitive manner at different fixed points on the sample surface.

(vi) The Force-Displacement curves are analyzed to calculate required parameters (Modulus of Elasticity).

Indentation parameters used in each test are given in Table 1.2. Force-Displacement curves were analyzed using the Olivar-Pharr method. The relaxation phase is used to dissipate and eliminate viscoelastic deformations. Therefore, equation (1.1) for the hardness and equation (1.6) for elasticity were used.

▪ Set I: where 𝜇 𝐽 is the arithmetic mean of all N values of each phase, and sJ, is the standard deviation that measures the dispersion of these values:

𝜇 𝐽 = 1 𝑁 𝐽 ∑ 𝐸 𝐾 𝑁 𝐽 𝑘=1
(1.14)

𝑠 𝐽 2 = 1 𝑁 𝐽 -1 ∑ (𝐸 𝐾 -𝜇 𝐽 ) 2 𝑁 𝐽 𝑘=1
(1.15)

Two cases can be considered:

•

The case of a single phase that corresponds to the case of homogeneous material, where the known mean and standard deviation are calculated to determine the fit distribution.

•

The case of n phases, where the overall frequency distribution of the elastic modulus obeys to the following theoretical probability density function:

𝑃(𝑥) = ∑ 𝑓 𝐽 𝑝 𝐽 (𝐸) 𝑛 𝐽=1
(1.16)

where fJ represents the volume fraction of each phase in the total volume. Therefore, the sum of all volume fractions should be equal to 1 i.e. representing 100% of the total volume. That reduce the number of unknowns to 3n-1.

To determine all these unknowns, one uses an optimization algorithm consisted of minimizing the standard error:

𝑚𝑖𝑛 ∑ (𝑝 𝑖 -𝑝(𝐸 𝑖 )) 2 𝑚 𝑚 𝑖=1
(1.17)

where p i is the experimental value of E, p(Ei) is obtained from equation (1.16) and m is the number of bins chosen to construct the histogram.

Figure 1.9 showed the presence of three families of calcite due to the presence to three peaks. This conclusion is supported by the previous SEM observations, and together allow a better understanding of the composition of the Lavoux limestone:

• Phase 1 corresponds to the first peak: Phase of inner layers of oolites (micro-calcite).

• Phase 2 corresponds to the second peak: Phase of outer layers of oolite (micro-calcite).

• Phase 3 corresponds to the third peak: Phase of inter-oolitic sparitic calcite (macrocalcite).

The calculated indentation moduli mean, standard deviation and volume fractions of each phase are given in Figure 1.9

▪

Set II:

We focused on oolite shown in Figure 1.10 thanks to a set of 36 indentation tests located on one line crossing an oolite. These tests cover different layers of this oolite. We used the same indentation protocol described before in order to obtain Load-Displacement curve and thus to calculate elastic modulus.

▪ Set III:

We focused in this set of indentation tests on the surface of oolites, cement micro-calcite and cement macro-calcite.

This set was used to verify the results of the first two set. Three conclusions can be made:

-Load-Displacement curves for the surface of oolites (Figure 1.11) give elastic moduli varying from 22 to 54 GPa, with a mean equal to 35 GPa. These high values verify the conclusion that these layers are the hardest layers of the oolites.

-The analysis on a cement micro-calcite gives an elastic modulus between 6 and 18 GPa.

These low values confirm that cement micro-calcite correspond to the first family of the three-phase distribution obtained from set I.

-The indentation on macro-calcite cement gives the highest elastic modulus values varying between 30 and 70 GPa. All these results are coherent with those of set I and verify the presence of the three main components in the Lavoux limestone sample. Over the last decade, the technique of X-Ray tomography has been used to analyze the microstructure of porous materials (Cloetens et al., 2002;1997;Maire et al., 2004;Farber et al., 2003;Mees et al., 2003;Taud et al., 2005). By definition, the tomography is a non-destructive method used to study the composition and the internal structure of an opaque object. It allows to characterize one or more physical parameters thanks to three dimensional images obtained by measuring the waves or radiations, emitted, transmitted or reflected by the examined object.

Several tomography methods are available such as ultrasound, magnetic, electrical, neutron, X-ray or gamma. The method used in this work and widely used in the medical and industrial field is Xray tomography (X-ray). Note that the tomographic data are digital images and then they provide access to many techniques of image analysis.

X-Ray tomography is a non-destructive technique that allows the reconstruction of the internal structure of a three-dimensional object thanks to cross-sectional images (virtual slices). These slices are reconstructed from the measurement of the attenuation of X-rays passing through the studied object. The X-ray attenuation measurement is mainly proportional to the density of the local electron. If the chemical composition of the object is uniform, the measurement is proportional to the local mass density of the object (Bossi et al., 1990). Thus, one can say that X-Ray tomography is a non-destructive measure of three-dimensional density fields. Moreover, X-Ray slices represent usually a finite thickness of the object and they are composed by voxels (voxel is the equivalent of a three-dimensional pixel). The size of the voxel corresponds to the spatial resolution of the measurement.

I.3.2 Data acquisition and processing

To characterize the microstructure of the Lavoux limestone and the iron ore, we used 3D X-Ray computed nanotomography. Tomographic data were obtained using a Phoenix nanotomogaph at the GeoRessources Laboratory (University of Lorraine, France). The specifications of the nanotomograph is given in Table 1.3. The scanned sample has a size of 10 mm of diameter for the Lavoux limestone and 4.67 mm for the iron ore. X-Ray images were recorded while the sample is rotating, step by step, over 360̊ along its vertical axis. The settings of the scan for the Lavoux limestone and the iron ore are given in table 1.4. These settings depend on the density of the considered material so they were chosen to guaranty first the best possible resolution, and second a good contrast between the components of the sample. Reconstruction of 3D volume data was done with system supplied DATOS_rec © software based on filtered back-projection algorithm developed by GE Measurement & Control Solutions (Wunstorf, Germany). Indeed, after the acquisition of the projections (X-ray images) which are the source images of the scanner, the first and the last images are compared by applying a displacement coefficient to correct the defects of the sample following its probable displacement during the scan. Then, the reconstruction zone is determined by an ROI (region of interest) and this zone depends on the part of the sample to be studied (Figure 1.12). The last step is to apply volume filters to the images in order to optimize their quality, so we applied a "Gauss Radius" filter and an adaptive contrast correction filter to make as large as possible the histogram of gray values that forms the images. The reconstruction is finally done using a computing cluster that exports a file with a ".vol" extension. The resulting isotropic voxel edge length or spatial resolution was 5.00 μm for the Lavoux limestone and 2.33 µm for the iron ore. This histogram has a range from 0 (Black) to 65535 (White) and 2 peaks: the first is representing the voids and the second the solid phase. However, the presence of two peaks does not mean necessarily the presence of 2 components only. That is expected since the Lavoux limestone is made from almost entirely of calcite (98%), therefore all the components other than the voids will have the same range of density or in other words same range of gray values in the histogram.

The same process was applied on the iron ore where two peaks were also observed on a total volume of 50 mm 3 .

The analysis of 2D and 3D slices showed that the Lavoux limestone and the iron ore consists of three main components (Figure 1.13; 1.14; 1.16): oolites having more or less spherical shape, interoolitic pores having irregular shapes, and interoolitic crystals (cement). It is shown in these figures the difference between the inter-oolitic porosity of the two materials where for the case of iron ore this porosity is lower and less represented in tomographic images. An important parameter that may be calculated from X-Ray images analysis is the porosity.

The determination of the porosity is based on the fact that the porous medium can be considered as the superposition of two parts: the solid part or white voxels, made up of various materials, and the empty part (pore space) or black voxels. Indeed, the porosity is the ratio between the number of black voxel and the total number of voxels. However, the visualization and the characterization process of 3D-images required high computational capacity. Therefore, in order to decrease the time of processing, one has to define a representative elementary volume from the total sample. By definition, the representative elementary volume (REV) of a heterogeneous material is an element of volume which is statistically representative of the properties of the material (composition, morphology, ...). In other words, the REV has to be big enough to represent the characteristics of the sample but as small as possible in comparison to the total volume (Kanit et al., 2003(Kanit et al., , 2006)). In our work, the fact to choose a REV is very useful in the calculation of the porosity and in the characterization of the components. This REV will be used for the rest of the work. In application to porosity, the REV is interpreted as a representative part of the media so that, if one selects such a volume in any region of the sample, porosity value should not undergo a significant variation. To determine that REV for the Lavoux limestone sample, we performed porosity calculation tests for several cubic subvolumes having the following side length: 1000 μm, 1200 μm, 1400 μm, increasing in steps of 200 μm till reaching 3000 μm (11 subvolumes). In addition, we repeated this procedure in five Therefore, a subvolume where the value of porosity is constant everywhere has to be selected. This condition is achieved from a subvolume having a side length of 2200-2400 μm, and the porosity tends to be constant at the average value of 8.33% for the Lavoux limestone, and a subvolume having a side length of 800 µm with a porosity equal to 3 % for the iron ore (Figure 1.19).

different
𝐹 𝐼 = 2𝜋𝑟𝛾𝑐𝑜𝑠𝜃 (1.18)
where 𝐹 𝐼 is the necessary force for expelling mercury, r is the pore entry radius, θ = 140° is the contact angle of the mercury on the material, and γ = 0.485 N.m -1 is the surface tension of the mercury.

𝐹 𝐸 = 𝜋𝑟 2 𝑃 (1.19)
where 𝐹 𝐸 the force needed to push the mercury into the pore under external pressure P.

Thus, at the equilibrium under pressure P, one writes the expression of the radius of pore filled by mercury from the balance of the two forces FI and FE:

𝑟 = -2𝛾cos (𝜃) 𝑃 (1.20)
The distribution of the porosity oh the Lavoux limestone as function of entrance radii of pores is given in Figure 1.21.a. This figure shows a spectrum containing a bimodal distribution of pore access size corresponding to two kinds of pores: intraoolite (0.001 < r <2 μm) and interoolite pores

(2 < r <30 μm), where r represents the entrance radius of pores. In addition, the first mercury injection curve gives access to the total porosity accessible with this method. For the Lavoux limestone sample tested here, the total mercury porosity (Nt, Hg) is 23.4%. However, before comparing the mercury porosity and tomographic porosity, several points should be clarified:

• Mercury porosimetry determines radius of access to the pores and not the radius of the pores.

Thus, we have an overestimation of the small radii and an underestimation of the larger accesses (Figure 1.20).

• Since we work with mercury under very high pressure, we are able to take into account very low radii, which is not possible with software like VGStudio especially when we have a resolution of 5 micrometers.

Figure 1. 20 : The notion of a radius of access to the pores (Bousquié, 1979).

To compare the two porosities for the Lavoux limestone, it is considered that we have access to only porosity of more than 5 micrometers (tomography spatial resolution). Mercury porosimetry gives a porosity of 7% for 5 micrometers access radius, which is on the same order to the porosity calculated by tomography.

For the case of the iron ore, the distribution of the entrance radii obtained from mercury porosity (Figure 1.21.b) showed that the intra-oolitic pores are dominating and few of inter-oolitic pores exist. This distribution explains the low value of porosity obtained by X-Ray tomography where the resolution of the sample was equal to 2.33. Same as in Lavoux limestone, porosity of tomography is underestimated since the intrapores cannot be represented in X-Ray images and then evaluated with an algorithm based on gray values. If we take only the average of the interpores obtained from mercury porosity, the value is equal to approximately 3% which verifies the porosity obtained from tomography. The importance of the study of the microstructure of considered materials is due to the effect of shapes of microstructural components of these materials on the homogenization schemes used later to evaluate effective elastic and thermal properties. In fact, the influence of the microstructure of heterogeneous porous rocks on macroscopic behaviour is taken into account by micromechanical models developed within the framework of the homogenization theory of random heterogeneous media. We will use Maxwell homogenization scheme detailed in chapter IV to evaluate effective properties in parallel to a simplified model that requires the approximation of the irregularly shaped pores by tri-axial ellipsoid and oolites by spheres. It is common practice in evaluating effective properties for three-dimensional pores to assume that the pores have ellipsoidal shapes. The main reason is that only ellipsoidal shapes have the property of uniform eigenstrain under remotely applied loading, so that the analytical solutions for strains and stresses around them can be utilized.

In this chapter, a statistical description of the components of the Lavoux limestone and the iron ore is presented, and a method to approximate these components by tri-axial ellipsoids is detailed within the framework of the simplified model.

II.1 Approximation of irregular shapes by tri-axial ellipsoids II.1.1 Introduction and theory

In general, a regular shape can be specified by geometrical parameters in relation with size, for example the diameter of a sphere, the aspect ratio and the length of one axis of a prolate or oblate spheroid in two dimensions or the semi-axis lengths of an ellipsoid in three dimensions. In other words, the quantification of such types of shapes is based on the analytical relation between their dimensions and their geometrical properties like volume and surface (Taylor et al., 2006).

The characterization of irregular shapes requires an important number of parameters. For example, Russ (1999) presented 10 parameters to quantify an irregular shape and Mather (1966) listed even more. The dimensions of 3D irregular shapes are usually defined using equivalent shape method.

To construct the equivalent regular shape, one can follow these steps: (i) one or more geometrical parameters are selected (volume, surface, inertia moments, etc.), (ii) these parameters are evaluated by some means for the irregular particle, (iii) a regular shape is selected, and its dimensions are determined by equating the geometrical properties of the irregular shape to the analytically known geometrical properties of the regular shape and solving the dimensions of the regular shape (Taylor et al., 2006).

Approximating 3D irregularly shaped pores by ellipsoidal shapes is very common practice in evaluating effective properties. Indeed, ellipsoids have the property of uniform eigenstrain under remotely applied loading, so that the analytical solutions for strain and stresses can be used (Eshelby, 1957;Mura, 2013). The approximation of irregularly shaped pores by ellipsoids can be based on several parameters depending on the application domain. In fact, two issues may arise:

the first is the choice of the best approximation of pores shape by ellipsoid (orientations and lengths of the principal axes) and the second is the accuracy of the selected approximation (Drach et al., 2011). In this work, we focus on mechanical applications that require the conservation of inertia moments of initial pores shape. That's the reason why we chose PCA method to approximate 3D irregularly shaped components of Lavoux limestone and iron ore by ellipsoids defined by their three semi-axes. PCA method can be also called singular value decomposition (SVD) (Golub and Van loan, 1983) or even eigenvalue decomposition (EVD). An eigenvalue decomposition is needed for a data covariance matrix after mean centering and normalizing the data matrix for each observation.

Alternatively, this method is based on moments of inertia and it provides the length and the orientation of the three principal axes of the ellipsoid. In order to apply this method in the case of irregular shapes, the following procedure was used:

1. 3D acquisition of real shape.

2.

Surface reconstruction and extraction of surface points.

3.

Ellipsoidal approximation based on PCA mathematical formulas.

II.1.2 Approximation of irregularly shaped constituents using principle component analysis (PCA)

Principal component analysis is probably the most popular multivariate statistical technique that can be used in different scientific domains. In fact, its origin can be traced back to Pearson (1901) or even Cauchy (1829), or Jordan (1874) and also Cayley, Silverster, and Hamilton (see Stewart, 1993;Boyer, 1989 for more details) but its modern instantiation was formalized by Hotelling (1933) who also coined the term principal component. Hence, the central idea of PCA is to reduce the dimensionality of a data set consisting of a large number of interrelated variables, while retaining as much as possible of the variation present in the data set (Joliffe, 2002).

Given an array containing p variables, this corresponds to p axes (p dimensions) under which we must represent the variables. The principle of the PCA is to project the data in a space of smaller dimensions allowing easy analysis of the data. It is therefore necessary to replace the p old axes with new ones. The new variables are expressed by the principal components Y1, Y2, …, Yp as a linear combination of the old variables.

Let ℂ be the covariance matrix associated with random vector X = [X1, X2, …., Xp] T . Given (λ1, e1), (λ2, e2), …, (λp, ep) as respectively the eigenvalues and the eigenvectors of ℂ, then the h th principal component writes:

Y h = e h T X = e h1 X 1 + e h2 X 2 … … + e hp X p (2.1)
where h = 1, 2, …, p and eh = [eh1, eh2, …, ehp].

eh1, eh2, …., ehp are the coefficients of these linear combinations that can be used to calculate the new points coordinates in the new coordinate system.

In order to obtain the principle components of a table data set, two important parameters were defined:

• The inertia of a column of the table: defined as the sum of the squared elements of this column:

𝛾 𝑗 2 = ∑ 𝑥 𝑖,𝑗 2 𝑛 𝑖 (2.2)
where the sum of all the 𝛾 𝑗 2 is called the total inertia.

•

The center of gravity (centroid, barycenter): defined as the mean of all the points in the column.

In general, the first principal component Y1 must have the largest possible variance or inertia.

However, this doesn't mean to choose large values for e11, e12, …, e1p because the choice of these coefficients should be restricted by the following relation:

e 11 2 + e 12 2 + ⋯ + e 1𝑝 2 = 1 (2.3)
The second principal component is calculated under the condition of being orthogonal to the first component and have obviously the second large inertia.

From a geometrical point of view, the new coordinate system or the new axes e1, e2, …, ep represented by these linear combinations are obtained thanks to an orthogonal transformation of the original system. In addition, given the p dimensional ellipsoid 𝑋 𝑇 ℂ -1 𝑋 = 𝑐 2 , the principal components define the axes of this ellipsoid. In our case, p = 3.

To find the three semi axes of the ellipsoid, first, we determine the mean of each variable, which are the surface points coordinates of pores in our case. Second, we subtract this mean from the values of the dataset to center the data around the origin. Third, we compute the covariance matrix ℂ of the surface points given as follows:

ℂ = (2.4)
This matrix is symmetric. By applying normal vectors decomposition, we rewrite this matrix as a function of eigenvectors and eigenvalues (Bronstein et al., 2012) as follows:

ℂ = ℚΛℚ (2.5)
where ℚ is the matrix of eigenvectors representing the direction of the three axes of the ellipsoid and Λ is the eigenvalues matrix having , as diagonal terms where 2  are the lengths of the three major semi-axes a, b and c of the ellipsoid.

II.1.3 Surface reconstruction algorithms

3D surface reconstruction refers to the technique used to obtain a three-dimensional representation of an object from a point cloud or a set of images taken from different points of view of the object. The reconstruction of a surface from a cloud of unorganized points can be stated as follows: at the input of the reconstruction process we have a set of unorganized points acquired from the surface of an object, and the aim is to produce at the outlet a surface which is as close as possible to the shape of the physical surface of the sampled object. In the literature, one lists several methods that lead to reconstruct 3D surfaces and allow to pass from a set of points with normals to surface, such as implicit methods (Hoppe et al., 1996;Curless and Levoy, 1996;Carr et al., 2001;Ohtake et al., 2003;2004), Delaunay methods (Boissonnat, 1984;Amenta et al., 1998;2000;2001;Dey and Goswami, 2003;Mederos et al., 2005), deformable methods (Kass et al., 1998;Whitaker, 1998;Zhao et al., 2001;Sharf et al., 2006;Osher and Sethian, 1988) Kazhdan et al. (2006). This method has several advantages:

•

Poisson reconstruction presents a global solution that considers all the data at once, and therefore no heuristic partitioning or blending will be faced.

•

Poisson reconstruction Poisson reconstruction creates very smooth surfaces that robustly approximate noisy data.

• Many methods based on implicit function, constrain this function only near the points forming the surface of the sample which contribute to a reconstruction that contains spurious surface sheets away from the sample. However, in Poisson reconstruction, the implicit function is constrained at all the points which lead to more realistic surface reconstruction.

Indicator function (Figure 2.1) is defined as follows: 1 at points inside the model and 0 at points outside. This function is useful thanks to its integral relationship with the oriented points sampled from the surface of the object. Moreover, since the indicator function is constant everywhere, its gradient is a vector field that is equal to zero also almost everywhere, except et points near the surface, where χ can be found by ensuring its gradient is as close as possible to the normal field N, in a least square sense:

‖∇χ -N ⃗⃗ ‖ (2.6)
We solve the Poisson problem for the function χ by applying the divergence operator in order to form the standard Poisson equation:

∇. ∇χ = Δχ = ∇. N ⃗⃗ (2.7)
Therefore, the surface is obtained via χ with a suitable iso-value, usually one takes the average value of χ at all input points. N ⃗⃗ is the smoothed normal field defined throughout the volume.

II.2 Approximation of oolites by spheres: shape study and statistical analysis II.2.1 Calculation of sphericity of oolites using selection tools

In this section, we will study the geometry of oolites which must be approximated by spheres if we use the simplified model (Giraud et al., 2015) within the framework of Maxwell homogenization scheme. To verify the approximation of oolites by spheres, we calculated their sphericity ratio. Thus, shape irregularity has a significant influence on this ratio. This irregularity in the case of the oolites affects not only the sphericity but it is manifested at three main scales (Krumbein and Sloss, 1963): Roundness R (cf. angularity), Sphericity Sph (cf. platiness) and smoothness SM (cf. roughness). In the literature, roundness is quantified as the ratio between the average radius of curvature of surface features and the radius of the maximum sphere that can be inscribed in the particle (Santamarina et al., 2004). Wentworth (1919) its analogue in three dimensions is called sphericity. Defined by Wadell (1935), the sphericity of a particle, is the ratio between the surface area of a sphere of the same volume as the particle and the actual surface area of the particle:

(2.8)

where Vp is the volume of the particle and Ap is the surface area of the particle.

Krumbein (1941) developed a fast method to estimate the sphericity based on a chart created by his own. Many authors like Sneed and Folk (1958) and Aschenbrenner (1956) proposed different

1/3 2/3 (6 ) p p V S A  
mathematical definitions for the sphericity ratio. The common point between all these formulas is their relation with the surface of the particle.

II.2.1.1 Influence of the composition of the Lavoux limestone on sphericity

The first method to calculate the sphericity of oolites is to extract them from X-Ray tomography images using an irregular shape selection tool in VGStudio MAX applied on different X-Ray 2D views of the material. Hence, one has to select the voxels that belong to one oolite and extract them as a new volume in order to calculate geometrical properties such as the volume, surface and dimensions. An example of 20 oolites (Figure 2.4) are extracted using this technique and sphericity was calculated using equation (2.8). For instance, a fast observation of oolites in MEB images or in tomography images of the REV shows that oolites are quasi-spherical. However, calculated sphericity leads to another conclusion. Figure (2.3) shows that sphericity varies between 0.58 and 0.79 which does not correspond to the real shape of oolites. The reason behind the underestimating of the sphericity is that the Lavoux limestone is a mono-mineral (calcite) geomaterial, which makes the identification of the true edges of oolites surface using selection method in tomographic images a difficult process since it is difficult to distinguish oolites from inter-oolitic cement (oolites and cement have the same gray scale values). Figure 2.2 shows the effect of this method on the surface of the oolite where some voxels are missed and the resultant shape does not represent accurately real oolite.

II.2.2 Calculation of sphericity of oolites using grinding process

To solve the problem raised in the previous section and to obtain more realistic values of sphericity, oolites should be liberated naturally from others components (i.e., calcite cement) and only in this case, real shape of oolites can be determined. To separate oolites from other components, we grinded a Lavoux limestone sample into different levels of grain size. Several grinding devices can be classified according to the type of stress they could apply in order to For particles larger than 1 mm in size, crushers (jaw crushers, gyratory crushers, etc.) and roller crushers are used, in which the fragmentation takes place by compression between two solid surfaces. They can be used for grinding hard to very hard materials. For the Lavoux limestone sample, we used two types of crushers:

-Jaw crusher (Figure 2.6; Left): It consists of a fixed jaw and a movable jaw reciprocated about a horizontal axis. The product to be treated is introduced into the upper part of the apparatus. When the moving jaw approaches the fixed jaw, it crushes the solid fragments.

When it spreads these downwards into a narrower part or they are again crushed and so on until they reach the exit orifice. We used this type of crusher to grind the large pieces of Lavoux limestone (Figure 2.7).

-Gyratory crusher (Figure 2.6; Right): A gyratory crusher is similar in basic concept to a jaw crusher, consisting of a concave surface and a conical head; both surfaces are typically lined with manganese steel surfaces. The inner cone has a slight circular movement, but does not rotate; the movement is generated by an eccentric arrangement. As for the jaw crusher, material travels downward between the two surfaces being progressively crushed until it is small enough to fall out through the gap between the two surfaces. Particles of each level were observed under a binocular microscope. For small levels of grinding (smaller than 200 μm), the shape of oolites was affected by crushing process so that they couldn't be used to evaluate the sphericity. However, complete oolites were found within the particles larger than 200 μm.

Moreover, the specific surface of the particle which defined as the total surface of an oolite depends on the smoothness factor of the surface which is a component of surface texture. Two independent geometric properties are the basics of surface texture: the degree of surface relief (rugosity) and the amount of surface area per unit of projected area. More the surface of the particle is rough, more it is overestimated and then the sphericity is underestimated. The surface from tomographic data is simply calculated voxel-wise, i.e. without a triangulation of the surface, so that the roughness will have a major influence on the values of specific surface, and then on the values of sphericity. Figure 2.8 illustrates the importance of the specific area in the calculation of sphericity ratio. The smoother this area, the smaller is the specific surface and then, the smaller is the error in sphericity ratio. Thus, to obtain a smooth surface for oolites, we proposed the following procedure:

1.

The oolites are extracted manually under the binocular microscope.

2.

The 40 selected oolites are put into a gel and a nanotomography scan is performed for the total sample (gel + oolites) (Figure 2.9a).

3.

Oolites are then approximated by ellipsoids using the PCA method (Figure 2.10) to obtain a smooth surface area. This surface area is used to estimate the sphericity. Sphericity ratio is recalculated using geometrical parameters of approximated ellipsoids. Figure 2.13a shows that adjusted sphericity values for ellipsoids approximating the oolites vary between 0.975 and 0.995. These values are more coherent with MEB and X-Ray images of oolites.

Furthermore, within the framework of proposed simplified model, these values are close to 1, and then, oolites can be reasonably approximated by spheres.

For iron ore material, the problem of density of materials is less affecting the selection of oolites from tomographic images, but the same procedure as the Lavoux limestone was applied because of its simplicity and accurate results. Oolites were then extracted in order to evaluate their sphericity ratio (Figure 2.9b). Figure 2.13b shows sphericity ratio of iron ore oolites. We observe that iron ore oolites have less spherical shapes than the Lavoux limestone ones. However, with values between 0.54 and 0.7 close to 1, one can still model oolites as spherical shapes in the simplified homogenization model.

Furthermore, a statistical analysis of oolites provided a detailed description on the following parameters: volume, radius, sphericity, center of every oolite, and surface. The radius of each oolite was calculated first by using the real volume of the oolite obtained from tomographic database and then by considering the oolite as a sphere so that:

(2.9) Thus, Figure 2.11 shows that Lavoux limestone oolites radii vary between 0.05 and 0.45 μm and they follow a logistic distribution with a mean μ = 0.218 and a standard deviation σ = 0.0348. Iron ore oolites radii (Figure 2.12) vary between 0.112 and 0.38 and they follow a normal distribution with a mean µ = 0.204 and σ = 0.053. Values of the mean of radius match with a good agreement with those found using MEB images. Therefore, we conclude that selected oolites are representative for all other oolites of the sample. -Diameter of the circumscribed sphere of the pore.

-Center X/Y/Z (voxel): Position of the center of the circumscribed sphere of the pore.

-Volume: total volume of all porosity and volume of each pore.

-Voxels: Number of voxels constituting the pore.

-Surface of each pore.

-Minimum gray value in the pore. the coordinate system. Figure 2.15 shows a possible anisotropic orientation due to the statistical small number of selected pores. However, SEM images shows an isotropic orientation and therefore, we considered that the pores inside the REV are randomly oriented and then the material is macroscopically isotropic. The same conclusion is made for the iron ore rock.

II.3.2 Extraction of 3D irregularly shaped pores from nanotomography data

The 3D pore structure inside the REV is used to select pores to be studied. Pores are selected randomly with one condition to not be too small due to meshing requirements. Therefore, a volume filter (V > 10^-4 mm 3 ) is applied. Figure 2.18 shows pore structure after applying the filter: For further evaluation of compliance contribution tensors, 10 pores were selected and numbered from 1 to 10 respectively. The geometrical properties of the pores determined using VGStudio 52 MAX are given in Table 2.1 for the Lavoux limestone and the same process is applied on iron ore sample; iron ore pores are presented in Appendix A. In order to apply the PCA method, these pores are exported as surface point cloud generated from X-Ray images using an extraction tool in VGStudio Max. This tool creates a triangular surface mesh every 10 µm and provides the coordinates of points of triangles constituting this mesh (Figure There are several parameters that affect the resolution and the quality of the reconstruction:

• Reconstruction depth: This integer is the maximum depth of the tree that will be used for surface reconstruction. The tree is a graphic parameter that helps the computer to organize the points of a 3D object very efficiently. Running at depth d corresponds to solving on a voxel grid whose resolution is no larger than 2^d x 2^d x 2^d. The first 2 parameters have a great influence on the generated surface:

•

• The higher the value for the octree-depth is chosen the more detailed results are obtained.

Thus, with noisy data, one keeps vertices in the generated mesh that are outliners but the algorithm doesn't detect them as such. So, a low value provides a smoothing surface but with less details. The higher the depth-value is, the higher is the resulting amount of vertices of the generated surface.

•

The samples per node parameter defines how many points the marching cubes algorithm puts into one node of the resulting octree. A high value like 10 means that the algorithm takes 10 points and puts them into the node of the octree. If we have noisy data, a high sample per node value provides a smoothing with loss of detail while a low value (between 1.0 and 5.0) keeps the detail level high. A high value reduces the resulting count of vertices while a low value remains them high.

So, a comparison of the effects of the depth value, the amount of samples per node and the effect of the offset parameter gives an idea on the best combination that should be used to obtain the best result. For the comparison of some parameter sets, we selected a random pore shown in Figure (2.19). This pore is presented from tomographic images as a point cloud. First, we studied the effect of the Depth value on the reconstruction of the surface (Figure 2.20). By reducing the depth value to 6, we obtained a model with less detail but a smoother surface. For a depth value equal to 8, we obtained a more detailed surface closer to the original one. From the value 8, no important changes affected the surface. Therefore, an acceptable level of detail and a reasonably smoothed surface were produced by using a depth value of 8. Second, by increasing the value of the samples per node with a constant depth value 6, it is observed that the surface gets even smoother (Figure 2.21). From these different figures, a value of 2-3 for the samples per nodes gives an acceptable result.

The best choice of these parameters is that which provides the closest pore shape to the original one. However, too many details on the surface of the pore make this pore unreadable by Comsol multiphysics software mesh module used later for a mesh process detailed in the following chapter.

For conclusion, the parameters set to reconstruct the surface using Poisson algorithm and make the pore shape are 8 for the reconstruction depth and 2 for the samples per node. These parameters ensure a good representation of the original shape and the availability of constructed shapes in further mesh process. 0.000377 0.00058 4.9 10 cov( , ) cov( , ) cov( , ) 0.000516 4.9 10 0.00124 The same process of Poisson surface reconstruction is performed on 10 randomly selected irregular pores. There is no important difference between different surface reconstruction parameters for the 10 pores. A reconstruction depth of 8 and a sample-per-node parameter between 1.5 and 2 was used was used for considered pores. Resultant reconstructed surface are presented in table 2.2.

x x x y x z y x y y y x z x z y z z                           0.
Table 2.2 shows also ellipsoidal approximation of the 10 selected pores usinf PCA method. For each pore, we obtained a correspondant tri-axial ellipsoid with the values of the three semi-axes.

We observe in this table that ellipsoidal approximations provide a good representation of all original irregular pores. Orientations of the three semi-axes of ellipsoids are the same of those of the original shapes. This was expected since the idea of PCA method is based on inertia moments and eigen vector decomposition. In addition, it seems that PCA is applicable in all the cases od considered shapes but it overestimated the volume of irregular pores. The influence of PCA on the volume of the pores and the influence of this volume on further mechanical applications is investigated in following section.

III.1 The first and second Eshelby problems The solution of the problem of inclusion with free deformation can be extended to the case where the medium M is subjected to a deformation 0 E or a stress 0  homogeneous at infinity. The deformation in the inclusion is therefore given as follows:

0

I E   𝕊 𝐸 : L I  (3.3) III.1.2
The second problem of Eshelby for ellipsoidal inhomogeneities III.1.2.1 Theory Consider a two-phase composite submitted to a remotely applied uniform loading at infinity and containing a region Ω having different elastic properties than the solid matrix. When the elastic moduli of this subdomain Ω differ from those of the remainder (matrix), the subdomain is called an inhomogeneity. Moreover, a material containing inhomogeneities is subjected to an internal stress (eigenstress) field, even if it is free from all external tractions. This problem with the presence of an inhomogeneity is called the second Eshelby problem or "inhomogeneity problem".

Our interest in this work is the elastic deformations on the boundary of the inhomogeneity due to the loading at the infinite of the matrix. Thus, the presence of this inhomogeneity will provide an extra strain or stress which is related to the said deformation by certain equations of integrals over V detailed below.

For ellipsoidal inhomogeneities (Ω), the first and the second Eshelby problem lead mathematically to the same solution, so one can use the analytical solution from the first Eshelby problem to verify the results of the second one.

However, for 3D-irregularly shaped inhomogeneities (real pores for example), the two Eshelby problems do not provide the same solutions and one has to carefully distinguish the difference between the two. Therefore, in the context of effective problems, second problem of Eshelby will be used and property contribution tensors characterizing contribution of individual inhomogeneity on the overall effective properties will be investigated.

III.1.2.2 Compliance contribution tensor of an inhomogeneity embedded in a matrix

Consider a REV with total volume V, bounded externally by surface While applying a uniform stress   on the boudary V  , the presence of an inhomogeneity disturbs the uniform stress and strain field, producing the variable stress field () x   and strain field () x   in M. Furthermore, the average strain is not only equal to   :

        (3.4)
where   is an extra strain due to the presence of inhomogeneities. Introducing n as the outward unit vector to the boundary ∂Ω of the inclusion and using the gradient formula of Gauss theorem, volume integral of strain tensor may be replaced by a surface integral and the average strain can be written as follows (Nemat-Nasser and Hori (1993); Sevostianov and Giraud (2013):

    11 2 T V u u dV V            11 . 2 n u u n dS V          (3.5)
The extra strain   due to inhomogeneity is then given by:

  11 . 2 n u u n dS V         (3.6)
or in component form:

  11 . 2 ij i i i i n u u n dS V       (3.7)
where u and n denote respectively displacements of the inhmogeneity boundary and a unit normal on  . V is the total volume (matrix + inhomogeneities) and , un nu denote a product of two vectors.

Furthermore, the relation between the extra strain is defined as function of the applied stress via a fourth rank tensor H called compliance contribution tensor:

* : ij ijkl kl V H V    (3.8)
where σ ∞ represents the applied stresses at the infinity, which is assumed to be uniform in the absence of the inhomogeneity.

III.2 Validation of the approximation of irregular pores by ellipsoids using compliance contribution tensors

To validate mechanically the approximation of pores by ellipsoids, we compared compliance contribution tensors of pores to those of approximated ellipsoids. Then, the results for the ellipsoids will be validated analytically (Eshelby solution for ellipsoidal pores) and an error between the real pores tensor and their approximations will be finally estimated.

Compliance contribution tensors for irregularly shaped pores need to be evaluated numerically.

Therefore, the elastic contribution tensor of irregular pores in the Lavoux limestone and their corresponding ellipsoids were determined using the finite element method (FEM).

We consider a homogeneous elastic matrix of volume V containing an inhomogeneity of volume V * composed by a different material. The compliance tensor of the matrix is 0 S while that of the inhomogeneity is 1 S . Since the calculation of compliance contribution tensors is independent of the Young's modulus of the matrix E0, we consider a value equal to 1 N/m 2 to simplify the calculations. However, this determination is not independent of the Poisson's ratio ν0, therefore Poisson's ratio value of Lavoux limestone should be taken into account. Thus, we took ν0 = 0.3 (Giraud et al., 2012). Indeed, to solve this elastic problem using finite element method, two parameters should be found:

• Displacement components denoted   , kl
 calculated at the nodes of the finite element mesh.

• Owtward unit vector n normal to the boundary of the inhomogeneity.

III.2.1 Numerical method for compliance contribution tensor of irregular pore and its ellipsoidal approximation

To evaluate the integral (3.15), 3D finite element simulation was performed thanks to Code-Aster, a free software for numerical simulation in structural mechanics:

III.2.1.1 Presentation of Code_Aster

Code_Aster is a general code for the study of the behavior of structures or materials by finite elements developed since 1989 by EDF (Electricity of France). From mechanical point of view, the behavior of the material leads to a link between the field of stresses and the field of deformations, directly according to Hooke's law in the simple case of an elastic behavior.

The main field of application of Code_Aster is the mechanics of deformable solids in static or dynamic form, but it also allows the modeling of hydraulic, thermal and acoustic phenomena whose effects can induce mechanical deformations. A phenomenon is a family of physical problems based on the same type of unknowns and associated with a type of conservation equation. Mechanical phenomenon, for example, uses the unknowns of displacement.

In addition to Code_Aster, 3 complementary software were used:

1-Gmsh (www.gmsh.info): is an open source three-dimensional finite element mesh generator with built-in pre-and post-processing facilities. Gmsh can produce a geometry and convert it to a usable mesh. Comsol Multiphysics (www.comsol.com): is a numerical software based on advanced numerical methods for modeling and simulating physics-based problems. We used this software to prepare the mesh, thanks to its advanced mesh functionality able to transform mesh surface into geometry element. Furthermore, the advantage of this software in meshing is in its automatic and semi-automatic meshing tools which create the mesh by defining a so-called mesh sequence. This mesh sequence allows to create a mix between several types of mesh elements (3D tetrahedral, 2D triangular or others). In addition, one controls the size of mesh elements to study the effect of the mesh size on the results. it can transform the mesh so that it can be readable by Code_Aster.

2-

III.2.1.2 Preparation of pore geometry and mesh

Finite element discretization of the pore surface was created after surface reconstruction.

Comsol Multiphysics software was used to generate the mesh of the pore. An example pore (volume = 0.00138 ) shown in Figure 3.3 is extracted from nanotomography images and then approximated by PCA method. The pore surface is provided by meshlab software as a mesh file. However, this mesh is formed by non-organized elements which cannot be read by Comsol workspace (Figure 3.3a). Therefore, a repair and remeshing of the surface is required. We used Gmsh to repair and remesh the surface to obtain an organized mesh for the pore that constitutes the input geometry element for Comsol (Figure 3.3b). Figure 3.4 illustrates the process of the mesh of the pore surface from tomographic image to Comsol geometry workspace. Then, irregular pore and its approximated ellipsoid are placed separately into a cube with limits ten times larger than the largest dimension of the pore. The exclusion of the irregular pore/ellipsoid from the cube is guaranteed thanks to "Difference" Boolean and partition operation.

II.2.1.3 Loading cases and command file in Code_Aster

Before running the numerical simulation, Code_Aster requires the writing of a command file that controls all the parameters of the simulation. For instance, a typical command file contains usually the following functions:

• Reading and modification of the mesh;

• Assignation of finite elements to the mesh;

• Definition of the properties of used materials;

• Assignation of the materials to the model;

• Set boundary conditions and loads;

• Choice of the adequate analysis type and solving;

• Calculation of the forces, stress, strain or more;

• Writing of the results in files, in ASCII and binary format.

Moreover, assigning the boundary conditions and materials to the model requires the naming of surfaces known in the code as group of elements. In our Cube-Pore model, adopted nomenclature and assigned materials properties are given by Figure 3.6.

From the Walpole matrix representation for Compliance contribution tensor and due to the symmetry, 21 independent components should be found. As a rule, to obtain all these 21 components, six independent loadings were considered: three uniaxial compressive stresses    (Table 3.2). For each case, 6 values were found. "a" is the side length of the reference volume cube, " " is the value of an applied strain and dx, dy, dz are the displacements in x, y, z directions.

Simulations are then run using Code_Aster, and displacement components are calculated at the nodes of the finite element mesh.

III.2.2 Numerical integration on the inhomogeneity surface (3D modelling) III.2.2.1 Integration on the surface of a quadratic 6-node triangular finite element Irregular pore is discretized with quadratic isoparametric 6-node Lagrange triangular finite elements (T6). The posed problem consists therefore of performing numerical integration on the curved surface of these finite elements. The surface integral 3.15 is then the sum of all the contributions of the integration on each element. We followed Zienkiewicz and Taylor (2000), Gnuchii (1978;1979) to calculate the integral on the surface of one T6 element. We define a reference plane on element surface given by 2 curvilinear coordinates ( , )

 . Thus, shape functions of the quadratic isoparametric 6-node Lagrange triangular finite element allow expressing global coordinates 1 2 3 ( , , ) z z z as functions of the two local (curvilinear) coordinates ( , )

 of the reference plane. Indeed, if z is the position vector on a T6 surface element, by using 0 Gauss Legendre quadrature rules, position vector coordinates i z may be expressed as function of local curvilinear coordinates ( , )

 : 6 11 1 ( , ) ( , ) i i i z N z        ; 6 22 1 ( , ) ( , ) i i i z N z        ; 6 33 1 ( , ) ( , ) i i i z N z        (3.12)
where () i j z denotes the cartesian coordinate j z of node () i .

By using basic tools of differential geometry, one defines metric tensor

ij g related to curvilinear coordinates 12 ( , );( , 1, 2) ij        : kk ij ij zz g     (3.13)
The integral on the whole surface is obtained by summing the contribution of all the finite elements. By denoting el N the total number of T6 finite elements on the pore boundary, the surface integral writes:

s zz zz N            ; 33 11 2 s zz zz N            ; 1 2 2 1 3 s z z z z N              (3.18)
The relation of this normal vector and the metric function ij g writes: Using the surface integration detailed above, and displacement components calculated from numerical simulations, equation (3.7) and (3.8) give the compliance contribution tensor for the 10 selected pores in Table 2.2. For instance, the results for pore shown in Figure 3.3 and its corresponding approximated ellipsoid are given respectively in the following tensors:

1

0 0 0 0 0 3.1471 0 0 0 0 0 0 2.7673 0 0 0 0 0 0 2.057                        (3.22)
The resulting matrices should be symmetric for mechanical considerations. This is the case for the approximated ellipsoid. For the irregular pore, numerical values of the non-diagonal components of each matrix obtained by processing FEA data were less than 1% of each other. Therefore, the matrices were symmetrized by calculating the average of each 2 corresponding non-diagonal terms:

2 ijkl ijlk ijkl ijlk HH HH   (3.23)
The terms of these two tensors (pore and ellipsoid) are very close, with a maximum relative error of 3%; the origin of this relative error is investigated below. Note that non-diagonal components of ellipsoidal tensors are equal to zero because it's the case of orthotropic symmetry which is not the case for irregular pores.

III.2.3 Analytical evaluation of compliance contribution tensor of approximated ellipsoids

In this section, we will present two methods to calculate the compliance contribution tensor using analytical method.

III.2.3.1 Analytical calculation of normal unit vector using implicit function

This first method is similar to the numerical method so that the same procedure is followed.

However, the normal unit vector is calculated thanks to an analytical solution. Hence, if the inhomogeneity i  is given by its implicit equation 1 2 3 ( , , ) 0 f z z z  , the outward unit normal vector to the boundary i  of the inhomogeneity is then proportional to the gradient of scalar function 1 2 3 ( , , ) f z z z . The formal proof of this result involves the theorem of implicit functions (Dini, 1878). Therefore, the outward unit normal vector is known on every point   1 2 3 ,, z z z belonging to the surface of the inhomogeneity. Thus, the components   , kl nn of this unit normal vector are known as functions of the spatial coordinates i z . They can be calculated at every point, including

Gauss integration points of a finite element, if Gauss-Legendre quadrature is used (for more details, one can refer to Barthélémy et al., 2016). As a result, compliance contribution tensor of the pore shown in Figure 3.3 using this analytical solution is given as follows:

1 2013) and we consider the case of an ellipsoidal inclusion phase embedded in a matrix. Hence, the relation between Eshelby and Hill tensors is given by the following equation within the framework of the first Eshelby problem:

ℙ = 𝕊 𝐸 : ℂ 0 -1 (3.25)
ℂ is the elastic stiffness tensor and index "0" refers to the matrix phase and index "1" refers to inclusion part. Thus, the importance of the first Eshelby problem is that one can determine the compliance contribution tensor as function of Hill and Eshelby tensors, i.e. for ellipsoidal inclusion the solution of compliance tensors may be obtained analytically. Indeed, compliance contribution tensor depends on the parameters of the ellipsoid but it is independent of the initial stress-strain state on the matrix. The relations between Hill's tensor and its counterpart for the compliance problem ℚ are given as follows:

ℙ = 𝕊 0 ∶ (𝕀 -ℚ ∶ 𝕊 0 ) (3.26) ℚ = ℂ 0 ∶ (𝕀 -ℙ ∶ ℂ 0 ) (3.27)
Thus, the compliance contribution tensor is given as follows:

ℍ = -𝕊 0 ∶ ℕ ∶ 𝕊 0 (3.28)
where N is a constant shape dependent fourth-order tensor given by:

ℕ -1 = (ℂ 1 -ℂ 0 ) -1 + ℙ (3.29)
For more details on the calculation of tensors ℕ, ℙ, and ℚ, one can refer to Appendix C.

• Numerical application

We applied this analytical method on the irregular pore shown in Figure 3 We used the Eshelby analytical solution for verification of the numerical method and a relative error indicating the distance between these two methods will be calculated for each pore.

III.2.4 Fobenius norm: relative error between numerical and analytical methods

To evaluate the relevancy of the ellipsoidal approximation, different calculated tensors should be compared. Thus, we chose to calculate the distance between these tensors thanks to

Frobenius norm defined as:

2 11 mn ij F ij Aa    (3.34)
Finally, to evaluate the accuracy of approximations, we defined Δ the distance between two matrices as:

exact approximation FF exact F AA A   (3.35)
We took analytical solution based on Eshelby problem as reference to calculate relative error because it is independent from the mesh. We applied presented procedure on 10 pores where they were extracted from the tomographic images of the chosen REV. These selected pores have different geometric parameters, different volumes and surfaces.

This method is equally extended to iron ore material (For more details, refer to Appendix C).

Table 3.3 shows the chosen irregularly shaped pores and their ellipsoidal approximations as well as the accuracy of this approximation for the Lavoux limestone pores. Table 3.3 illustrates also the distance between analytical and numerical matrices for these same pores. The distance between analytical and numerical solutions for iron ore material pores are presented in Appendix D. The values of the distance between ellipsoidal analytical (Eshelby) and numerical (Irregular pores) solutions for both rocks produce discrepancy of max 4.5 %. The distance between compliance contribution tensors of ellipsoids obtained via numerical and analytical solution based on implicit function, compared to Eshelby analytical solution, does not exceed 3% for both rocks. Even if this error reaches 4.5% in the case of ellipsoidal pore (Eshelby solution) in comparison to the original irregular pore, it is still in the acceptable range of difference from a mechanical point of view.

Analyzing these values, we note that there is no major difference between the highest and the lowest relative error for all the calculations. Therefore, proposed computational method is suitable for all the shapes of analyzed irregular pores. As a consequence, we can reasonably consider that the approximation of an irregularly shaped pore by an ellipsoid is relevant within the framework of the simplified model of Maxwell homogenization scheme. Since the numerical solutions are kind of approximation, the analysis of the quality of the mesh and its influence on the results is a fundamental point that deserves special attention. The geometry of the model used to evaluate compliance contribution tensors via 3D simulations is formed by the irregular/ellipsoidal pore and the cube. Thus, the pore zone has severe gradients of stresses and deformations. It is therefore crucial to choose carefully the optimal size of the discretization in order to obtain reliable results with a reasonable computation time.

III.3.1.1 Mesh generation

In our model, quadratic tetrahedral Lagrange finite elements were chosen for 3D simulations.

Theoretically, the mesh error must decrease with increasingly finer mesh, until it asymptotically reaches values independent of the mesh size. Thus, to analyze the quality of the mesh and its influence on the solution, 3 meshes have been tested for irregular and ellipsoidal pores. In all these cases, we generated mesh from Comsol where predefined element sizes levels were used.

Available mesh size levels vary from extremely coarse to extremely fine. We took into consideration the mesh around the pore to be finer than the mesh in the remote matrix. The number of nodes and elements forming each mesh are shown respectively in Tables 3.4:

• Mesh A (Figure 3.7a): Matrix: Normal; Pore surface: Fine.

• Mesh B (Figure 3.7b): Matrix: Finer; Pore surface: Extra Fine.

• Mesh C (Figure 3.7c): Matrix: Extra Fine; Pore surface: Extremely Fine.

Due to data storage limitations, no finer mesh could be generated. We calculated compliance contribution tensors for both irregularly shaped pore and its ellipsoidal approximation using Meshes A, B and C.

III.3.1.2 Discussion

The best mesh choice will be the mesh with less relative error in comparison to Eshelby analytical solution with a reasonable storage data and time computation:

• Case of irregular pore:

The Frobenius distance between compliance contribution tensors evaluated using Mesh A, B and C shows a maximum discrepancy of respectively 4.5%,4.9%,5.1% in comparison to Eshelby analytical solution. However, the computation of Mesh B and Mesh C took significantly more time and data storage.

This difference is understandable because the volume of the pore is small relatively in comparison to matrix size. Thus, by increasing the mesh size in the matrix, time calculation will increase but the effect of the matrix mesh on the results is negligible. Therefore, if we take into consideration the 2 parameters, precision and computing time, no important effect of the size mesh on H tensor components and Mesh A seems to be the most reasonable choice.

•

Case of ellipsoidal pore:

The Frobenius distance between ellipsoidal numerical solution of compliance contribution tensor from the analytical solution for Mesh A does not exceed 0.61%, Mesh B 0.44% and this distance decreases to 0.35% for Mesh C. As in the case of irregular pore, the time computation increased from 1 to 4 hours for one calculation and data storage increased also to hundreds of gigabytes.

Therefore, Mesh A was chosen as well for ellipsoidal case.

In conclusion, we choose to use quadratic elements with the Mesh A with finer mesh in the zone of the pore. This mesh size represents a good compromise between the accuracy of the results and the calculation time.

III.3.2 Effect of pore volume and surface area on compliance contribution tensor Porous materials are known to have combinations of several properties such as mechanical, thermal and many others. Effective properties of these materials depend on the intrinsic parameters of the pores. These parameters include for example the volume, the surface, the spatial arrangement, orientation and shapes of pores. The micro-scale geometric configuration of pores has an important influence on the macroscopic response. The influence of these parameters on effective properties were presented by many authors (Beyerlein, 2002). Indeed, when computing the compliance contribution tensors, the choice of the approximation method is controlled by the physical phenomena involved (mechanical, thermal, permeability, absorption, heat transfer etc ...).

Thus, each of these methods is based on different geometrical parameters to approximate the pores (Ordonez-Iranda and Alvarado-Gil, 2012;Kumar and Han, 2005).

III.3.2.1 Three-parameter equivalent shape models of irregularly shaped pores

For mechanical phenomenon presented before, inertia moments were the most important parameter to study. Thus, PCA method was used to provide the principle directions of the three semi-axes of the ellipsoids. In this section, we will investigate the influence of approximations methods on volume and surface of pores and then compliance contribution tensors of resulting ellipsoids will be compared.

In general, one-parameter equivalent shapes model is defined based on one geometrical property, two-parameter models by two geometrical properties and so forth. However, finding the new approximated ellipsoids, means finding 3 new semi-axes a, b, and c with the same orientation found using PCA method. This means that the one parameter surface or volume is not enough and one should use three-parameter equivalent models (Taylor et al., 2006). The volume of a tri-axial ellipsoid is:

4 3 ellipsoid V abc   (3.36)
This formula shows that the volume of the ellipsoid is simple in its composition where it is a product with a constant of the lengths of the semi-axes a,b and c. However, the surface area of an ellipsoid is more complex but Thomsen presented an approximated formula for the ellipsoid surface which accuracy is 1.061%  :

  or the surface of the shape. One can use it to first find the directions of the tri-axial ellipsoid and then adjust the volume since it is known from X-ray CT images. Formula of the volume of an ellipsoid makes it easy to define 3 parameters (or 3 equations) to find the new approximated ellipsoid that has the same volume of the original extracted irregular pore. This ellipsoid will have the same orientation of the one approximated by PCA but also the same semi-axes lengths ratios.

The three new semi-axes are denoted respectively ,, ab  and :

III.3.2.3 LWT method

One practical method to find the approximate tri-axial ellipsoid is to directly measure three dimensions denoted L for length (x direction), W for width (y direction) and T for thickness (z direction). A well known procedure consists to measure the longest line within the body and assign obtained value to the length L; then measure another longest line in y direction in condition to be orthogonal to L, and call it W; and finally a same procedure is applied to find T as the longest line orthogonal to L and W in z direction (Figure 3.9). Parameters for each of the 10 selected pores were measured using 2D sections in CT images provided by VGStudio MAX (Figure 3.9). Note lengths of semi-axes ,, abc    of tri-axial ellipsoid are the half of L, W and T. For pore in Figure 3.3, we obtained from 2D sections: L = 0.33, W = 0.13 and T = 0.08. r x y z    , and ij  is the kronecker delta function that is equal to unity if i = j and zero otherwise. The use of this inertia moment method and its relation with the covariance matrix and eigenvectors decomposition leads to PCA method that we have already presented and used.

Therefore, we had to define another order n of moment to find another realtions between the parameters of the tri-axial ellipsoid. For instance, one can choose xyz  or 4

x  but obtaining three unique values from them is not easy. Hence, another moment relations will be considered, defined by the three absolute first moments AFM as follows:

|| iV AFM x    (3.41)
These kind of moments presented in Press et al. (1989), computes the average deviation of the data from the mean. For the case of irregular pore, we take the center of inertia in the center of the shape, so in general if we calculate the average mean, this moment will be zero. However, here the moment is not zero since we take the absolute values of the coordinates. The Illustration of ellipsoidal approximation of the original irregular pore using AFM method (Yellow wireframe ellipsoid).

III.3.2.5 Discussion

Using the three choices of dimensions that have been defined, we can now define 3 threeparameter equivalent shape models for an irregular object: an ellipsoid with same volume of the original pore; en ellipsoid with semi-axes equal to 1/2L, 1/2W and 1/2T; and finally an ellipsoid defined from AFM values. The influence of these three-paramter models on the volume and surface will be considered. The slope factor is an indication about the precision of the correlation. Indeed, the closer this slope to unity is, the more useful and physical is the approximation of the tri-axial ellipsoid to the original shape. For volume conservation, it is easy to conclude that the method based on the same volume and same inertia moments as the original pore provides best results since the slope of this method is unity. However, this method affects largely the surface area where the slope is only equal to 0.76. PCA method is the second best method for volume conservation, with a slope equal to 1.19, and it presents a good surface area conservation with a slope of 0.86. Same for AFM method where the slope for the volume is 1.3 and for the surface area 0.91. LWT is a better solution than PCA and AFM if conservation of volume is needed, but it is less good for surface area.

Another important parameter that can be used to investigate the goodness of these correlation is the R 2 value. It gives an idea of how many data points fall within the results of the line formed by the regression equation. Applying this definition to our data, we conclude that for volume conservation, PCA+ Same volume method is confirmed as the best one with an R 2 = 1, then PCA method (R 2 = 0.9483), AFM (R 2 = 0.887) and finally LWT (R 2 = 0.875). Results of surface area, agree with the slope parameter where the same ranking of method is obtained: AFM (R 2 = 0.974), PCA (R 2 = 0.95), LWT (R 2 = 0.95) and finally PCA + Same volume (R 2 = 0.89).

To understand the influence of volume and surface area of pore on compliance contribution tensors, we evaluate these tensors for 10 selected ellipsoidal approximation pores. For ellipsoidal pores, we computed compliance tensor via numerical or analytical method. Then, new compliance contribution tensors were compared with those of original pores and a relative error was calculated.

For approximation methods based on volume conservation, we chose the same volume method, and for surface area conservation, AFM method was chosen.

Table 3.5 shows a maximum discrepancy of 4.5% for the influence of the volume on compliance contribution tensors and 5.1% for that of surface area. Therefore, proposed approximation methods for conservation of volume and surface area is suitable for all the shapes of analyzed irregular pores. In addition, relative error between used methods and PCA is relatively small and does not exceed 3%. As a consequence, we conclude that the use of PCA is enough for ellipsoidal approximation from a mechanical point of view. The homogenization method is used to determine the behavior of a heterogeneous medium at macroscopic scale thanks to the description of its constituents at the microscopic scale (Bornert et al., 2001;Zaoui, 1999). For heterogeneous materials, the homogenization method is considered applicable if these three characteristic lengths verify certain conditions:

• D -size of the heterogeneous structure.

• L -size of REV.

• d -size of heterogeneities.

The first condition is l << D so that the structure can be treated as a continuous medium. The second one is d << l, so that the behavior of the VER can be characterized by a homogeneous law (Figure 4.1). Then, to homogenize heterogeneous materials, the homogeneization methodology in three steps was proposed by Zaoui (1999):

-The stage of representation which allows to describe the material on the microscopic scale: definition of the phases and description of their spatial distribution and their mechanical behavior.

-The localization step that allows to model the mechanical behavior and to determine the relations between the quantities defined at the microscopic scale and the quantities of defined macroscopic equivalents.

-The homogenization step allows to analyze the previous results and put them in the expected global form to obtain macroscopic parameters.

Microstructure observations of the Lavoux limestone via MEB and X-Ray computed tomography,

showed that this rock is composed of quasi-spherical (oolites) grains cemented by a calcitic cement. The solid phase (micritic grains inside oolites and sparitic cement between oolites) is mainly made of pure calcite. Therefore, the reference elastic material used to study the solid phase is the pure calcite. This mineral belongs to trigonal anisotropic crystal system (Winkler, 2013).

However, as the distribution and the orientation of the mineral solid inclusions is random, the behavior and the shape of these inclusions have been simplified (Giraud et al., 2012) and isotropised (Ulm et al., 2005). So, we model the calcite "matrix" as the isotropic background. This simplification is used by Giraud and Sevostianov (2013), and it is relevant because of reasonably accurate experimental data on the macroscopic elastic constants of the oolitic rocks indicating its isotropic macroscopic behavior (Grgic, 2011;Nguyen et al., 2011).

Calculated contribution tensors H will be used as an input within the framework of the homogenization scheme to determine these effective properties. Three scales were identified (Figure 4.2):

-The microscopic scale that corresponds to the intra-oolitic pores -The mesoscopic scale that corresponds to the oolites, the inter-oolitic cement and pores.

-Third is the largest scale and corresponds to the representative elementary volume (REV)

which is considered large compared to the oolites size and inter-oolitic pores. The volume fraction of each phase is done as the ratio between the volume of each phase and the total volume  as follows:

o b c        (4.1) o o f    , b b f    , c c f    (4.2)
where 0 ,, bc f f f represent the volume fraction of each phase and -The resulting far-field is the sum of all the far-fields generated by all individual inhomogeneities within a homogenized region Ω.

-The resulting field at far points is equal to the far-field generated by the large homogenized region Ω with unknown effective properties. Moreover, the single pore micromechanical modeling is based on the compliance contribution tensors of the inhomogeneities calculated thanks to finite element analysis. These compliance contribution tensors describe the far-field generated by the inhomogeneities as well as their contribution into effective elastic properties. As a consequence, Sevostianov and Giraud (2013) showed that Maxwell homogenization scheme can be reformulated in terms of compliance contribution tensor to evaluate effective elastic properties. For this purpose, the same procedure as Maxwell was followed by considering a REV of volume V * from the composite and placed into the matrix material (Figure 4.4) and then, the resulting far-field is calculated thanks to same two methods described by Maxwell. A general solution of this problem has been recently presented by Sevostianov et al. (2005) for various inhomogeneities (including spherical inhomogeneity).

By equating the two quantities of the far-field, one obtains the general equation for the Maxwell scheme:

97 * 1 eff i i i V H V H VV   (4.7)
H-tensors of inhomogeneities are known and obtained by finite element analysis. However, the left hand of equation (4.7) is not known, it is described by the overall properties of the REV and its shape. The solution of this equation and the influence of the shape parameter is discussed by several authors (Sevostianov and Giraud, 2013;Sevostianov and Kachanov, 2012), where equation (4.7) is rewritten for ellipsoidal shape of inhomogeneities as follows:

1 1 0 * 1 eff i i reg i S S V H Q V               (4.8)
where reg Q is a tensor that reflects the shape of the REV. Equation (4.8) is valid for inhomogeneities of diverse shapes. The application of this homogenization scheme is described later to evaluate effective elastic properties of the Lavoux limestone containing irregularly shaped pores.

IV.1.2.2 Two-Step homogenization: Application to elasticity Using this model, the homogenization of the Lavoux limestone is carried out in two steps:

• First step: homogenization of micropores and solids grains inside the oolites:

The first homogenization step concerns intra-oolitic pores of spherical shape within oolites using self-consistent approximation. In this scheme, the inhomogeneity is embedded into material with effective properties that can be found from the solution of a single inclusion problem (Sevostianov and Giraud, 2013). This scheme was first developed for effective conductivity by Bruggeman (1935); by Kroner (1958) for elastic properties of polycrystals and by Skorohod (1961), Hill (1965) and Budiansky (1965) for effective elastic properties of matrix composites.

Self-consistent approximation employed in our micro-macro model allows the transition from the microscopic to the mesoscopic scale (Figure 4.5). In the framework of this first step homogenization, effective elastic properties of oolites at the meso scale, denoted o k (bulk modulus), o  (shear coefficient), are obtained considering spherical shape for solid grains and intra oolitic pores Compliance contribution tensors previously calculated by finite element method, are used to determine effective elastic properties of Lavoux limestone. Compliance contribution tensors are calculated considering an isolated pore inside a cube with respect to the reference system of the cube, so, the isotropy is not taken into consideration even if the material is macroscopically isotropic. Therefore, these tensors need to be replaced by corresponding isotropised tensor ℍ ̅ 𝑏

1 1 ( ) / oa ss o o o o o kf k k k k     , 1 1 ( ) / oa ss o o o o o f          (4.9) 3 34 o o oo k k     , 6 (
cc o c o c c c c o c c k H kk                 (4.18)
The reformulation of Maxwell homogenization scheme in terms of compliance contribution tensors gives the effective compliance tensor 𝕊 𝑒𝑓𝑓 as follows:

𝕊 𝑒𝑓𝑓 = 𝕊 𝑐 + {[𝑓 0 ℍ 0 + 𝑓 𝑏 ℍ ̅ 𝑏 ] -1 -ℚ 𝛺 } -1 = 1 3𝑘 𝑒𝑓𝑓 𝕁 + 1 2𝜇 𝑒𝑓𝑓 𝕂 (4.19)
where 𝕊 𝑐 is the compliance tensor of the matrix, and ℚ represents the counterpart of Hill tensor P for the compliance problem for domain  :

ℚ = ℂ 𝑐 ∶ (𝕀 -ℙ ∶ ℂ 𝑐 ) (4.20)
In the case of isotropic effective medium, the shape of the REV can be approximated by a sphere (Giraud and Sevostianov, 2013): and shear moduli for a selected irregular pore are compared to the solution of the corresponding numerical ellipsoidal approximation, where dashed lines correspond to the irregularly shaped pore and solid line corresponds to the ellipsoidal shape. Effective elastic parameters for all the other pores and ellipsoids were also compared. In all cases, the deviations of the effective elastic constants are small and do not exceed 5%. We conclude that for the most shapes, the PCA approximation of actual pores is convenient and its use in the context of Maxwell homogenization scheme is reasonable. Fourier law for macro level defined the effective conductivity tensor ** ij j K k i  as follows:

ℚ 𝛺 =
* . q K T       (4.22)
where T   and q  are the macroscopic temperature gradient and heat flux vector, respectively.

The evaluation of this conductivity tensor and its relation with the dipole moments were presented in Kushch and Sevostianov (2014). Recently, most of the results on conductivity of heterogeneous media are related to macroscopically isotropic materials. An application to thermal conductivity using a composite sphere assemblage model for porous rocks has been presented by Chen et al.

(2017). In addition, the thermal conductivity of porous media with spheroidal pores has been studied by Pabst and Gregorov (2014). The most known equations in the problem of heat conduction in materials write:

() .0 i q  (4.23) ( ) ( ) . ii i q K T    (4.24)
where () i T is the temperature, () i q is the heat flux and i K is the conductivity tensor of the i th  phase (i = 0 for matrix, and 1 for inhomogeneities).

By adopting notations similar to the elastic problem, one denotes  the temperature gradient vector T   (4.25)

In the same way, we introduce the second order resistivity contribution tensor r H related to an isolated inhomogeneity. Hence, one can assume a linear relation between temperature gradient T  and the heat flux q  (Fourier law).

T   and q  are taken as the volume-averaged values of the following local fields:

1 ;

V T Tdx V      1 V q qdx V     (4.26)
where V is the volume of the representative elementary volume of considered material.

106 Moreover, resistivity tensor may be calculated thanks to a numerical integration on the surface.

Hence, the gradient formula allows equally to rewrite the volume integral of thermal conductivity as a surface integral:

1 i i T T n dS xV           (4.27)
The second order resistivity contribution tensor r H is then defined as follows:

  * . R V T H q V     (4.28) Or in terms of components     * i R mn n m V T H q V     (4.29)
where q  denotes the uniform thermal flux vector, imposed on the infinite boundary, depending on the considered loading case. q  is also referred as the remotely applied thermal flux vector field that, in the absence of the inhomogeneity, would have been uniform.

IV.2.1.2 Loading cases, boundary conditions and simulations in Code_Aster

In general, resistivity tensor is characterized by 9 components where six of them are independent since the tensor is symmetric. To evaluate these components, the same procedure as the elastic problem is followed. However, only three loading cases are needed instead of six (Figure 4.12). In other words, three separate uniform thermal fluxes should be applied on the matrix cube: The three loading cases are described in table 4.1 as follows:

• Loadcase 1 (uniform thermal flux in x direction).

of the finite element mesh was obtained. Computation of resistivity tensor requires also numerical integration on the inhomogeneity surface and numerical calculation of normal unit vector to the surface of triangular finite elements. Hence, same procedure as mechanical problem was used to evaluate these two quantities via Gauss-Legendre quadrature. Therefore, one can use 2 methods for the case of ellipsoidal pores (analytical one by using implicit function to calculate the unit normal vector and numerical one) and 1 method for irregularly shaped pore (numerical method). Using the surface integration detailed before and temperature components calculated from numerical simulations, equation (4.30) and (4.31) give the resistivity contribution tensor for the selected pore shown in Figure 3.3: In the same way, the numerical result for the ellipsoidal approximation, where the normal unit vector is calculated thanks to vector product, is given as follows:

IV.2.1.3 Results and discussion

1 To evaluate the accuracy of the approximation of irregularly shaped pore with tri-axial ellipsoid, we calculated Frobenius distance between resistivity tensor of irregular pore and that of ellipsoidal one. In addition, to validate the numerical method, we evaluated the distance between two solutions (analytical: from implicit function and numerical: from vector product) for the same ellipsoidal pore. Table 4.2 shows the relative distance between different calculated tensors for the pore shown in Figure 3.3: Irregular pore 0

The same procedure was followed to calculate resistivity tensors for the 10 selected pores in Table 2.2. Two methods (numerical and analytical) are used to calculate these tensors in the case of ellipsoidal pore and the numerical method was used for the case of irregularly shaped pores. For resistivity tensors, we took the analytical method as a reference to verify the numerical ones since it is the least depending on the mesh. Table 4.3 illustrates also the distance between analytical and numerical resistivity tensors for the 10 selected pores: The extension of the evaluation of resistivity tensors for iron ore material is presented in Appendix C. Results for both materials show a maximum discrepancy of about 12 % between the numerical solution of irregularly shaped pores and numerical and analytical solutions of ellipsoidal pores. In addition, Frobenius distance between resistivity contribution tensors for ellipsoids obtained via numerical and analytical solution does not exceed 1.5 %. Even if these values are higher than the mechanical case, they still show a good agreement and an acceptable range of difference. Thus, analyzing the values of distance, we also conclude that proposed computational method is valid in thermal case for all considered shapes. As a consequence, we can once again reasonably consider that the approximation of an irregularly shaped pore by an ellipsoid is relevant within the framework of Maxwell homogenization scheme. where P  is the second-rank Hill's tensor for domain  .

Following homogenization scheme application, the reformulation of Maxwell homogenization scheme in terms of resistivity contribution tensor will be considered to evaluate effective conductivity of The Lavoux limestone.

IV.2.2.1 Two-step homogenization: Application to thermal conductivity Using the model presented above for the Lavoux limestone, the homogenization is carried out via two steps:

•

First step: homogenization of micropores and solid grains inside oolites:

This step represents the transition from the microscopic scale to the mesoscopic scale. Therefore, oolites are homogenized via self-consistent scheme and the result is the porous oolite. The selfconsistent scheme was chosen due to the random microstructure of oolites. One denotes respectively 0 ,, I oa    conductivity of the oolite at the mesoscale, conductivity of the intra oolitic porous phase, conductivity of micritic solid grains () oc   . We followed Torquato (2002) to define the self-consistent approximations for a two-phase material with spherical particles as the positive root of the quadratic equation: Values of thermal conductivity for irregular pores were compared to those for ellipsoidal ones and results for a selected pore are presented in Figure 4.16. Solid blue line corresponds to the irregularly shaped pore and dashed red line corresponds to the approximated ellipsoid. Effective conductive properties of all other studied pores were also compared. In all cases, very good agreement between irregular pores and ellipsoidal approximation was found and the deviations of the effective conductivity values do not exceed 4%.

As a result, we conclude that for considered shapes, PCA approximation of irregular pores is relevant and it can be reasonably used within the framework of Maxwell homogenization scheme to determine effective conductive properties. The approximation of inter-oolitic pores by tri-axial ellipsoids is verified by evaluating property contribution tensors of an isolated heterogeneity. Two problems were considered: Elastic problems where compliance contribution tensors were used and stationary linear diffusion problem where resistivity contribution tensors were used.

Due to high irregularity of pore, these tensors have to be evaluated numerically by finite elements method based on the second Eshelby problem. For the case of ellipsoidal pores, we presented three methods to calculate compliance contribution tensors: (i) Numerical method similar to the irregular pores case (ii) Analytical method based on implicit function of the ellipsoid (iii) Analytical method based on first Eshelby problem. For the resistivity contribution tensors only two methods were presented: numerical and analytical based on implicit function.

The relative error (Frobenius distance) between numerical compliance contribution tensors for irregularly shaped pores and analytical tensors for ellipsoidal approximations based on Eshelby problem, produces a deviation of the order of 4.5% in the elastic problem. This error between analytical solution based on implicit function and analytical solution based on Eshelby problem does not exceed 3%. The error between resistivity numerical contribution tensor for irregular pores and analytical solution based on implicit function was maximum 12%. For conclusion, from all these values of relative error, the approximation of the irregular pores by the ellipsoids is reasonable using the presented procedure.

The effective elastic properties (elastic modulus and shear coefficient) and thermal conductivity (thermal conductivity) were evaluated as a function of the porosity using contribution tensors calculated for the two cases: irregular pores and ellipsoids. Obtained results showed an important influence of the porosity on effective properties. In addition, a comparison between effective properties of a material containing irregularly shaped pores and the same material containing their ellipsoidal approximation showed a good agreement with a maximum deviation of 5%.

Among the perspectives, we mention the development of micromechanical models by introducing additional microstructural information obtained by micromechanical tests (nanoindentation for example).

Then, the increase of the number of selected pores and the consideration of pores with higher irregularity levels.

The current work is restricted to analysis of isotropic effective media. Extension to anisotropy could be also an important issue.

Numerical results presented for elastic problem and thermal conductivity for iron ore need to be compared with experimental data. Hence, it would be convenient to extend presented numerical method to other limiting cases, dry case and wet case for example, and to compare obtained results with experimental and predicted effects of porosity on the effective thermal conductivity.

Finally, the employment and comparison of other homogenization schemes and the study of their effects on effective properties would be interesting.

• Doerner, M.F., Nix, W.D. (1986). A method for interpreting the data from depth-sensing indentation instruments. J. Mater. Res. 1, 601-609. • Drach, B., Tsukrov, I., Gross, T., Dietrich, S., Weidenmann, K., Piat, R., Böhlke, T. (2011).

Numerical modeling of carbon/carbon composites with nanotextured matrix and 3D pores of irregular shapes. International Journal of Solids and Structures, 48(18), 2447-2457. • Dullien, F. A. L. (1981). Wood's metal porosimetry and its relation to mercury porosimetry. Powder technology, 29(1), 109-116. • Eshelby, J. D. (1957, August). The determination of the elastic field of an ellipsoidal inclusion, and related problems. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences (Vol. 241, No. 1226, pp. 376-396). The Royal Society. • Farber, L., Tardos, G., Michaels, J. N. (2003). Use of X-ray tomography to study the porosity and morphology of granules. Powder Technology, 132(1), 57-63. • Fischer-Cripps, T. (2002). Newnes interfacing companion: computers, transducers, instrumentation and signal processing. Newnes. • Giraud, A., Nguyen, N. B., Grgic, D. (2012). Effective poroelastic coefficients of isotropic oolitic rocks with micro and meso porosities. International Journal of Engineering Science, 58, 57-77. • Giraud, A., Nguyen, N. B., Grgic, D. (2012). Effective poroelastic coefficients of isotropic oolitic rocks with micro and meso porosities. International Journal of Engineering Science, 58, 57-77. • Giraud, A., Nguyen, N. B., Grgic, D. (2012). Effective poroelastic coefficients of isotropic oolitic rocks with micro and meso porosities. International Journal of Engineering Science, 58, 57-77. • Ohtake Y., Belyaev A., Alexa M., Turk G., Seidel H. P. (2003). Multi-level partition of unity implicits. In ACM SIGGRAPH, vol. 22, pp. 463-470. • Ohtake Y., Belyaev A., Seidel H. P. (2004). 3D scattered data approximation with adaptive compactly supported radial basis functions. In Shape Modeling International, pp. 153-164. (1992). Relationship of porosity and permeability to various parameters derived from mercury injection-capillary pressure curves for sandstone (1). AAPG bulletin, 76(2), 191-198. • Press, W.H., Flannery, B.P., Teukolsky, S.A., and Vetterling, W.T. (1989). Numerical Recipes:

The Art of Scientific Computing (Cambridge University Press, Cambridge). • Qiu, Y. P., Weng, G. J. (1990). On the application of Mori-Tanaka's theory involving transversely isotropic spheroidal inclusions. International Journal of Engineering Science, 28(11), 1121-1137. • Roussillon, T., Piégay, H., Sivignon, I., Tougne, L., Lavigne, F. (2009). Automatic computation of pebble roundness using digital imagery and discrete geometry. Computers Geosciences, 35(10), 1992-2000. • Russ J.C. (1999). The Image Processing Handbook, 3rd ed, CRC Pres, Boca Raton, Florida, p. 771.

• Santamarina, J. C., Cho, G. C. (2004, March). Soil behaviour: the role of particle shape.

In Advances in geotechnical engineering: The skempton conference (Vol. 1,. Thomas Telford, London.

• Sevostianov, I., & Giraud, A. (2012). On the compliance contribution tensor for a concave superspherical pore. International journal of fracture, 177(2), 199-206. • Sevostianov, I., Giraud, A. (2013). Generalization of Maxwell homogenization scheme for elastic material containing inhomogeneities of diverse shape. International Journal of Engineering Science, 64, 23-36. One may simplify this expression as follows:

kk ij i j i j i j i j x x x x x x x x g                             (B.

Figure 0 . 1 :

 01 Figure 0. 1: Les différents constituants du calcaire de Lavoux observés au MEB : oolithes sphériques avec des couches concentriques, pores inter-oolithiques de formes différentes et ciment inter-oolithique. ................................................................................................................................................................. XIX Figure 0. 2 : Les différents constituants du minerai de fer observés au MEB : oolithes avec des couches concentriques, pores inter-oolithiques de formes différentes et ciment inter-oolithique. ........................ XIX Figure 0. 3: Illustration en 3D des composants du calcaire de Lavoux obtenus par analyse des images tomographiques : (a) une oolite ; (b) ciment inter-oolitique. .................................................................... XX Figure 0. 4: Illustration en 3D des composants du minerai de fer obtenus par analyse des images tomographiques : (a) une oolite ; (b) ciment inter-oolitique ..................................................................... XX Figure 0. 5: Réseau de porosité d'un VER obtenu par analyse de données voxels tomographiques du calcaire de Lavoux (a) et du minerai de fer (b). .................................................................................. XXI Figure 0. 6: : Répartition de la fréquence expérimentale du module élastique montrant une distribution à trois phases avec 3 pics : pic 1 : couches intérieures des oolithes ; pic 2 : couches extérieures des oolithes ; pic 3 : macro-calcite interoolithique. ....................................................................................................... XXI Figure 0. 7: (a) Une oolithe obtenue à partir des images tomographiques ; (b) approximation de l'oolithe par une ellipsoïde par ACP. .................................................................................................................. XXIII Figure 0. 8: Le maillage de la surface et du volume du pore irrégulier et ellipsoïdal dans un cube de référence à l'aide des éléments quadratiques pour les simulations numériques élastiques et thermiques. .............................................................................................................................................................. XXIV Figure 0. 9: Modèle utilisé pour la deuxième étape d'homogénéisation basée sur le schéma d'homogénéisation de Maxwell. ........................................................................................................... XXVII Figure 0. 10 : Illustration du module de compressibilité élastique et du coefficient de cisaillement de 10 pores de forme irrégulière en fonction de la porosité en utilisant le schéma de Maxwell : chaque courbe correspond à un modèle contenant un des pores irréguliers orienté aléatoirement. ........................... XXVIII Figure 0. 11: Module de compressibilité élastique et coefficient de cisaillement (schéma de Maxwell) pour un pore sélectionné: vérification de l'approximation du pore de forme irrégulière par un ellipsoïde. .............................................................................................................................................................. XXIX Figure 0. 12: Gauche : Illustration de la conductivité thermique de 10 pores de forme irrégulière en fonction de la porosité en utilisant le schéma de Maxwell : chaque courbe correspond à un modèle contenant un des pores irréguliers orienté aléatoirement. Droite : La conductivité thermique pour un pore sélectionné: vérification de l'approximation du pore de forme irrégulière par un ellipsoïde. ................................... XXIX

Figure 1 . 9 :

 19 Figure 1. 9: Experimental frequency distribution of E showing the three-modal distribution with 3 peaks. .................................................................................................................................................................... 18 Figure 1. 10: : (a) Nanoindentation tests on oolite layers with corresponding Elastic moduli of each indentation point (b) Force-Displacement curves in: outer and inner layers. ............................................. 20 Figure 1. 11: Illustration of the indentation points of: Set II:indentation on inner and outer layers of the oolite; .......................................................................................................................................................... 21 Figure 1. 12: Region of interest selected for the volumetric reconstruction of X-Ray images. .............. 23 Figure 1. 13: Illustration of components of the the Lavoux limestone: (a) one example oolite; (b) Interoolitic cement .............................................................................................................................................. 25 Figure 1. 14: Illustration of components of the Iron ore: (a) one example oolite; (b) Inter-oolitic cement. .................................................................................................................................................................... 25 Figure 1. 15: Histogram evaluating pixel gray values of the sample ...................................................... 26 Figure 1. 16: Porosity network of selected REV of the Lavoux limestone and the Iron ore .................. 26 Figure 1. 17: Approximate location of the five subvolumes created for the determination of REV in section 2-D of the samples ...................................................................................................................... 27 Figure 1. 18: Study of porosity stability with variation of cube side length in different regions of Lavoux limestone sample. ........................................................................................................................................ 28 Figure 1. 19: Study of porosity stability with variation of cube side length in different regions of Iron ore sample. ........................................................................................................................................................ 28

Figure 0 . 1 :

 01 Figure 0. 1: Les différents constituants du calcaire de Lavoux observés au MEB : oolithes sphériques avec des couches concentriques, pores inter-oolithiques de formes différentes et ciment inter-oolithique.

Figure 0

 0 Figure 0. 2 : Les différents constituants du minerai de fer observés au MEB : oolithes avec des couches concentriques, pores inter-oolithiques de formes différentes et ciment inter-oolithique.

Figure 0 . 4 :

 04 Figure 0. 4: Illustration en 3D des composants du minerai de fer obtenus par analyse des images tomographiques : (a) une oolite ; (b) ciment inter-oolitique En utilisant la base de données d'un volume élémentaire représentatif (VER) fournie par nanotomographie et un algorithme de segmentation basé sur les niveaux de gris, la structure poreuse est obtenue (Figure 0.5). La porosité calculée (6.9% pour le calcaire de Lavoux et 3% pour le minerai de fer) représente le volume des pores interoolithiques (mésopores) dans le volume total.Une deuxième partie de pores qui est la microporosité n'est pas accessible par tomographie X car la résolution des images est limitée à 5 µm pour le calcaire de Lavoux et 2.5 µm pour le minerai de fer. L'analyse du VER de l'échantillon du calcaire de Lavoux et du minerai de fer fournit aussi une base de données statistiques concernant les composants du matériau (oolithes, pores interoolithiques, etc…). Plusieurs paramètres sont obtenus tel que le rayon, le volume, la distribution spatiale, la surface, ….

Figure 0 A

 0 Figure 0. 6: : Répartition de la fréquence expérimentale du module élastique montrant une distribution à trois phases avec 3 pics : pic 1 : couches intérieures des oolithes ; pic 2 : couches extérieures des oolithes ; pic 3 : macro-calcite interoolithique.

XXV•

  Méthode analytique basée sur le premier problème d'Eshelby Une autre solution analytique peut être obtenue pour le cas élastique seulement en considérant le cas d'une inclusion ellipsoïdale intégrée dans une matrice. Par conséquent, la relation entre les tenseurs d'Eshelby et de Hill est donnée par l'équation suivante dans le cadre du premier problème Eshelby: le tenseur de Hill, 𝕊 𝐸 est le tenseur d'Eshelby et ℂ est le tenseur de rigidité élastique.

(

  micrite) + micro-pores (b) Matrice poreuse : calcite (sparite) + macro-pores. Une homogénéisation à deux étapes est donc appliquée (Figure 0.9). La première étape d'homogénéisation concerne les pores intra-oolithiques de forme sphérique dans les oolites en utilisant un schéma auto-cohérant (

Figure 0

 0 Figure 0. 9: Modèle utilisé pour la deuxième étape d'homogénéisation basée sur le schéma d'homogénéisation de Maxwell.

  A l'échelle mésoscopique, nous considérons que le milieu hétérogène est formé par trois phases : les oolites poreuses (o) approximées par des sphères, des macro-pores interoolithiques (b) approximés par des ellipsoïdes et du ciment inter-oolithique (ou calcite syntaxiale) (c) constitué de grains de calcite pure. La Figure 0.10 montre les différents paramètres élastiques calculés en fonction de la porosité. Dans la Figure 0.11, les valeurs des modules de compressibilité élastique et de cisaillement pour un pore irrégulier sont comparées à la solution de l'approximation ellipsoïdale correspondante. Les propriétés élastiques effectives pour tous les autres pores et ellipsoïdes ont également été comparées. Dans tous les cas, les écarts des paramètres élastiques sont faibles et ne dépassent pas 5%.

Figure 0 .

 0 Figure 0. 10 : Illustration du module de compressibilité élastique et du coefficient de cisaillement de 10 pores de forme irrégulière en fonction de la porosité en utilisant le schéma de Maxwell : chaque courbe correspond à un modèle contenant un des pores irréguliers orienté aléatoirement.

Figure 0 .

 0 Figure 0. 12: Gauche : Illustration de la conductivité thermique de 10 pores de forme irrégulière en fonction de la porosité en utilisant le schéma de Maxwell : chaque courbe correspond à un modèle contenant un des pores irréguliers orienté aléatoirement.Droite : La conductivité thermique pour un pore sélectionné: vérification de l'approximation du pore de forme irrégulière par un ellipsoïde.

Figure 1

 1 Figure 1. 1: JSM-7600F SEM at GeoRessources Laboratory (Lorraine University, Nancy, France)

Figure 1

 1 Figure 1. 2: SEM (scanning electron microscope) observations on the components of the Lavoux limestone: spherical oolites with concentric layers, interoolitic pores of different shapes and inter-oolitic cement.

Figure 1

 1 Figure 1. 3: SEM observations on the components of the iron ore: oolites with concentric layers, interoolitic pores of different shapes and inter-oolitic cement.

  equipment allows testing under varying levels of saturation. Nanoindentation measurements controlled by temperature and humidity were used to determine the mechanical properties of the argillite phases. The nanoindentation tests were carried out under controlled temperature and humidity conditions where the viscoelastic properties were calculated by Vandamme and Ulm and the Fischer-Cripps models(Arnold et al., 2015).

Figure 1 . 4 :

 14 Figure 1. 4: Berkovich indenter shape; Berkovich imprint and residual deformations after a typical nanoindentation test on the Lavoux limestone.

Figure 1 .

 1 Figure 1.4 shows a Berkovich imprint and residual deformation after a typical nanoindentation test. The angle measured between an edge and the opposite face is 63.3 °. The hardness number is obtained as the ratio of the load and the impression depth:

Figure 1

 1 Figure 1. 5: Typical Load-Displacement curve.

I. 2 . 4

 24 Figure 1. 6: Schematics of the measurement head and composition of the CSM nano-indentation tester.

  Figure 1. 7: Grid indentation on a heterogeneous system where the probed microvolume (cross-hatched regions below the indenter) is either (a) smaller or (b) larger than the characteristic length scale, D, of the phase of interest (from chen et al., 2010)

Figure 1 .

 1 Figure 1.8 shows an example of a typical Force-Displacement curve obtained from the set I of the tested material.

Figure 1

 1 Figure 1. 8: Typical Force-Displacement curve from a nanoindentation test on Lavoux limestone

  Figure 1.10.b shows a typical Force-Displacement curve for an example of outer and inner layers of one oolite under the same experimental settings. The evolution of elastic modulus for different indentation points are given in Figure 1.10 for this oolite. Elastic modulus of different layers of oolite vary between 8 and 28 GPa. Figure 1.10.b shows that innerlayers are weaker than the outer ones and the maximal displacement hmax corresponding to the maximal load 40 GPa is equal to 2.29 μm. The maximal displacement for the same maximal load is 1.04 mm. Figure1.10 gives some indications on the hardness of each layer. This figure shows that inner layers have lower elastic modulus than outer ones. In fact, the porosity of layers has a major influence on their hardness. The more the layer is porous, the less the value of the elastic modulus and hence the less the hardness of the layer. Therefore, since the inner layers are more porous than the outer ones, they have lower elastic modulus and therefore they are weaker from hardness point of view.

Figure 1 .

 1 Figure 1. 10: : (a) Nanoindentation tests on oolite layers with corresponding Elastic moduli of each indentation point (b) Force-Displacement curves in: outer and inner layers.

Figure 1 .

 1 Figure 1. 11: Illustration of the indentation points of: Set II:indentation on inner and outer layers of the oolite; Set III: indentation on the surface of the oolite I.3 Computed X-Ray nanotomography:

Figure 1 .

 1 Figure 1. 13: Illustration of components of the the Lavoux limestone: (a) one example oolite; (b) Inter-oolitic cement

Figure 1 .

 1 Figure 1. 15: Histogram evaluating pixel gray values of the sample

  regions of the sample: upper left corner, upper right corner, lower left corner, lower right corner and center of the sample (Figure 1.17). The variation of the porosity for different subvolumes is given in Figure 1.18. By analyzing the graph in this figure, if a 1000 μm-side cubic volume is selected, we obtain 15.8% porosity in the upper right corner and 3% porosity in the lower left corner.This porosity represents the volume of inter-oolitic pores (meso and macro pores) over the total volume. Micropores inside oolites are not taken into account because they are not accessible due to the limitation of the resolution (5 micrometers) of X-ray tomography images.

Figure 1 .

 1 Figure 1. 17: Approximate location of the five subvolumes created for the determination of REV in section 2-D of the samples This means that for the same 1000 µm-side cubic volume, the difference in the values of porosity is significant. The reason behind this difference is that the choice of the subvolume may overestimate or underestimate the volume fractions of the constituents (pore and solid phases).

Figure 1 .

 1 Figure 1. 18: Study of porosity stability with variation of cube side length in different regions of Lavoux limestone sample.

Figure 1 .

 1 Figure 1. 19: Study of porosity stability with variation of cube side length in different regions of Iron ore sample.

Figure 1 .

 1 Figure 1. 21: Mercury porosity (1st injection): Distribution of entrance pore radii (porous spectrum for total porosity) for: (a) Lavoux limestone (b) iron ore(Grgic, 2011;2013)

Figure 2 . 1 :

 21 Figure 2. 1: Intuitive illustration of Poisson reconstruction in 2D (From Kazhdan, 2006).

  first defined roundness, and thenWadell (1932) presented the Wadell Roundness index that compares the outline of twodimensional projection of a particle to a circle in the same projection plane. For more details on the estimation of the roundness index, one may refer toCailleux (1947), Diepenbroek et al. (1992) andPissart et al. (1998) or toHayakawa and Oguchi (2005) andRoussillon et al. (2009) for computational methods of measuring true roundness. Hence, roundness applies in two dimensions;

Figure 2

 2 Figure 2. 2: 2D (a) and 3D (b) views of an oolite extracted from the Lavoux limestone sample using selection tool.

Figure 2

 2 Figure 2. 3: Values of sphericity ratio for 40 Lavoux limestone oolites extracted using selection tool based on gray scale values.

Figure 2

 2 Figure 2. 4: 3D view of selected oolites in the sample of the Lavoux limestone.

  Figure 2.5.

Figure 2 . 5 :

 25 Figure 2. 5: Classification of grinding processes according to the size of the obtained material.

Figure 2

 2 Figure 2. 6: Left: Jaw crusher; Right: Gyratory crusher.After grinding the material, we applied a sieving process in order to separate the particles according to their size. Three levels were obtained (Figure2.7): Particles larger than 200 μm; particles between 100 and 200 μm; and particles between 40 and 100 μm.

Figure 2 . 7 :

 27 Figure 2. 7: Different levels of grinfing of Lavoux limestone in order to separate oolites from other components.

Figure 2

 2 Figure 2. 8: Effect of the specific surface on sphericity ratio; S1 for smooth surface and S2 for real oolite surface.

Figure 2

 2 Figure 2. 9: Tomographic 3D images of oolites inside the gel after a grinding process; (a): Lavoux limestone; (b): Iron ore.

Figure 2 .

 2 Figure 2. 10: (a) (c) 3D view of an oolite respectively of Lavoux limestone and iron ore extracted after a grinding process from tomographic images (b) (d) Approximation of oolites by an ellipsoid using PCA to obtain a smooth surface for the calculation of sphericity ratio.

  Figure 2. 11: Distribution of 40 Lavoux limestone oolites radii obtained after grinding process; Logistic distribution with a mean = 0.218 mm and a std dev = 0.0348.

Figure 2 .

 2 Figure 2. 12: Distribution of 40 Iron ore oolites radii obtained after grinding process; Normal distribution with a mean = 0.112 mm and a std dev = 0.053.

Figure 2 .

 2 Figure 2. 13: Sphericity ratio of oolites after approximation by ellipsoids usining PCA method; (a): Lavoux limestone (b): Iron ore.

Figure 2 .

 2 Figure 2. 14: 3D views of the porous network of the REV of the Lavoux limestone obtained after applying a threshold algorithm Several geometrical and morphological parameters are obtained:

Figure 2 .

 2 Figure 2. 16: Distribution of pores radius obtained after a threshold algorithm based on gray scale values: Top: Lavoux limestone; Bottom: Iron ore.

Figure 2 .

 2 Figure 2. 17: Distribution of pores volume obtained after a threshold algorithm based on gray scale values: Top: Lavoux limestone; Bottom: Iron ore.

Figure 2 .

 2 Figure 2. 18: Porous network of the REV of the Lavoux limestone before and after the application of the volume filter.

  2.19a). 3D surface is then reconstructed using Poisson algorithm thanks to Meshlab software. The first step is to calculate normal of each point since the Poisson filter requires every point of the cloud to have an assigned vector. Therefore, this vector is the normal vector for a plane through that point. It describes the orientation and the direction of the plane (Figure 2.19b). Then, Poisson algorithm takes as input the set of 3D points with their oriented normals in order to reconstruct the surface. The output is a surface representation of the original shape (Figure 2.19c) exported into a variety of file formats such as .STL, .OBJ, .PLY, .3DS and can be imported into a 3D modeling program.

Figure 2 .

 2 Figure 2. 19: (a) Pore surface point cloud generated from X-Ray tomography (b) Computed normals on every point from point cloud (c) Reconstructed surface of the pore using Poisson algorithm.

  Samples per node: This floating point value specifies the minimum number of sample points that should fall within an octree node as the octree construction is adapted to sampling density. For noise-free samples, small values in the range [1.0 -5.0] can be used. For more noisy samples, larger values in the range [15.0 -20.0] may be needed to provide a smoother, noise-reduced, reconstruction. • Surface offsetting: This floating point value specifies a correction value for the isosurface threshold that is chosen. Values <1 means internal offsetting and values >1 means external offsetting.

Figure 2 .

 2 Figure 2. 20: Influence of reconstruction depth parameter on the reconstructed surface of a pore from the porous network of the Lavoux limestone.

Figure 2 .

 2 Figure 2. 21: Influence of sample per node parameter on the reconstructed surface of a pore from the porous network of the Lavoux limestone.

III. 1 . 1

 11 The first problem of Eshelby for ellipsoidal inclusionEshelby (1957) solved the problem of the mechanical equilibrium of an ellipsoidal inclusion, immersed in an infinite matrix, with the same mechanical properties as the matrix and subjected to free deformation. Let M be an infinite linear elastic medium, called matrix, of rigidity C, initially subjected to no deformation or stress. Let I be an ellipsoidal domain in M, called inclusion. The problem is illustrated in Figure3.1.

Figure 3 . 1 :

 31 Figure 3. 1: Illustration of inclusion problem.Eshelby showed that if a domain I of ellipsoidal shape undergoes a homogeneous free deformation L I  , the deformation

V

  . This volume contains an inhomogeneity  of volume V * , with the boundary  . The remainder of the REV, i.e. the volume excluded the inhomogeneity, is called the matrix denoted by M (Figure3.2).

Figure 3 . 2 :

 32 Figure 3. 2: Illustration of a volume V containing an inhomogeneity Ω.

  (www.salome-platform.org): is an open source software for pre-and post-Processing of numerical simulation. It is developped by CEA (French agency for nuclear power), EDF and OPEN CASCADE. This software contains several modules such as a mesh module, post-processing module and many more. It is often used with Code_Aster because

3 mm 68 Figure 3

 683 Figure 3. 3: Surface reconstruction and surface mesh repair of one selected irregularly shaped pore from the RVE: (a) (b) mesh surface before and after the repair; (c) Ellipsoidal approximation using PCA method of repaired pore. The resulting ellipsoid, obtained by PCA approximation, has the following semi-axes lengths: a = 0.1614 mm, b = 0.07 mm and c = 0.0467 mm. The volume of the approximated ellipsoid is:

Figure 3 . 4 :

 34 Figure 3. 4: Process of preparation of pore mesh surface to be imported to Comsol workspace for model mesh generation.

Figure 3

 3 Figure 3. 6: Reference volume, pore surface mesh, coordinate plane notation and material properties values.

  details on the integration on the surface of a quadratic 6-node triangular finite element using Gauss-Legendre quadrature rule, refer to Appendix B) III.2.2.2 Numerical calculation of normal unit vector to the surface of triangular finite element We recall equations (3.7) and (3.8) used to determine compliance contribution tensor. In equation (3.7), displacement components were calculated from numerical simulations, it only remains to compute the normals on the surface of quadratic finite elements to find the compliance contribution tensors.The position vector z on a T6 surface element writes: One calculates the normal vector to the surface defined by 2 parameters thanks to the vector product or also called cross product as follows(Zienkiewicz and Taylor, 2000

  ℚ and ℝ forming Hill tensor ℙ are calculated respectively:

Figure 3 . 7 :

 37 Figure 3. 7: Pore and ellipsoid volume mesh inside the reference cube using respectively Mesh A, B and C.

(

  , we consider the selected irregularly shaped pore shown in Figure33 semi-axes: a = 0.1614; b = 0.07; c = 0.0467 mm. By solving the system of 3 equations with 3 unknowns, the new semi-axes of the ellipsoid method on volume and surface area is presenetd below.

Figure 3

 3 Figure 3. 8: Left: Illustration of ellipsoidal approximation of irregularly shaped pore by PCA (Red ellipsoid); Right: Illustration of ellipsoidal approximation having the same volume of the original irregular pore (Green ellipsoid).

Figure 3

 3 Figure 3. 9: Left: Measurement of L, W and T values from 2D sections in VGStudio MAX of an irregular pore; Right: Ellipsoidal approximation of the irregular pore using LWT method (Magenta ellipsoid).III.3.2.4 Absolute first moment method (AFM)This is another method based on moments. In general, all bodies possess a center of mass and a center of volume. If the density of the body is homogeneous, then, these 2 centers are the same. Therefore, the center of volume of a body writes:

Figure 3

 3 Figure 3. 10: Left: Illustration of ellipsoidal approximation (Red wireframe ellipsoid) of irregularly shaped pore by PCA; Right:

  Figure 3.11 and 3.12 contain 8 curves (4 volume and 4 surfaces) of approximate volume and approximate surface area obtained from 4 methods vs. volume and surface area of the original pore. The 4 methods used are: PCA, PCA + Same volume, LWT and AFM. Figure 3.11 and 3.12 show an acceptable linear corralation between the various approximations and original values. To investigate the goodness of these correlations, one calculates the slope of curves and the R 2 values.

Figure 3 .

 3 Figure 3. 11: Influence of approximation methods on volume of irregular pore: linear correlation between approximation values and volumes obtained from CT images.

Figure 3 .

 3 Figure 3. 12: Influence of approximation methods on surface of irregular pore: linear correlation between approximation values and surfaces obtained from CT images.

Figure 4 . 1 :

 41 Figure 4. 1: Three characteristic length scales involved in homogenization scheme

Figure 4

 4 Figure 4. 2: Three scales on the level of the material from left to right: Macroscopic scale: scale of REV, Mesoscopic scale: scale of oolites and interoolitic pores, Microscopic scale: porous oolite, scale of intraoolite pores. Furthermore, the presence of pores in the oolites and in the cement (matrix) allows us to propose a geometric model for the homogenization of oolitic Lavoux limestone as follows:

Figure 4

 4 Figure 4. 3 :Proposed model for the second homogenization step based on Maxwell's homogenization scheme.

Figure 4 . 4 :

 44 Figure 4. 4: Scheme of Maxwell's homogenization method. Effective properties of a composite (a) are calculated by equating effects produced by a set of inhomogeneities embedded in the matrix material (b) and by fictitious domain having yet unknown effective properties (c).



  denote bulk modulus and shear coefficient of solid phase constituting oolite o.

Figure 4 . 5 :

 45 Figure 4. 5: First homogenization step: micropores inside oolite core are homogenized by using self-consistent method (2D representation of a 3D microstructure).

Figure 4 . 6 :

 46 Figure 4.6: Dependencies of bulk modulus of pores of irregular shape as a function of the porosity using the Maxwell scheme: each curve corresponds to a model containing one of the irregular pores randomly oriented.

Figure 4 . 6 :

 46 Figure 4. 6: Dependencies of shear coefficient of pores of irregular shape as a function of the porosity using the Maxwell scheme: each curve corresponds to a model containing one of the irregular pores randomly oriented.

Figure 4

 4 Figure 4. 7: Dependencies of bulk modulus of approximated ellipsoidal pores as a function of the porosity using the Maxwell scheme: each curve corresponds to a model containing one of the irregular pores randomly oriented.

Figure 4 . 8 :

 48 Figure 4. 8: Dependencies of shear coefficient of approximated ellipsoidal pores as a function of the porosity using the Maxwell scheme: each curve corresponds to a model containing one of the irregular pores randomly oriented.

Figure 4

 4 Figure 4. 9: Bulk modulus (Maxwell scheme) for a selected pore: verification of the approximation of the irregularly shaped pore by an ellipsoid.

Figure 4 .

 4 Figure 4. 10: Shear coefficient (Maxwell scheme) for a selected pore: verification of the approximation of the irregularly shaped pore by an ellipsoid.

  schemes, conductivity is defined by the resistivity tensor, the counterpart of compliance contribution tensor in thermal problems. This tensor denoted ℝ was introduced by Sevostianov andKachanov (2002) in the context of cross-property connection between elastic and conductive properties of heterogeneous materials.

Figure 4 .

 4 Figure 4.13 shows the field of temperature obtained from loading case 1, where the temperature is given in Kelvin. This figure illustrates the thermal flux generated inside the reference cube due to the application of two opposite temperatures on both south and north sides.

Figure 4 .

 4 Figure 4. 12: Temperature field result from loading case 1 in Kelvin.

IV. 2 . 2

 22 Maxwell homogenization scheme reformulated in terms of resistivity contribution tensors Another application of Maxwell homogenization scheme is the evaluation of the effective conductivity of the Lavoux limestone. This method has been used in the context of effective conductivity by Kushch et al. (2014); Sevostianov et al. (2014). Maxwell's homogenization model has been used inZimmerman (1989) to predict the thermal conductivity in a wide range of rockmicrostructure type of sedimentary by considering randomly distributed spheroidal pores. In fact, the same idea as the mechanical problem is employed to solve this problem. The far-field generated by inhomogeneities is calculated via two methods and by equating the results, effective conductivity is then obtained. For the thermal conductivity problem, the result has the

•

  Second step: transition from the mesoscopic scale to the macroscopic scale with ellipsoidal pores: At the mesoscale, we considered three-phase heterogeneous medium, which is composed of porous oolites (volume fraction 0.7 o f  from statistical study in chapter II), pure solid calcite referred as spiritic cement, and irregularly shaped mesopores randomly distributed and approximated by ellipsoids. To obtain conductive effective properties of the Lavoux limestone, we followed the reformulation of Maxwell homogenization scheme in terms of resistivity contribution tensors by Sevostianov and Giraud (2012) described above. This reformulation provides the simple scalar relation established under assumption of macro isotropy: tensor b H of irregular pore is calculated thanks to FEM. A random orientation distribution of irregular pores conduces to replace b H by its isotropized tensor b on the influence of the porosity on overall thermal conductivity is considered. Figure 4.14 illustrates the dependencies of the thermal conductivity on the volume fraction b f for a material containing irregularly shaped pores with random orientation and Figure 4.15 illustrates the dependencies of the thermal conductivity on the volume fraction considering the same material containing 10 ellipsoidal approximations of original irregular pores. Each curve corresponds to one model of the 10 studied pores. For a volume fraction of meso porosity equal to 0 i.e. in the absence of any meso porosity the thermal conductivity is equal to 2.62. The values of thermal conductivity for both cases (irregular and ellipsoidal pores) decrease with the increase of the porosity fb since the thermal conductivity of pores is much lower than that of the solid. As in the case of elastic effective properties, the values of thermal conductivity for the 10 selected pores are relatively close (less than 3% from the mean values).

Figure 4 .

 4 Figure 4. 13: Effective thermal conductivity as a function of mesoporosity: each curve corresponds to one of the 10 irregular pores with random orientation.

Figure 4 .

 4 Figure 4. 14: Effective thermal conductivity as a function of mesoporosity: each curve corresponds to one of the 10 irregular pores with random orientation.

Figure 4 .

 4 Figure 4. 15: Effective thermal conductivity for a selected pore: verification of the approximation of one irregularly shaped pore by an ellipsoid.

•

  Giraud, A., Sevostianov, I. (2013). Micromechanical modeling of the effective elastic properties of oolitic limestone. International Journal of RockMechanics and Mining Sciences, 62, 23-27. • Giraud, A., Sevostianov, I., Chen, F., Grgic, D. (2015). Effective thermal conductivity of oolitic rocks using the Maxwell homogenization method. International Journal of RockMechanics and Mining Sciences, 80, 379-387. • Gnuchii, Y. B. (1978). Numerical integration over a surface in the finite-element method. Strength of Materials, 10(6), 704-706. • Gnuchii, Y. B., Tsybenko, A. S. (1979). Numerical integration over the surface of plane and axisymmetric finite elements. Strength of Materials, 11(6), 649-651. • Golub G.H., Van Loan C.F. (1983). Matrix Computations, North Oxford Academic. • Grgic, D. (2001). Modélisation du comportement à court et à long terme des roches de la formation ferrifère lorraine, PhD Thesis, I.N.P.L. -Nancy -France. • Grgic, D. (2011). Influence of CO2 on the long-term chemomechanical behavior of an oolitic limestone. Journal of Geophysical Research: Solid Earth, 116(B7). • Grgic, D., Giot, R., Homand, F., Giraud, A. (2005). Effect of suction on the mechanical behaviour of iron ore rock. International journal for numerical and analytical methods in geomechanics, 29(8), 789-827. • Grgic, D., Giraud, A., Auvray, C. (2013). Impact of chemical weathering on micro/macromechanical properties of oolithic iron ore. International Journal of Rock Mechanics and Mining Sciences, 64, 236-245. • Hayakawa, Y., Oguchi, T. (2005). Evaluation of gravel sphericity and roundness based on surfacearea measurement with a laser scanner. Computers & geosciences, 31(6), 735-741.

•

  Oliver, W. C., Pharr, G. M. (1992). An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments. Journal of materials research, 7(06), 1564-1583. • Oliver, W. C., Pharr, G. M. (2004). Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology. Journal of materials research, 19(01), 3-20. • Ordonez-Miranda, J., Alvarado-Gil, J. J. (2012). Effect of the pore shape on the thermal conductivity of porous media. Journal of Materials Science, 47(18), 6733-6740. • Osher, S., Sethian, J. A. (1988). Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. Journal of computational physics, 79(1), 12-49. • Pabst, W., & Gregorová, E. (2014). Conductivity of porous materials with spheroidal pores. Journal of the European Ceramic Society, 34(11), 2757-2766. • Pearson, K. (1901). LIII. On lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11), 559-572. • Pineau, J.L. (1978). Contribution à la caractérisation géomécanique des roches. Application au minerai de fer lorrain. PhD thesis, I.N.P.L -Nancy -France. • Pissart A, Duchesne F, Vanbrabant C. (1998). Determination of confidence limits characterizing 563 gravel petrographic counts and roundness values. Comparison of roundness values by the Cailleux and the 564 Krumbein methods. In: Geomorphologie: relief, processus, environnement, 195-214. • Pittman, E. D.

  derivatives may be found in standard finite element textbooks (see among othersZienkiewicz and Taylor (2000)). One defines metric tensor ij g

  of this tensor were used to evaluate the contribution of one quadratic triangular element T6 as follows:

  the whole surface obtained by summing the contribution of all the finite elements.By denotingel N the total number of T6 finite elements on the pore boundary, the surface

Table 1 .

 1 1: Technical specifications of nanoidenter tester ...................................................................... 15 Table 1. 2: Indentation series and parameters. ........................................................................................ 16 Table 1. 3: Technical specifications of the nano-tomograph. ................................................................. 22 Table 1. 4: Tomographic scan settings. .................................................................................................. 23 Table 2. 1: Properties of 10 selected pores from the Lavoux limestone porous network. ...................... 52 Table 2. 2: Surface reconstruction and ellipsoidal approximation for 10 selected irregularly shaped pores..................................................................................................................................................................... 57

Table 3

 3 Frobenius distance in % between numerical compliance contribution tensors for irregular and ellipsoidal pores compared to analytical solution based on Eshelby's theory; Conductivity problem: Frobenius distance in % between numerical resistivity contribution tensors for both irregular and ellipsoidal pores compared to analytical solution based on implicit function of the ellipsoid. .............. 143Ce travail de thèse se place dans le cadre de l'étude de l'influence de la microstructure des roches poreuses hétérogènes sur le comportement à l'échelle macroscopique. En général, les comportements court et long terme des roches sont modélisés soit par des approches macromécaniques (i.e., phénoménologiques) soit par des approches micromécaniques. Les approches micromécaniques ont un avantage particulier et pratique puisqu'un nombre relativement faible d'hypothèses est nécessaire. En effet, dans ce cas la complexité du comportement macroscopique simulé ne dérive pas de la complexité du formalisme introduit, mais plutôt de quelques hypothèses 'objectif de la thèse est la caractérisation expérimentale micro/macromécanique multiéchelle des roches poreuses hétérogènes oolithiques. Deux types de roches sont étudiés : le calcaire de Lavoux et le minerai de fer. Ces roches sont composées par un assemblage de grains poreux plus ou moins sphériques (oolithes), de pores et de cristaux inter-granulaires. On s'intéresse à la caractérisation de la géométrie des constituants de la microstructure de ces roches à l'aide des techniques d'imagerie tridimensionnelles comme la microscopie électronique à balayage (MEB) et la nanotomographie 3D aux rayons X. Le MEB permet d'obtenir plusieurs informations caractéristiques du matériau comme la topographie, la morphologie et la composition chimique de

	La microstructure de deux types de roches oolitiques hétérogènes, calcaire de Lavoux
	(Figure 0.1 ; 0.3) et minerai de fer (Figure 0.2 ; 0.4), est étudiée par le biais de la microscopie
	électronique à balayage (MEB) et la nano-tomographie 3D aux rayons X.

. 1: Volumetric and surface mesh information of Cube-IrregularPore and Cube-Ellipsoid models. .................................................................................................................................................................... 69 Table 3. 2: The six independent loadings considered to obtain the 21 independent components of compliance contribution tensor. .................................................................................................................. 71 Table 3. 3: Frobenius distance in % between analytical and numerical compliance contribution tensors for irregular and ellipsoidal pores compared to the analytical solution based on Eshelby's theory. .......... 78 Table 3. 4: Volumetric and surface mesh information of Cube-IrregularPore and Cube-Ellipsoid models using Mesh A. ............................................................................................................................................. 80 Table 3. 5: Frobenius distance in % between numerical compliance contribution tensors for irregular pore and ellipsoidal approximations using PCA + Same volume and AFM methods. ............................... 88 Table 4. 1: Loadings and boundary conditions considered to obtain the components of resistivity contribution tensor .................................................................................................................................... 107 Table 4. 2: Frobenius distance in % between analytical and numerical resistivity contribution tensors .................................................................................................................................................................. 109 Table 4. 3: Frobenius distance in % between numerical contribution tensors for irregular and ellipsoidal pores compared to the analytical solution. ................................................................................................ 110 Table A. 1: Properties of 10 selected pores from the iron ore porous network. ................................... 133 XVI Table D. 1: Elastic problem: XVII Résumé de la thèse 0.1 Contexte général et problématique : sur la microstructure et de considérations statistiques sur les éléments constitutifs (grains, cristaux, pores, fissures). Le comportement macroscopique thermomécanique des matériaux est affecté par la microstructure. La relation micro-macroscopique est décrite par le biais des modèles micromécaniques développés dans le cadre de la théorie de l'homogénéisation des milieux hétérogènes (méthode d'homogénéisation de Maxwell par exemple). Les méthodes d'homogénéisation (changement d'échelle) sont utilisées pour définir des propriétés effectives (élastiques, conductivité thermique), à l'échelle "macroscopique" où le matériau hétérogène peut être considéré comme un matériau homogène équivalent. Ces méthodes utilisent les informations microstructurales comme les propriétés physiques des constituants, l'effet de forme, la taille et la distribution pour calculer ces propriétés effectives. Les applications envisagées dans ce travail concernent principalement les matériaux hétérogènes de type roches poreuses. On s'intéressera en particulier aux roches poreuses oolithiques comme le calcaire de Lavoux et le minerai de fer largement étudiées au laboratoire GeoRessources dans le cadre de l'étude de faisabilité du stockage géologique du dioxyde de carbone CO2 (Sterpenich et al. 2009, Grgic 2011) et de l'étude des effondrements miniers de Lorraine (Grgic et al. 2013, Dagallier et al. 2002, Maitte et al. 2015), et donc une description détaillée des propriétés minéralogiques, pétro-physiques et hydromécaniques de ces roches existe et elle peut être utilisée comme référence dans ce travail. XVIII 0.2 Objectifs de la thèse Lla surface alors que la nanotomographie permet la caractérisation du réseau poreux (taille, répartition spatiale et fraction volumique), des oolithes et des cristaux inter-oolithiques. Les modèles micro-macroscopiques utilisés sont enrichis par les données microstructurales obtenues à partir des observations microscopiques, d'où l'originalité du travail effectué pendant cette thèse qui est due à la combinaison unique d'observations microstructurales (MEB, Nano-tomographie 3D aux rayons X), d'expériences micro-mécaniques (nanoindentation) et de modélisations micro/macroscopiques. En plus, ce travail est innovant car il s'adresse aux matériaux naturels (roches) avec des microstructures réalistes (hétérogènes) contenant des inhomogénéités (pores, inclusions) de différentes formes, orientations et propriétés.

La plupart des méthodes d'homogénéisation des milieux aléatoires se basent sur les approches d'Eshelby et donc sur la solution fondamentale de l'inhomogénéité ellipsoïdale isolée en milieu infini obtenue par

Eshelby (1957Eshelby ( , 1961))

. Un modèle simplifié dans le cadre de la méthode d'homogénéisation de Maxwell est considéré où les oolithes poreuses sont approximées par des sphères, et les pores de formes irrégulières sont approximés par des ellipsoïdes. Cette approximation est réalisée grâce à la méthode d'analyse en composantes principales (ACP) qui fournit les propriétés géométriques telles que la longueur des demi-axes et l'orientation des ellipsoïdes résultants. Cette approximation est ensuite vérifiée par le biais de la reformulation du modèle de Maxwell en fonction des tenseurs de contribution de souplesse et de résistivité thermique

(Sevostianov and Giraud, 2013)

. XIX 03. Observations microscopiques et propriétés micromécaniques du calcaire de Lavoux et du minerai de fer

Table 1

 1

. 1: Technical specifications of nanoidenter tester

Table 1

 1

	. 2: Indentation series and parameters.
	Set	I	II	III
	Maximal load (mN)	75	40	125
	Loading rate (mN/min)	100	100	150
	Unloading rate (mN/min)	100	100	150

Table 1

 1 It is used to calculate the various properties of the voxels in the selected volume or ROI. A fast calculation of the porosity can be made using this tool based on black gray values of the volume.I.3.4 Microstructural components of the Lavoux limestone and the iron oreIn order to distinguish different components of the Lavoux limestone, 3D description of the material is needed thanks to 2D and 3D X-Ray images obtained from tomography. This process applied first on the Lavoux limestone is done using the pixel gray value histogram where each gray value corresponds to a specific value of density. Therefore, the first step is to perform a segmentation of different components after the removing of the noise background. The total volume obtained was 600 mm 3 . Figure1.15 shows the histogram evaluating pixel gray values.

	. 3: Technical specifications of the nano-tomograph.
	Max Voltage/power	Geom. Magnification (3D)	Min voxel size	Rotation	Max sample diameter
	180 Kv / 20 W	1.5 × -300 ×	Down to 0.3 microns	0° -360° × n	< 1 mm to 240 mm

Figure 1. 12: Region of interest selected for the volumetric reconstruction of X-Ray images.

Table 1. 4: Tomographic scan settings.

Timing rate (ms per image) Voltage (Kv) Current (mA) The CT-Scan delivers the internal structure of the scanned sample in a series of 2D slices encoded as 16 bits per pixel data. Tomographic images are represented by pixel gray value histogram where each gray value represents the attenuation of the X-ray beam from the scanner.

To analyze the tomographic images, we used VGStudio MAX 2.2, an image processing software developed by Volume

Graphics (www.volumegraphics.com)

. This software package is employed for visualizing and analyzing voxel data. It is used in a wide range of applications such as industrial CT, medical research, life sciences, and many more. In order to distinguish the different

  constructing a function from the sample called indicator function χ in order to extract the most appropriate iso-surface. The most relevant work describing Poisson surface reconstruction is presented in

	cov( , ) cov( , ) cov( , ) x x x y x z 
	cov( , ) cov( , ) cov( , ) cov( , ) cov( , ) cov( , ) y x y y  y x z x z y z z   
	i 	i 	1, 2,3

and many others. In our work, we used "Poisson's algorithm", one of the most famous methods in 3D surface reconstruction based on implicit function method. Poisson reconstruction method consists of

Table 2 .

 2 1: Properties of 10 selected pores from the Lavoux limestone porous network.

	Pore	Volume (mm 3 )	Diameter (mm)	Voxels	Surface (mm 2 )
	1	0.0010789	0.29	8631	0.112
	2	0.0023486	0.3	18789	0.188
	3	0.0034797	0.41	27838	0.246
	4	0.0002641	0.06	2113	0.037
	5	0.00138	0.35	10770	0.143
	6	0.0001166	0.18	933	0.033
	7	0.0006767	0.19	5414	0.072
	8	0.0005477	0.22	4382	0.065

Table 2 .

 2 2: Surface reconstruction and ellipsoidal approximation for 10 selected irregularly shaped pores.

	9933 0.0878 0.0742  
	0.0623 0.1313 0.09893 0.0966 0.9874 0.1249     

Figure 2. 22: Example of ellipsoidal approximation for a selected 3D irregular pore.

Table 3 .

 3 2: The six independent loadings considered to obtain the 21 independent components of compliance contribution tensor.

	uniaxial compression
	Loadcase	direction	south	north	east	west	top bottom	xy	mid-plan	xz	mid-plan	yz	mid-plan
	1	x	𝜎 11	𝜎 11						dz=0		dy=0		dx=0
	2	y			𝜎 22	𝜎 22				dz=0		dy=0		dx=0
	3	z					𝜎 33	𝜎 33		dz=0		dy=0		dx=0
	shear deformation
	Loadcase	plane	displacement	south	north		east		West		top		bottom

  Analytical compliance contribution tensor based on first Eshelby problemWe follow the paper ofKachanov et al. (

	.3754  0.3774 0.4341 0 0 0 
	0.3774 1.9807  0.6171 0 0 0  
	0.4341 0.6171 3.0447 0 0 0  0 0 0 3.1473 0 0   	(3.24)
	0  0 0 0 2.7669 0 0 0 0 0 0 2.082   
	III.2.3.2

Table 3 .

 3 3: Frobenius distance in % between analytical and numerical compliance contribution tensors for irregular and ellipsoidal pores compared to the analytical solution based on Eshelby's theory.

	Pore	Ellipsoidal approximation	Volume (mm 3 )	Ellipsoid (implicit function)	Ellipsoid (numerical)	Pore (numerical)
			0.0004952	2.6	2.8	3.9
			0.0023486	2.92	2.93	2.96
			0.0017	2.85	2.97	3.71
			0.0005874	1.974	2.21	4.52
			0.0006579	0.18	0.61	2.04
			0.0002641	1.21	2.35	1.12
			0.0006767	0.98	2.17	3.12

III.3 Effect of geometrical parameters on compliance contribution tensors III.3.1 Effect of mesh size

Table 3 .

 3 4: Volumetric and surface mesh information of Cube-IrregularPore and Cube-Ellipsoid models using Mesh A.

	Mesh	Nodes	Surface elements	Volumetric elements
	A
	Cube + IrregularPore	59069	7266	41404
	Irregular pore	2102	1050	-
	Cube-Ellipsoid	92104	11052	64313
	Ellipsoidal pore	7046	3521	-
	B
	Cube + IrregularPore	191662	13406	137067
	Irregular pore	7030	3514	-
	Cube + Ellipsoid	212110	25453	186657
	Ellipsoidal pore	16970	7635	-
	C
	Cube + IrregularPore	932482	53023	676611
	Irregular pore	23638	11818	-
	Cube + Ellipsoid	1252300	150276	1102024
	Ellipsoidal pore	100184	48089	-

  Ellipsoidal approximation of pores using PCA with conservative volume PCA method, used before in mechanical problem, is one of the best methods to approximate irregular shapes. However, this method based on inertia moments does not conserve the volume

	S	ellip		4 	1/   3 p p p p p p a b a c b c   	p	with	p 	1.6075	.	(3.37)
	III.3.2.2

Table 3 .

 3 5: Frobenius distance in % between numerical compliance contribution tensors for irregular pore and ellipsoidal approximations using PCA + Same volume and AFM methods.

	Pore	Hsame_V	HAFM
	1	3.5	3.9
	2	4.5	4.7
	3	2.9	3.3
	4	3.7	4.2
	5	3.8	4
	6	4.4	5.1
	7	4.1	4.8
	8	4.5	4.9
	9	2.6	3
	10	3.6	4

  The dependencies of the bulk modulus eff k and the shear coefficient eff  on the volume fraction  for the material containing (randomly oriented) irregularly shaped pores are presented in Figures 4.6 and 4.7. In addition, we estimated effective properties to the same material containing the 10 approximated ellipsoids, also randomly oriented (Figure 2.15). Each line corresponds to a model containing one of the 10 considered pores randomly oriented. For both cases (irregular and ellipsoidal pores), the volume fraction  of pores varies from 0 to 0.3, so that the effective bulk modulus eff k is equal to the bulk modulus of the solid matrix m k in the absence of pores and it decreases with the increase in porosity. Shear coefficient also decreases with the increase in porosity. Figures 4.6 to 4.9 show then that the influence of the volume fraction on elastic stiffness coefficients is significant. Furthermore, for every imposed volume fraction, most values of bulk and shear moduli are very close, so the difference (less than 4% from the mean values) in effective properties between the 10 selected pores is not very significant. This is due to similarities in geometry shape parameters of the pores. In Figures 4.10 and 4.11, the values of bulk

	12𝑘 𝑐 𝜇 𝑐 3𝑘 𝑐 +4𝜇 𝑐	𝕁 +	2 5	𝜇 𝑐 (9𝑘 𝑐 +8𝜇 𝑐 ) 3𝑘 𝑐 +4𝜇 𝑐	𝕂
	IV.1.2.3 Results and discussion:	(4.21)

Table 4

 4

	. 2: Frobenius distance in % between analytical and numerical resistivity contribution tensors
	of pore shown in Figure 3.3.
	H R tensor	Ellipsoid (implicit function)	Ellipsoid (vector product)	Irregular pore
	Ellipsoid (implicit function)	0	0.52 %	8.18
	Ellipsoid (vector product)		0	8.12

Table 4 .

 4 3: Frobenius distance in % between numerical contribution tensors for irregular and ellipsoidal pores compared to the analytical solution.

	Pore	Ellipsoidal approximation	Volume (mm 3 )	Ellipsoid (numerical)	Pore (numerical)
			0.0004952	0.78	4.4
			0.0023486	0.5	6.27
			0.0017	0.55	8.5
			0.0005874	0.4	9.7
			0.00138	0.52	8.18
			0.0002641	0.52	9.6
			0.0006767	0.33	9.27

  denotes the number of Gauss integration points in a T6 finite element set to 4 (TableB.1).

	op h	( , )  	11 22 g g  ( , ) op f  	2 12 g	(B.10)
	The total surface integration is then calculated thanks to standard Gauss Legendre quadrature rules
	1 00 d   1      with el op I G n el i op G op i I w h  ( , ) op h d    i i G G    that writes for one element 1 ( , ) where G i i i G G G w  n ,, represent respectively the local coordinates of Gauss integration points and (B.9) 137 (B.11) corresponding weight.

  Table B.1: Coordinates (ξG; ηG) and weights ωG of Gaussian integration points for the 6-nodes isoparametric Lagrange finite element

Acknowledgements

Firstly, I would like to express my sincere gratitude to my advisor Dr. Dragan Grgic, my co-advisor Dr. Christophe Auvray and Prof. Albert Giraud for the continuous support of my PhD study and related research, for their patience, motivation, and immense knowledge. Their guidance helped me in all the time of research and writing of this thesis. I could not have imagined having a better advisors and mentors for my PhD study.

I would like to thank the rest of my thesis committee: Prof. Dashnor Hoxha, Prof. Jian Fu Shao,

The parameters of the pores obtained are presented in Figure 2.16 and Figure 2.17. Particular attention is given to the radius and the volume of pores to understand their size as well as the angle that they form with the vertical to analyze their distribution inside the REV.

The distribution of porosity for the Lavoux limestone shows the presence of interoolitic pores while the intraoolitic pores are not represented because of the limited resolution of the tomographic scan (5 µm). Interoolitic pores radius varies between 6 and 130 µm with a concentration between 6 and 25 µm with more than 70% of the radii and more than 80% of the pores have a volume between 1000 and 11000 µm 3 . For the iron ore, the pores are clearly smaller than those of the Lavoux limestone and they vary between 1 and 166 with a strong presence of small pores until 10 µm. The tomography scan detected smaller pores than the Lavoux limestone since the resolution was 2.4 µm. However, intraoolitic couldn't been detected and therefore they are not represented in the distribution graph. More than 90% of iron ore pores have a volume between 10 and 150 µm 3 .

To study the distribution of pores inside the REV, we calculated the angle with the vertical for 30 pores randomly selected from binary data of the Lavoux limestone. These angles are calculated thanks to the formulas of Euler angles which describe the orientation of the pore with respect to

Chapter III

Compliance contribution tensor of 3D irregularly shaped pores Cube-IrregularPore and cube-Ellipsoid models were discretized using quadratic isoparametric Lagrange finite elements (Figure 3.5):

• 10-nodes tetrahedrons have been used for volume discretization of the matrix.

• 6-nodes triangles (T6) have been used to discretize the surface of the pore and the ellipsoid.

These 6 nodes triangles represent the borders of volumetric finite elements composing the pore boundary.

Figure 3. 5: Irregularly shaped pore and approximated ellipsoid volume and surface mesh inside a reference cube using quadratic elements.

Special care was taken during meshing to avoid distorted elements or orphan nodes (nodes do not belong to any mesh) because otherwise, simulations cannot be performed due to mesh errors.

As a result, one refers to Table 3.1 for details about the number of nodes and elements for both meshes (irregular pore and ellipsoid):

Chapter IV

Determination of effective properties with Maxwell homogenization scheme

• Loadcase 2 (uniform thermal flux in y direction).

• Loadcase 3 (uniform thermal flux in z direction). Txx -Txx q = 0 q = 0 q = 0 q = 0 2 q = 0 q = 0 Tyy -Tyy q = 0 q = 0 3 q = 0 q = 0 q = 0 q = 0 Tzz -Tzz

where Tii is the imposed temperature on cube faces to generate a gradient and q is the heat flux. We used the same mesh generated for the mechanical problem (Figure 3.5) and simulations were performed using Code_Aster. Thus, a new command file has to be written so that the phenomenon was changed to thermal one instead of mechanical. In the result file, the temperature on each node

Conclusion and Perspectives

The presented work is focused on the influence of the microstructure of heterogeneous porous rocks on the macroscopic behavior.

Oolitic porous rock such as Lavoux limestone and iron ore were chosen as a reference due to their simple composition.

Microscopic observations by SEM and computed nantomography showed that these rocks are composed of three main components: oolites, interoolitic pores and interoolitic cement.

Statistical information on the components of the rocks are obtained by using a segmentation algorithm based on gray scale values to distinguish different components of the material. The distribution of the volume and the radius showed that the Iron ore material is more porous than the Lavoux limestone. However, the size of pore in the Lavoux limestone is clearly larger which means that iron ore material contains more intra-oolitic pores. These intra-oolitic pores are not represented in the distribution due to the limitation of the tomography scan.

Micromechanical properties such as hardness and elastic modulus were obtained using nanoindentation tests. These tests showed that the Lavoux limestone is composed of three families of calcite with different sizes: The inner layers of oolites (microcalcite), the outer layers of oolites (microcalcite) and the inter-oolitic sparry calcite (macro-calcite). The values of elastic modulus of each family showed that the macro calcite is the hardest one while the inner layers od oolites are the weakest. This conclusion is due to the effect of porosity where the inner layers of oolites have higher porosity than the other two families.

Tomographic images showed that oolites have quasi-spherical shapes while the pores have irregular shapes. Therefore, analytical ellipsoidal homogenization method based on 1 st Eshelby solution cannot be used. Hence, a simplified model within the framework of Maxwell homogenization scheme is applied where oolites are approximated by spheres and irregularly shaped pores are approximated by tri-axial ellipsoids. PCA method was used to perform this approximation.

The approximation of oolites by spheres is verified by calculating their sphericity ratio. Calculated sphericity of oolites was close to 1 for Lavoux limestone and close to 0.6 for iron ore and therefore the approximation of oolites by spheres is reasonable.

• Amenta, N., Bern, M., Kamvysselis, M. (1998, July). A new Voronoi-based surface reconstruction algorithm. In Proceedings of the 25th annual conference on Computer graphics and interactive techniques (pp. 415-421). ACM. • Amenta, N., Choi, S., Dey, T. K., Leekha, N. (2000, May). A simple algorithm for homeomorphic surface reconstruction. In Proceedings of the sixteenth annual symposium on Computational geometry (pp. 213-222). ACM. 1829). Sur l'équation à l'aide de laquelle on détermine les inégalités séculaires des mouvements des planètes, vol. 9. Oeuvres Complètes (IIème Série); Paris: Blanchard. • Chen, F., Giraud, A., Grgic, D., Kalo, K. (2017). A composite sphere assemblage model for porous oolitic rocks: Application to thermal conductivity. Journal of Rock Mechanics and Geotechnical Engineering, 9(1), 54-61. • Cloetens, P., Ludwig, W., Boller, E., Helfen, L., Salvo, L., Mache, R., Schlenker, M. (2002, January). Quantitative phase contrast tomography using coherent synchrotron radiation.

In International Symposium on Optical Science and Technology (pp. 82-91). International Society for Optics and Photonics. • Cloetens, P., Pateyron-Salomé, M., Buffiere, J. Y., Peix, G., Baruchel, J., Peyrin, F., & Schlenker, M. (1997). Observation of microstructure and damage in materials by phase sensitive radiography and tomography. Journal of Applied Physics, 81(9), 5878-5886. • Constantinides, G., Chandran, K. R., Ulm, F. J., Van Vliet, K. J. (2006). Grid indentation analysis of composite microstructure and mechanics: Principles and validation. Materials Science and Engineering: A, 430(1), 189-202. • Constantinides, G., Ulm, F. J. ( 2007). The nanogranular nature of C-S-H. Journal of the Mechanics and Physics of Solids, 55(1), 64-90. • Curless, B., Levoy, M. (1996, August). A volumetric method for building complex models from range images. In Proceedings of the 23rd annual conference on Computer graphics and interactive techniques (pp. 303-312). ACM. • Dey, T. K., Goswami, S. (2003, June). Tight cocone: a water-tight surface reconstructor.

In Proceedings of the eighth ACM symposium on Solid modeling and applications (pp. 127-134). ACM.

• Diepenbroek, M., Bartholomä, A., Ibbeken, H. (1992). How round is round? A new approach to the topic 'roundness' by Fourier grain shape analysis. Sedimentology, 39(3), 411-422. • Dini, U. (1878). Fondamenti per la teorica delle funzioni di variabili reali. Tip. T. Nistri. (1961). Calculation of the effective isotropic moduli of disperse solid systems.

Poroshkovaja Metallurgija (Powder Metallurgy), #1, 50-55 (in Russian). • Sneddon, I. N. (1965). The relation between load and penetration in the axisymmetric Boussinesq problem for a punch of arbitrary profile. International journal of engineering science, 3(1), 47-57. • Sneed, E. D., Folk, R. L. (1958). Pebbles in the lower Colorado River, Texas a study in particle morphogenesis. The Journal of Geology, 66(2), 114-150. • Sterpenich, J., Sausse, J., Pironon, J., Gehin, A., Hubert, G., Perfetti, E., Grgic, D. (2009).

Experimental ageing of oolitic limestones under CO2 storage conditions. Petrographical and chemical evidence, Chem. Geol.,265 (1-2):99-112. • Stewart, G. W. (1993). On the early history of the singular value decomposition. SIAM review, 35(4), 551-566.

• Suzuki, T., Ohmura, T. (1996). Ultra-microindentation of silicon at elevated temperatures. Philosophical Magazine A, 74( 5 • Thompson, A. H., Katz, A. J., Raschke, R. A. (1987) • Zhao, H. K., Osher, S., Fedkiw, R. (2001). Fast surface reconstruction using the level set method.

In Variational and Level Set Methods in Computer Vision, 2001. Proceedings. IEEE Workshop on (pp. 194-201). IEEE. • Zhu, Y. T., Beyerlein, I. J. ( 2002). Bone-shaped short fiber composites-an overview. Materials Science and Engineering: A, 326(2), 208-227. • Zienkiewicz, O. C., Taylor, R. L. (2000). The finite element method: solid mechanics (Vol. 2).

Butterworth-heinemann.

• Zimmerman, R. W. (1989). Thermal conductivity of fluid-saturated rocks. Journal of Petroleum Science and Engineering, 3(3), 219-227.

Appendix A: Geometrical properties of 10 selected pores for iron ore material.

Appendix B: numerical integration on the surface of irregularly shaped pore

Irregular pore is discretized with quadratic isoparametric 6-node Lagrange triangular finite elements (T6). The posed problem consists therefore of performing numerical integration on the curved surface of these finite elements. We define a reference plane on element surface given by 2 curvilinear coordinates ( , )

 .

For the 6 -nodes triangular element, i N denotes the global numbering of the node in the 3D finite element mesh,

correspond to vertices of triangular elements, and nodes

))

correspond to midedge nodes.

Thus, shape functions of the quadratic isoparametric 6-node Lagrange triangular finite element allow expressing global coordinates 1 2 3 ( , , ) z z z as functions of the two local (curvilinear) coordinates ( , )  of the reference plane. Indeed, if z is the position vector on a T6 surface element, by using Gauss Legendre quadrature rules, position vector coordinates i z may be expressed as function of local curvilinear coordinates ( , )

where () i j z denotes the cartesian coordinate j z of node () i .

First partial derivative of spatial coordinates writes:

;

Appendix C: Analytical solution for compliance contribution tensors of ellipsoidal pores.

Compliance contribution tensor depends on the parameters of the ellipsoid but it is independent of the initial stress-strain state on the matrix. The relations between Hill's tensor and its counterpart for the compliance problem ℚ are given as follows:

Thus, the compliance contribution tensor is given as follows:

where N is a constant shape dependent fourth-order tensor given by:

Strain concentration tensor 𝔻 introduced by Wu (1966), is defined as the dilute concentration tensor of a homogeneous ellipsoidal inhomogeneity of aspect ratio 1  embedded in an infinite isotropic matrix:

Relations between 𝔻 and ℕ tensors write:

By inverting the two members of the equality, one can verify:

However, for the case of an ellipsoidal pore embedded in a matrix, the stiffness tensor for the pore is equal to zero (ℂ 1 = 0):

And the compliance contribution tensor is then given by:

For an isotropic matrix: Assuming that 𝜌 1 ≥ 𝜌 2 ≥ 𝜌 3 and renaming 𝜌 1 = 𝑎, 𝜌 2 = 𝑏, and 𝜌 3 = 𝑐 because these notations will be used in the rest of the work for the semi axes of the tri-axial ellipsoid, one can write (Kachanov et al., 2003): where ℱ(𝜃, 𝑘) is the elliptic integral of the first kind of amplitude and parameter: Keywords: Property contribution tensor, Microstructural analysis, Finite element method, 3D irregular shape inhomogeneities, Effective properties, X-ray 3D Nanotomography, Shape analysis, Homogenization.

Abstract:

The aim of this work is to study the influence of the microstructure of heterogeneous porous rocks on the behavior at the macroscopic scale. Thus, we characterized the microstructure and micromechanical properties (thanks to nano-indentation tests) of two porous oolitic rocks (Lavoux limestone and iron ore) to calculate their effective mechanical and thermal properties.

Oolitic rocks are constituted by an assemblage of porous grains (oolites), pores and inter-granular crystals. Scanning electron microscopy and X-ray 3D Computed Tomography were used to identify the different components of these rocks. Particular attention was given to X-Ray computed tomography since this analytical method allows the characterization of the porous network (size, spatial distribution, and volume fraction), and the shapes of oolites and inter-oolitic crystals. The novelty of this work lies in taking into account the 3D real shape of pores. Hence, we approximated porous oolites by spheres and irregularly shaped pores by ellipsoids. This approximation was performed thanks to the principal component analysis (PCA), which provides the geometrical properties such as length of semi-axes and orientation of resulting ellipsoids. The sphericity of the approximated oolites was calculated and the values close to 1 allowed us to consider oolites as spheres.

To verify the approximation in the case of pores, we evaluated the contribution of these irregularly shaped three-dimensional pores to the overall elastic properties. Thus, compliance contribution tensors for 3D irregular pores and their ellipsoidal approximations were calculated using the finite element method (FEM). These tensors were compared and a relative error was estimated to evaluate the accuracy of the approximation. This error produces a maximum discrepancy of 4.5% between the two solutions for pores and ellipsoids which verifies the proposed approximation 146 procedure based on PCA. The FEM numerical method was verified by comparing the numerical solution for compliance contribution tensors of ellipsoids to the analytical solution based on

Eshelby's theory. The difference between these two solutions does not exceed 3%. The same numerical method was used to calculate thermal resistivity contribution tensors.

Calculated compliance and resistivity contribution tensors were used to evaluate effective elastic properties (bulk modulus and shear coefficient) and effective thermal conductivity by considering the two-step Maxwell homogenization scheme. The results showed an important influence of the porosity on effective properties. Finally, the results obtained for irregular pores were compared to those for ellipsoidal ones and they showed a good agreement with a maximum deviation of 4% which verifies once again the approximation of irregularly shaped pores by tri-axial ellipsoids.

Titre : Caractérisation microstructurale et modélisation numérique de roches poreuses oolithiques.

Mots clés : Tenseur de contribution, Analyse microstructurale, Méthode des éléments finis, Hétérégonénéités 3D de forme irrégulière, Propriétés effectives, Nanotomographie 3D aux rayons X, Analyse de forme, Homogénéisation.

Résumé :

Le but de ce travail est d'étudier l'influence de la microstructure de roches poreuses hétérogènes sur le comportement à l'échelle macroscopique. Ainsi, nous avons caractérisé la microstructure et les propriétés micromécaniques (grâce à des tests de nano-indentation) de deux roches oolithiques poreuses (calcaire de Lavoux et minerai de fer) pour calculer leurs propriétés mécaniques et thermiques effectives.

Les roches oolithiques sont constituées d'un assemblage de grains poreux (oolithes), de pores et de cristaux intergranulaires. La microscopie électronique à balayage et la tomographie 3D aux rayons X ont été utilisées pour identifier les différents composants de ces roches. Une attention particulière a été accordée à la tomographie aux rayons X car cette méthode analytique permet de caractériser le réseau poreux (taille, distribution spatiale et fraction volumique), ainsi que la forme des oolithes et des cristaux inter-oolithiques. La nouveauté de ce travail réside dans la prise en compte de la forme 3D réelle des pores. Par conséquent, nous avons approximé les oolites poreuses par des sphères et les pores de forme irrégulière par des ellipsoïdes. Cette approximation a été réalisée grâce à l'analyse en composantes principales (ACP), qui fournit les propriétés géométriques telles que la longueur des demi-axes et l'orientation des ellipsoïdes résultants. La sphéricité des oolites approximées a été calculée et les valeurs proches de 1 nous ont permis de considérer les oolithes comme des sphères.

Pour vérifier l'approximation dans le cas des pores, nous avons évalué la contribution de ces pores tridimensionnels de forme irrégulière aux propriétés élastiques et thermiques effectives.

Ainsi, les tenseurs de contribution de souplesse pour les pores irréguliers 3D et leurs approximations ellipsoïdales ont été calculés en utilisant la méthode des éléments finis (FEM). Ces tenseurs ont été comparés et une erreur relative a été estimée pour évaluer la précision de l'approximation. Cette erreur produit une distance maximale de 4,5% entre les deux solutions pour les pores et les ellipsoïdes, ce qui vérifie la procédure d'approximation proposée basée sur ACP.

La méthode numérique FEM a été vérifiée en comparant la solution numérique des tenseurs de contribution des ellipsoïdes à la solution analytique basée sur la théorie d'Eshelby. La différence entre ces deux solutions ne dépasse pas 3%. La même méthode numérique a été utilisée pour calculer les tenseurs de contribution de résistivité thermique.

Les tenseurs de souplesse et de résistivité calculés ont été utilisés pour évaluer les propriétés élastiques effectives (module élastique et coefficient de cisaillement) et la conductivité thermique effective en considérant le schéma d'homogénéisation de Maxwell en deux étapes. Les résultats ont montré une influence importante de la porosité sur les propriétés effectives. Enfin, les résultats obtenus pour les pores irréguliers ont été comparés à ceux des ellipsoïdes et ils ont montré un bon accord avec un écart maximal de 4% ce qui vérifie l'approximation des pores de forme irrégulière par des ellipsoïdes triaxiaux.