Je Remercie Tous Particulièrement Monsieur

Djimédo Kondo

Minh-Tan Nguyen

Quy-Dong To

Vincent Monchiet

Guy Bonnet

Computation of the Brinkmann coefficients with FFT methods

jury d'avoir accepté de juger ce travail.

Conclusion générale et perspectives

Bibliographie

Table des matières

Résumé 117

Abstract 119

Table des figures of α and β. The results are obtained for two materials with the following parameters ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 3. The optimal values by the two methods are respectively (α, β) = (0.355, 0.298) and (α, β) = (0.309, 0.349) 2.3 Different arrangements of non-overalapping spheres in cubic unit cell. From left to right : simple cubic (SC), body centered cubic (BCC), face centered cubic (FCC) and random distribution (RD). The random distribution is generated by Event Driven Molecular Dynamics method. . . . . . . . . . . 2.4 Normalized effective bulk modulus κ e /κ 0 vs inclusion volume fraction f of BCC array. Elastic properties of the constituents are ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 10. The results are computed by first order estimates which all coincide with HS estimates, second order estimates of the three methods (OR,ON and OD) and the numerical method FFT at convergence. . . . . . 2.5 Normalized effective shear moduli µ e /µ 0 vs inclusion volume fraction f of BCC array. Elastic properties of the constituents are ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 10. The results are computed by first order estimates, second order estimates of the three methods (OR,ON and OD) and the numerical method FFT at convergence. . . . . . . . . . . . . . . . . . . . . . . . . . .

Table des figures

2.6 Normalized effective bulk modulus κ e /κ 0 vs inclusion volume fraction f of FCC array. Elastic properties of the constituents are ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 10. The results are computed by first order estimates, second order estimates of the three methods (OR,ON and OD) and the numerical method FFT at convergence. . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Normalized effective shear moduli µ e /µ 0 vs inclusion volume fraction f of FCC array. Elastic properties of the constituents are ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Normalized effective bulk modulus κ e /κ 0 vs inclusion volume fraction f for random distribution of rigid spheres (µ 1 /µ 0 = κ 1 /κ 0 = ∞). The solutions of the present work (OD-1, OD-2) are compared with the results of [START_REF] Torquato | Effective stiffness tensor of composite media. I. Exact series expansions[END_REF][START_REF] Torquato | Effective stiffness tensor of composite media : II. Applications to isotropic dispersions[END_REF].

The first order estimates coincide with HS estimates. . . . . . . . . . . . . 2.9 Normalized effective bulk modulus κ e /κ 0 vs inclusion volume fraction f with random distributions of spherical voids (µ 1 /µ 0 = κ 1 /κ 0 = 0). . . . . . parison of the second order OS-2 estimations (with and without S (3) approximation) and Torquato's expression [START_REF] Torquato | Effective electrical conductivity of two ?phase disordered composite media[END_REF] (TO), Kim and Torquato [START_REF] Ic | Effective conductivity of suspensions of hard spheres by brownian motion simulation[END_REF] (KT), the simulation data of Bonnecaze and Brady [START_REF] Rt Bonnecaze | The effective conductivity of random suspensions of spherical particles[END_REF] (BB). The Clausius-Mossoti (CM) approximation coincides with first order OS approximation OS-1 (n = 1). .

Notations

• Tensory notation a scalar, a vector, A second order tensor, A fourth order tensor, I second order identity tensor, I fourth order identity tensor, .

contracted product of order one, : contracted product of order two, ⊗ tensor product, δ ij Kronecker symbol,

(A⊗B) ijkl = 1 2 (A ik B jl + A il B jk )
• Notations common to all chapters ξ wave vector, Γ 0 (ξ) Green operator for the strain, ∆ 0 (ξ) Green operator for the stress, σ tensor of the microscopic stresses,

• Notations specific to Chapter 1

T temperature field, e temperature gradient, E macroscopic temperature gradient, q heat flux, Q macroscopic heat flux, K conductivity tensor, R resistivity tensor, χ the characteristic function,

• Notations specific to Chapter 2 C stiffness tensor, C e effective stiffness tensor, C 0 stiffness tensor of reference materials, E 1 , .., E 6 tensors of Walpole base, S flexibility tensor, S 0 flexibility tensor of reference materials,

• Notations specific to Chapter 3 v field of microscopic speeds (velocities), V field of macroscopic speeds, p champ de microscopic pressure, G macroscopic pressure gradient, J macroscopic pressure gradient dimensionless, d rate of deformation, K static permeability tensor, f volume force,

• Notations specific to Chapter 4

C α elastic tensor of the phase α, I α characteristic functions describing volumes V α , S α compliance tensor of the phase α, U the strain elastic energy, W the complementary elastic potential,

Introduction générale

Micromechanics is a discipline that studies the relation between the overall behavior of heterogeneous materials and their constituents. Despite a long history [START_REF] Sd Poisson | Second mémoire sur la théorie du magnétisme[END_REF][START_REF] Maxwell | A treatise on electricity and magnetism[END_REF][START_REF] Sanchez-Palencia | Non-homogeneous media and vibration theory[END_REF][START_REF] Nemat-Nasser | Micromechanics : overall properties of heterogeneous materials[END_REF][START_REF] Torquato | Random heterogeneous materials : microstructure and macroscopic properties[END_REF], the research on this field is still very active with developments from both theoretical and numerical aspects. Theoretical contributions concern the derivation of estimates or bounds from the microstructure informations, for example the properties of the constituents, their spatial distribution, shape and orientation... Many theoretical estimates are based from Eshelby's single inclusion problems [START_REF] Jd Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] and taking account of surrounding interaction following Mori-Tanaka [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] (see for instance the book of or Christensen [START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF], Mura [START_REF] Mura | Micromechanics of defects in solids[END_REF], Nemat-Nasser [START_REF] Nemat-Nasser | Micromechanics : overall properties of heterogeneous materials[END_REF]). Variational principle and energetical theorems are also used to obtain rigorous bounds for the overall properties [START_REF] Hashin | A variational approach to the theory of the effective magnetic permeability of multiphase materials[END_REF][START_REF] Willis | Bounds and self-consistent estimates for the overall properties of anisotropic composites[END_REF][START_REF] Walpole | Elastic behavior of composite materials : theoretical foundations[END_REF]. With the fast development of computers, numerical approaches based on Finite Element Method (FEM) have been widely used to determine the effective properties of heterogeneous materials. Alternative methods based on Fast Fourier Transform (FFT) [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF][START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF][START_REF] Monchiet | A polarization-based fft iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF] are able to deal with complex problems involving realistic microstructures obtained from modern imagery techniques such as microtomography.

The present contribution is devoted to the determination of the effective properties of heterogenous materials applied to different microstructures : composites with reinforcements, polycristals and porous media with different physical problems : Heat transfer, elasticity and mass transfer. Classically, the homogenized properties are determined by solving unit cell problems with the periodic conditions on its boundary. The cornerstone of this work is to formulate the local problem with the Lippmann-Schwinger (LS) equation and to solve it with series expansion methods. Each term of the series can be computed analytically with some hypothesis about the microstructure and simplifications, this is the aim of chapter 1 and 2. Each terms of the series can be also computed numerically by means of the FFT algorithm, this is the subject of chapter 3 and 4.

In chapter 1, we study the estimates of the overall conductivity of two phase material. To this end we derive a class of LS equations in terms of polarization p. The latter is derived by combining linearly two LS equations associated to basic schemes and exploiting the relation between the temperature gradient e, the flux j and the polarization p. For each LS equation in the family, we obtain an estimate of the effective conductivity which contains both a truncated Neumann series and an approximation of its residual. The estimate can be further improved by minimizing the operator norm, i.e its maximal eigenvalue, to ensure the rapid convergence of the series. In the case of spherical suspension, the structure factors at different order have naturally appeared in the estimates. The structure factor is an important statistical information in Fourier space that reveals how the particles are arranged with respect to each other. Numerical examples show that estimates based on the first two order structure factors are close to exact results for a large range of volume fraction and material contrast.

Chapter 2 follows the same strategy as chapter 1 but deals with the elastic behavior of composites. We consider estimates based on the LS equations formulated with the eigenstress and the associated Neumann series. The operator of the latter must be optimized so that the estimates can yield good results. Despite the similarities with the conduction problem, the elasticity one involves fourth order (anisotropic) tensors which raise considerable difficulties. For example, the operator of the Neumann series is the product of characteristic function and a non self-adjoint operator which is an explicit isotropic transverse tensor in Fourier space. In order to bound the spectral radius of the former operator, we propose different optimization techniques : i) minimizing the spectral radius of the non self-adjoint operator (OR), minimizing its norm (ON) and minimizing directly the spectral radius of the former operator (OD). Different optimization method will lead to different Neumann series and estimates of the overall properties. Comparisons with numerical simulations show that the OD estimates always yield the best result. As for the conduction problem, we also show that the structure factors are present in the estimates and the use of the first two order structure factors are satisfactory in many cases.

In chapter 3, we deal with the Stokes flow problem in a porous microstructure. Classically, the effective behavior is depicted with the Darcy equation which involves the determination of a permeability tensor. The associated unit cell problem has been long far established by Auriault [START_REF] Jl Auriault | Study of the macroscopic behavior of a saturated porous elastic medium[END_REF], Sanchez-Palencia [START_REF] Sanchez-Palencia | Non-homogeneous media and vibration theory[END_REF], Levy [START_REF] Levy | Suspension of solid particles in a newtonian fluid[END_REF], etc, using the asymptotic series expansion method. The Darcy equation is retrieved by keeping only the first term of the asymptotic series. Later, Auriault [3] examined the importance of the higher order term of the series following Boutin in the case of elasticity [START_REF] Boutin | Microstructural effects in elastic composites[END_REF]. It found that the higher order terms imply a dependence of the macroscopic filtration law with the double, triple, etc, gradient of pressure. Moreover by making a variable change, Auriault [3] establish the connection with the Brinkman equation [START_REF] Hc Brinkman | A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles[END_REF]. The objective of chapter 3 is to compute numerically with the FFT method the corrective terms of the Darcy equation. The FFT algorithm is inspired by Monchiet et al. [START_REF] Monchiet | A fft-based method to compute the permeability induced by a stokes slip flow through a porous medium[END_REF] that is adapted to compute the solutions of the higher order cell problems. The corrective terms are compared to exact expression for the flow along parallel pores and computed for a regular array of rigid cylinders.

In chapter 4, we derive numerical bounds of the elastic properties of polycristals. The homogenized elastic coefficients are computed from Voronoi-type unit cells. The FFT method combined with the shape functions is used following the method introduced by Bonnet [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF]. The use of the shape function guarantee the bound character of the solutions comparatively to the original method of Moulinec and Suquet, which introduces some uncontrolled approximations. The shape functions for a 2d or 3d polycristals are numerically computed by taking advantage of the recent paper of Wuttke [START_REF] Wuttke | Form factor (Fourier shape transform) of polygon and polyhedron[END_REF]. Some examples are provided for 2d and 3d microstructures of cubic crystals.

Chapitre 1

Conductivity estimates of spherical-particle suspensions based on triplet structure factors

Introduction

Determining effective properties of heterogeneous materials from the constituents' has enjoyed a long history [START_REF] Maxwell | A treatise on electricity and magnetism[END_REF][START_REF] Landau | Electrodynamics of continuous media[END_REF] but is still a very active research area. Due to the variety of the microstructures and the physical phenomena, different approaches have been proposed to evaluate the overall behavior of the homogenized materials. The usual procedure is to solve the localization problem in a RVE (Representative Volume Element) and find the relation between the average quantities such as stress-strain in elasticity or heat fluxtemperature gradient in conduction which is the main focus of the present paper, etc.. The RVE is chosen to represent the material with all heterogeneity details, can be ranging from the unit cell for periodic media [START_REF] Sanchez-Palencia | Non-homogeneous media and vibration theory[END_REF], coated sphere models in micromechanics [START_REF] Gw Milton | The theory of composites[END_REF][START_REF] Nemat-Nasser | Micromechanics : overall properties of heterogeneous materials[END_REF] or a sufficiently large material volume for random media [START_REF] Torquato | Random heterogeneous materials : microstructure and macroscopic properties[END_REF]. For the latter, bounds on the effective conductivity can be established using variational principles [START_REF] Hashin | A variational approach to the theory of the effective magnetic permeability of multiphase materials[END_REF]. Since the localization problem can be written in the form of integral equations [START_REF] Dj Bergman | Dielectric constant of a two-component granular composite : a practical scheme for calculating the pole spectrum[END_REF], it is possible to expand it into Neumann series and truncate to derive the estimations [START_REF] Phan-Thien | New bounds on the effective thermal conductivity of n-phase materials[END_REF]. This type of approximation is termed weak contrast expansion due to the strict convergence condition. Another type of approximation is proposed by Brown [START_REF] Wf Brown | Solid mixture permittivities[END_REF] based on the polarization and cavity field integral equation. This approach have shown to estimate effectively the conductivity for some suspensions of spheres [START_REF] Torquato | Effective electrical conductivity of two ?phase disordered composite media[END_REF][START_REF] Ic | Effective conductivity of suspensions of hard spheres by brownian motion simulation[END_REF][START_REF] Rt Bonnecaze | The effective conductivity of random suspensions of spherical particles[END_REF][START_REF] Rt Bonnecaze | A method for determining the effective conductivity of dispersions of particles[END_REF]. However, this series expansion which is derived from the inversion of a conditionally convergent series, can be subject to some subtle issues such as singularities when the two phases are identical and unclear convergence condition which depends on microstructures, especially near the percolation limit [START_REF] Gw Milton | The theory of composites[END_REF].

In line with these works, the present paper considers two-phase composite materials constituted of identical spherical inclusions embedded in a matrix. The arrangement of the spheres can be either periodic like cubic crystals or randomly distributed like the case of hard sphere suspensions. In the periodic settings, a family of integral equation of polarization p is constructed by linearly combining E based and J based polarization equations [START_REF] Qd To | Closed-form solutions for the effective conductivity of two-phase periodic composites with spherical inclusions[END_REF][START_REF] Qd To | Conductivity of periodic composites made of matrix and polydispersed aggregates[END_REF]. When minimizing the integral operator norm, we obtain an optimal integral equation whose associated Neumann series always converges with fast rate for the whole contrast range. It turns out that the latter is rather similar in form with Eyre and Milton's accelerated scheme [START_REF] Dj Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF]. Consequently, new estimations of the conductivity based on the fast convergent series have been derived in this works. Our strategy, which is closely related to the estimation of the remainder of the series [START_REF] Qd To | A numerical-analytical coupling computational method for homogenization of effective thermal conductivity of periodic composites[END_REF], is to apply n times the recurrence relation to obtain an integral equation containing n leading terms of the Neumann series and make approximation. This provides an expression of the effective conductivity tensor involving correlation functions at different orders. Next, approximations can be made to estimate these correlation functions. Compared to previous works in this field, a new original feature is to use information on correlations in the Fourier domain, more particularly the structure factor, and to combine this information with a series expansion whose convergence is optimized. For n = 2 (two leading terms), numerical results show that such approach works very well for cubic crystal system at a large range of volume fraction and conductivity contrast. For random distribution of hard sphere, it is interesting to note that the derived estimation is explicitly related to structure factors at different orders. Our results compare relatively well with the other existing estimations in the literature.

To develop the above ideas, the content of our work is divided into 4 parts. After the Introduction, we present the theory of integral equation based on the polarization, how to use it to estimate the overall conductivity, the treatment of cubic crystal arrangement or random materials as well as the approximation of the triplet structure factor. After the section 2, section 3 is dedicated to numerical application of the theory and compare with the available theoretical and numerical results from the literature. Finally, the paper ends with the conclusions and remarks.

1.2 Homogenization of the periodic thermal conduction problem and approximation method

Class of integral equations based on polarization

In the framework of the homogenization theory, we consider the thermal conduction problem of the unit cell V = a × a × a, governed by the system of equations (The temperature gradient is defined as e(x) = -∇T where T is the temperature. Introducing temperature is not necessary since it does not affect the solution of our problem).

∇ × e = 0, ∇ • j = 0, j = ke (1.1)
where the heat flux j and temperature gradient e and conductivity k are V -periodic functions of coordinate x. For materials constituted of spherical particles (inclusion-i) embedded in a host material (matrix-m), the conductivity k is constant in each phase and equal to

k = k m in matrix, k = k i in inclusion (1.2)
The effective conductivity k e is then determined via the linear relation between the cell average, notation ... V , of flux J and of temperature gradient

E J = k e E, J = j(x) V , E = e(x) V . (1.3)
We note here that although the original setting of our problem is the thermal conduction phenomena, the same formulation (1.1)-(1.3) and the subsequent results can be applied to diffusion, electrostatic, dielectric, magnetic or electromagnetic wave in the quasistatic limit [START_REF] Gw Milton | The theory of composites[END_REF]. For the latter case, all constants and fields present in the above equations are complex and depend on the excitation frequency ω.

Due to the periodicity of the problem, Eq. (1.1) can be solved using the Fourier transform techniques. Defining the two projection operators G and H as

G(ξ) = ξ ⊗ ξ, H(ξ) = I -ξ ⊗ ξ, ξ = ξ ξ , (1.4)
we can rewrite the Fourier transform of the first two equations of (1.1) in the form

H(ξ)e(ξ) = 0, G(ξ)j(ξ) = 0, ∀ξ = 0. (1.5)
These two equations result from the fact that e(ξ) is collinear with the wave vector ξ of norm ξ and j(ξ) is normal to it. In the reciprocal lattice, ξ are the points whose coordinates are given by

ξ i = 2πn i a , n i = 0, ±1, ±2, ..., i = 1, 2, 3. (1.6)
Now adopting the matrix as the reference material, we define the polarization fields p and q j = k m e + p, e = j/k m + q, (1.7)

which are both vanished in the matrix and finite in the inclusion. Substituting (2.28) into (1.5) and accounting for (1.4) yields the relation

e(ξ) = G(ξ)q(ξ), j(ξ) = H(ξ)p(ξ) (1.8)
valid for all ξ = 0. Equivalently in the physical space, we have

e = E + Gq, j = J + Hp (1.9)
where Gq and Hp represent convolution integrals. In the real space, Eq. (1.9) can be written as

e(x) = E + ξ =0 G(ξ)q(ξ)e -iξ.x , j(x) = J + ξ =0 H(ξ)p(ξ)e -iξ.x (1.10)
Making use of the definitions of p and q in relation to e and j and the property p = -k m q, we obtain two dual integral equations in terms of p

p = αχ (-k m E + Gp) , p = βχ (J + Hp) , (1.11)
In (1.11), α and β are coefficients depending uniquely on the material ratio

α = 1 - k i k m , β = 1 - k m k i (1.12)
and χ the characteristic function

χ = 0 in matrix, χ = 1 in inclusion (1.13)
A class of integral equations based on polarization p can be generated by linearly combining the two equations of (1.11). Using two arbitrary coefficients λ 1 and λ 2 with λ 1 + λ 2 = 1, we obtain the general form

p = χA + χBp (1.14) in which A = (-λ 1 αk m + λ 2 βk e )E, B = λ 1 αG + λ 2 βH. (1.15)
As long as the operator norm B < 1, we can apply successively the relation (1.14) to derive the exact Neumann series expansion of

p p = χA + χBp = χA + χBχA + χBχBp = ... = χ ∞ 0 (Bχ) n A (1.16)
This series presents some features that are original. Like in several works using the conductivity tensor of the matrix as reference materials [START_REF] Nemat-Nasser | Micromechanics : overall properties of heterogeneous materials[END_REF][START_REF] Nemat-Nasser | On composites with periodic structure[END_REF], this series involves only the characteristic function χ. However, unlike in these previous works, the series can always be chosen to be convergent, as shown thereafter. Here, we note that the operator (Bχ) acting on φ yields

(Bχ)φ = ξ e iξ•x B(ξ) ξ ′ χ(ξ -ξ ′ )φ(ξ ′ ) (1.17)
while the Fourier transform of the constant A admits the simple form

A(ξ = 0) = 0, A(0) = A (1.18)
Using (1.17, 1.18), one can write explicitly all the terms in the series (1.16). Regarding the convergence rate, the optimal series can be obtained by minimizing B (hence Bχ ) in (1.15) with λ 1 and λ 2 (see Appendix A). In this case, we have

(λ 1 , λ 2 ) = β β -α , -α β -α , B = αβ β -α = k i -k m k i + k m (1.19)
and

A = αβ α -β (k e -k m )E, B = αβ α -β (I -2G). (1.20)
Certainly, by setting (λ 1 , λ 2 ) = (1, 0) and (λ 1 , λ 2 ) = (0, 1), we recover the two original equations in (1.11). However, the two associated Neumann series, called respectively Eseries (ES) and J -series (JS) can diverge if |α| > 1 or |β| > 1. In contrary, the optimized series (OS) associated to (1.20) having operator B in a similar form as Eyre and Milton's accelerated scheme [START_REF] Dj Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF][START_REF] Gw Milton | The theory of composites[END_REF], always converge with faster rate. These arguments suggests that estimations based on OS could yield better results than those based on ES or JS [START_REF] Qd To | Closed-form solutions for the effective conductivity of two-phase periodic composites with spherical inclusions[END_REF].

Estimation based on Neumann series expansion and integral equation

To determine the effective property k e , we need to find the average polarization p over the inclusion domain Ω = f V . Indeed, from Eq. (2.28), we have the relation

p Ω = f -1 (k e -k m )E (1.21)
Next, from the solution (1.16) and (1.15), we have

p Ω = (λ 2 βk e -λ 1 αk m ) ∞ j=0 C j E (1.22)
With the effective behavior being isotropic and f being the volume fraction, combining (1.21) and (1.22) yields

k e -k m λ 2 βk e -λ 1 αk m = f 3 ∞ j=0 tr(C j ) (1.23)
The explicit expressions of C n are given below

C 0 = I, C 1 = f -1 ξ χ(-ξ)B(ξ)χ(ξ), C 2 = f -1 ξ χ(-ξ)B(ξ) ξ ′ χ(ξ -ξ ′ )B(ξ ′ )χ(ξ ′ ), C j = f -1 ξ 1 ,ξ 2 ,...,ξ n χ(-ξ 1 )χ(ξ 1 -ξ 2 )...χ(ξ j-1 -ξ j )χ(ξ 1 )B(ξ 1 )...B(ξ j ), (1.24)
In numerical practice, we can only kept a finite number n of leading terms of the infinite series, or

k e -k m λ 2 βk e -λ 1 αk m = f 3 n j=0 tr(C j ) (1.25)
On the other hand, the average p Ω can be estimated by making approximation directly to the integral equation at order n

p = χ n-1 j=0 (Bχ) j A + (χB) n p (1.26)
This equation is obtained by repeating the recurrence relation at step n as shown in (1.16). By replacing p with its inclusion phase average p Ω p ≃ χ p Ω (1.27) in the last integral term of (1.26), we obtain an equation for p

Ω p Ω ≃ (λ 2 βk e -λ 1 αk m ) n-1 j=0 C j E + C n p Ω (1.28)
Different from the series truncation approximation, the present approximation includes the information of the remainder of the series and is expected to provide better results [START_REF] Qd To | A numerical-analytical coupling computational method for homogenization of effective thermal conductivity of periodic composites[END_REF]. Consequently, the final equation for the effective conductivity is the following

k e -k m λ 2 βk e -λ 1 αk m = f 3 tr (I -C n ) -1 n-1 j=0 C j (1.29)
In the numerical applications, we shall only integral equation approximation since it provides better results than the series truncation scheme.

Cubic and random distribution of spherical particles

We now assume that our unit cell V = a 3 contains N identical non-overlapping spherical particles of radius R and volume V s = 4πR 3 /3. The shape functions χ(ξ) become

χ(ξ) = V s V F (ξ)ρ(ξ), ρ(ξ) = N i=1 e -iξ.x i , F (ξ) = 3 sin η -η cos η η 3 , η = ξR, (1.30)
where x i is the location of the sphere numbered i in the cell. In statistical mechanics, the term F (ξ) correspond to the form factor of the spherical particle and ρ the microscopic density

ρ(x) = N i=1 δ(x -x i ). (1.31)
For cubic arrangements of sphere, the explicit formulae of χ(ξ) are known analytically.

-Simple Cubic (SC)

χ(ξ) = 3f sin η -η cos η η 3 (1.32)
-Body Centered Cubic (BCC)

χ(ξ) = 3f 2 sin η -η cos η η 3 [1 + (-1) n 1 +n 2 +n 3 ] (1.33)
-Face Centered Cubic (FCC)

χ(ξ) = 3f 4 sin η -η cos η η 3 [(-1) n 1 + (-1) n 2 + (-1) n 3 + (-1) n 1 +n 2 +n 3 ] (1.34)
These expressions allow us to compute numerically C n at any order n. It is also possible to derive analytical expressions for the first two orders C 0 = I and C 1 which is a lattice sum. For C 1 , it can be evaluated by keeping some leading fluctuating terms near the origin and approximate the remainder with a continuous integral [START_REF] Qd To | Closed-form solutions for the effective conductivity of two-phase periodic composites with spherical inclusions[END_REF]. Finally, we obtain a closed-form expression which is accurate for finite fraction f and exact at the infinite volume limit f = 0. However, for higher order n > 1, only numerical solution is available.

For random distribution (RD), at infinite volume limit (both N and V → ∞) and under ergodicity hypothesis, the tensors C j are identical to their ensemble average, notation ... . As a result, we obtain a statistical relations .35) with S (2) , S (3) being the structure factor and the triplet structure factors

C 1 = V s V ξ F (ξ)F (-ξ)B(ξ)S (2) (ξ) C 2 = V 2 s V 2 ξ,ξ ′ F (ξ)F (ξ ′ )F (ξ -ξ ′ )B(ξ)B(ξ ′ )S (3) (-ξ, ξ ′ ), etc... ( 1
S (2) (ξ) = 1 N ρ(ξ)ρ(-ξ) , S (3) (ξ, ξ ′ ) = 1 N ρ(ξ)ρ(ξ ′ )ρ(-ξ -ξ ′ ) , etc... (1.36)
In the present work, we are not interested in higher order structure factors which are complicated and hard to determine in practice. Furthermore, in many numerical applications, using only S (2) (ξ) and S (3) (ξ) have already yielded satisfactory results. Posing η = Rξ and η ′ = Rξ ′ , and considering that for large domains we can rewrite the lattice sums (2.93) as integrals in the Fourier space, we obtain :

C 1 = 1 6π 2 F (ξ)F (-ξ)B(ξ)S (2) (ξ)dη C 2 = 1 36π 4 F (ξ)F (ξ ′ )F (ξ -ξ ′ )B(ξ)B(ξ ′ )S (3) (-ξ, ξ ′ )dηdη ′ . (1.37)
For isotropic distribution of hard spheres, we have the following property

1 6π 2 F 2 (ξ)S (2) (ξ)dη = (1 -f ) (1.38)
regardless the local pair distribution S (2) . As a result, we can evaluate directly C 1

C 1 = 1 3 (2λ 2 β + λ 1 α)(1 -f )I (1.39)
With this expression, the first order approximation (1.29) with (n = 1) based on ES, JS or OS all coincide with the Clausius-Mossoti equation. Regarding S (3) or the three particle distribution g (3) associated to C 2 , due to the lack of accurate analytical solution, it is usually approximated from statistical information at lower levels, for example the radial distribution function g (2) or S (2) . For hard sphere system in equilibrium, both g (2) and S (2) are known from the solution of Ornstein-Zernike equation [START_REF] Ls Ornstein | Accidental deviations of density and opalescence at the critical point of a single substance[END_REF] and Percus-Yevick closure approximation [START_REF] Percus | Analysis of classical statistical mechanics by means of collective coordinates[END_REF][START_REF] Ms Wertheim | Exact solution of the Percus-Yevick integral equation for hard spheres[END_REF]. In literature, there are two popular approximations. The first one is the Kirkwood superposition approximation (KSA) [START_REF] Kirkwood | Statistical mechanics of fluid mixtures[END_REF][START_REF] Torquato | Random heterogeneous materials : microstructure and macroscopic properties[END_REF][START_REF] Hansen | Theory of simple liquids[END_REF].

g (3) (x, x ′ ) ≃ g (2) (x)g (2) (x ′ )g (2) (x ′ -x), (1.40)
and the second one for triplet structure factor [26]

S (3) (ξ, ξ ′ ) ≃ S (2) (ξ)S (2) (ξ ′ )S (2) (ξ -ξ ′ ). (1.41)
Despite the fact that the second approximation provide direct approximation of S (3) (ξ, ξ ′ ), it violates the core condition, i.e the spheres can overlap [START_REF] Hansen | Theory of simple liquids[END_REF]. This can result a bad estimation of the overall conductivity. The KSA scheme does not suffer from this limitation but it is formulated in physical space and must be translated to ξ space. From the link between S (3) and g (3) and the KSA scheme (see Appendix B), we propose the following alternative approximation

S (3) (ξ, ξ ′ ) ≃ S (2) (ξ)S (2) (ξ ′ ). (1.42)
1.2.4 Integral evaluation of C 2 for random distribution with S (3) approximation

Regarding C 2 , it is impossible to derive such a simple analytical formula and independent of S (3) . However, using approximation (1.42), we can convert the six-fold integral

C 2 = 1 36π 4 T (η)T (η ′ )F (|η -η ′ |)B(η)B(η ′ )dηdη ′ , (1.43)
with T (η) = F (ξ)S (2) (ξ) into a two-fold integral involving η and η ′ . Indeed, we fix first η, η ′ and l = |ηη ′ | and integrate over a circle C ϕ of radius R ϕ = η ′ sin ϕ (see Fig. 1.1). By posing m 1 = λ 1 α and m 2 = λ 2 β, we have

Cϕ B(η ′ )ds = πm 1 R ϕ [2 cos 2 ϕG(η) + sin 2 ϕH(η)] +πm 2 R ϕ [(2 -sin 2 ϕ)H(η) + 2 sin 2 ϕG(η)] (1.44)
Now, multiplying B(η) with (1.44) yields

B(η) Cϕ B(η ′ )ds = 2m 1 πR ϕ [m 1 cos 2 ϕ + m 2 sin 2 ϕ]G(η) + m 2 πR ϕ [m 2 (2 -sin 2 ϕ) + m 1 sin 2 ϕ]H(η) (1.45)
Next, integrating with η ′ running over the sphere surface S η ′ of radius η ′ and then with η running over the sphere surface S η of radius η, we obtain it is now ready to rewrite C 2 in a simple form

S η ′ Sη F (η -η ′ )B(η)B(η ′ )dSdS ′ = = 2πm 1 η ′2 Sη G(η)dS π 0 F (l)[m 2 sin 2 ϕ + m 1 cos 2 ϕ] sin ϕdϕ +πm 2 η ′2 Sη H(η)dS π 0 F (l)[m 2 (2 -sin 2 ϕ) + m 1 sin 2 ϕ] sin ϕdϕ, l = η 2 + η ′2 -2ηη ′ cos ϕ (1.46)
C 2 = 2I 27π 2 T (η)T (η ′ )R(η, η ′ )η 2 η ′2 dηdη ′ , (1.48)
with the function R(η, η ′ ) being defined by

R(η, η ′ ) = 1 -1 F ( η 2 + η ′2 -2ηη ′ t)[m 2 2 (1 + t 2 ) + 2m 1 m 2 (1 -t 2 ) + m 2 1 t 2 ]dt. (1.49)
Numerical examples show that the second order approximation (1.29) with (n = 2) and (1.43) based on OS provides improved result with respect to the CM estimates.

Sample generation and direct computation of C 2 for random distribution

As shown above, the calculation of the effective properties based on the integral form encounter difficulties due to the lack of accurate information of S (3) . Although many approximations of S (3) (or its dual functions g (3) , h (3) ) have been proposed in the literature [START_REF] Jl Barrat | Factorization of the triplet direct correlation function in dense fluids[END_REF][START_REF] Ar Denton | High-order direct correlation functions of uniform classical liquids[END_REF], each approximation scheme is subject to its own limitation and can only capture certain aspects of S (3) . In this work, it has been to determine directly C 2 from the ensemble average of real samples. To this aim, we use the ensemble average version of Eq.1.50.

C 2 = f -1 ξ χ(-ξ)B(ξ) ξ ′ χ(ξ -ξ ′ )B(ξ ′ )χ(ξ ′ ) (1.50)
Here, we use the expression Eqs. (2.93,2.94) and calculate the discrete convolution in Fourier space. The Event Driven Molecular Dynamics (EDMD) method [START_REF] Dc Rapaport | The Art of Molecular Dynamics Simulation[END_REF] is used to generate systems of 500 hard spheres in equilibrium at the unit temperature. The final results are then obtain by averaging over 50 samples and will be compared with other estimations. As an example of intermediate result, the radial distribution function (RDF) of the system with a snapshot of a typical sample is shown in Fig. 1.2.

Cubic lattice arrangements

We apply first the approximation schemes to cubic lattice arrangements of spheres, simple cubic (SC) and face centered cubic (FCC). The reference results for comparison are the exact numerical solutions obtained from a standard FFT numerical method. The details of the latter has been described elsewhere [START_REF] Monchiet | A polarization-based fft iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF][START_REF] Qd To | A numerical-analytical coupling computational method for homogenization of effective thermal conductivity of periodic composites[END_REF] and shall not be discussed further here.

From results for SC arrangement in Table 1.1, we find that the OS approximation performs significantly better at all order n = 1 or 2 and contrast ratio k i /k m when compared with ES and JS approximation. The second order approximation (n = 2) yields better . results than the first order (n = 1) except for ES series for k i /k m = 10, 1000 and JS series for k i /k m = 0.1, 0.001 (not shown in Table 1.1). This could be explained from the fact that the associated series diverge, making higher order approximation meaningless. For k i /k m = 0.1, 0.001, there is little difference between the approximation scheme and the exact results. At very high contrast ratio k i /k m = 1000 and very high volume fraction, the series approximation work less well. The best OS approximation is only acceptable at f ≤ 0.4. It means that higher order approximation n > 2 is needed to achieve the better precision, for example within 5% error.

Figure 1.3 for FCC arrangement shows that the OS approximation still perform very well for n as low as 2, k i /k m as high as 100 and volume fraction upto f = 0.6. While the second order approximation have significantly improved the first order. There is a very slight difference between the third order and the second order. This agrees with the very fast convergence rate of the OS series, i.e the final solution is not far from the sum of several initial terms in the series.

All the results obtained for cubic lattice arrangement have confirmed our belief on the performance of the coupling between the OS series and integral equation approximation.

We shall proceed with the random distribution of sphere in the next section. The approximation based on triplet structure factors and Neumann series is now applied to determine the effective conductivity for a randomly distribution of spheres in equilibrium. These results are then compared with previous works on the same system in the literature. The first is the Torquato expression [START_REF] Torquato | Effective electrical conductivity of two ?phase disordered composite media[END_REF] (TO) based on another expansion [START_REF] Wf Brown | Solid mixture permittivities[END_REF] which gives rise the microstructure parameter ζ. The latter is intrinsically linked to three point correlation function and can be evaluated with KSA approximation (1.40) or by simulation [START_REF] Miller | Effective conductivity of hard-sphere dispersions[END_REF]. As mentioned earlier, this series expansion can yield good result by keeping two or three leading terms but does not guarantee the convergence of the whole series, especially near the percolation limit. Regarding numerical simulations, Bonnecaze and Brady [START_REF] Rt Bonnecaze | A method for determining the effective conductivity of dispersions of particles[END_REF][START_REF] Rt Bonnecaze | The effective conductivity of random suspensions of spherical particles[END_REF] (BB) using multipole expansion of particle potential and Kim and Torquato [START_REF] Ic | Effective conductivity of suspensions of hard spheres by brownian motion simulation[END_REF] (KT) using random walk method have obtained the effective conductivity of the system.

It is clear that our OS-1 approximation coincides with the Clausius-Mossoti expression. This is a good starting point since one can expect higher order approximation will improve the latter. From results given in Table 1.2, we find that our OS-2 estimation with or without S (3) aproximation is in very good agreement with the numerical results issued from previous works at moderate contrast k i /k m < 10 for whole volume fraction f . For infinite contrast k i /k m = ∞, significant deviation can be observed at f = 0.5 for OS-2 with S (3) approximation. This can be explained from the fact that, near the jamming state f = 0.5, the difference between the theoretical PY solution and the RDF curves have been already significant [START_REF] Rt Bonnecaze | A method for determining the effective conductivity of dispersions of particles[END_REF][START_REF] Rt Bonnecaze | The effective conductivity of random suspensions of spherical particles[END_REF] while the OS-2 estimation with S (3) aproximation relies heavily on the PY solution. In this case, the direct computation based on the real samples should be used. Indeed, numerical results have shown that the direct calculation is in very good agreement with the literature result even for the case of infinite contrast and volume fraction as large as f = 0.5.

Final remarks and discussion

To summarize, we have developed an estimation scheme for conductivity of composite materials made of spherical particles and matrix. To do that, we derive first the optimal integral equation of polarization p from a family of integral equation of the same type, i.e the corresponding Neumann series must converge the fastest. Then, a new integral equation is constructed including the first n -1 terms of the series and an integral operator of power n acting on the polarization. By applying approximation to the latter, one can expect a good estimation of the average polarization and hence the effective conductivity. Indeed, numerical examples have provided supporting evidence for this conclusion.

Another interesting result of the present work concerning random media is that we have established the statistical connection with the structure factors, the triplet structure factors and higher order structure factors etc... These quantities are very important since they provide useful information of local distribution of particles. For the case of colloidal 3) approximation) and Torquato's expression [START_REF] Torquato | Effective electrical conductivity of two ?phase disordered composite media[END_REF] (TO), Kim and Torquato [START_REF] Ic | Effective conductivity of suspensions of hard spheres by brownian motion simulation[END_REF] (KT), the simulation data of Bonnecaze and Brady [START_REF] Rt Bonnecaze | The effective conductivity of random suspensions of spherical particles[END_REF] (BB). The Clausius-Mossoti (CM) approximation coincides with first order OS approximation OS-1 (n = 1). suspensions, they can be obtained experimentally via the scattering techniques [START_REF] Borsali | Soft-Matter Characterization[END_REF] or via particle simulation methods [1]. Such dual formulation in wave vector have thus given an alternative way to the use of correlation function [START_REF] Torquato | Random heterogeneous materials : microstructure and macroscopic properties[END_REF] in order to characterize the properties of the random media.

Based on the estimation of the remainder of the Neumann series that is closely related to FFT numerical method, the present method improves the limitations of the effective medium approximation. It is known that the latter fails when the inclusion strongly interact especially at high inclusion volume fraction and high contrast ratio. Without limiting to spherical inclusions, the present method is general enough to deal with complex microstructures, for example, those constituted of cylindrical, ellipsoidal shape inclusions, which will be subject to our future study.

ACKNOWLEDGMENTS We thank Mr. Van-Huyen Vu for his help to construct the random hard sphere samples.

Appendix Appendix A : Operator norm minimization

The norm of the self-adjoint operator B is the maximal eigenvalue λ satisfying

Bv = λv, (1.51)
Making use of (1.15), we have the following relations in Fourier space

[λ 1 α -λ] G(ξ)v(ξ) + [λ 2 β -λ] H(ξ)v(ξ) = 0 (1.52)
Since G(ξ) and H(ξ) are projection operator, there are two possibilities

G(ξ)v(ξ) = 0, λ = λ 2 β or H(ξ)v(ξ) = 0, λ = λ 1 α (1.53)
As a result, the operator norm B admit the expression

B = max {|λ 1 α| , |λ 2 β|} (1.54)
Optimizing B with λ 1 and λ 2 and accounting for λ 1 + λ 2 = 1, we have

λ 1 = β β -α , λ 2 = -α β -α , B = |αβ| |β -α| (1.55)
Regarding the operator Bχ in the Neumann series, we can estimate its norm via the inequality from the previous result

Bχv ≤ B χv ≤ B χ v ≤ B v , (1.56)
From this sequence of inequalities, we conclude that

Bχ ≤ B = |αβ| |β -α| = k i -k m k i + k m < 1, (1.57)
which ensures the convergence of the series.

Appendix

Appendix B : Triplet structure factor approximation

Using the definition of the triplet structure factor and after some mathematical manipulations, we obtain the expression for S (3)

(-ξ, ξ ′ ) S (3) (-ξ, ξ ′ ) = 1 N i e iξ.x i j e -iξ ′ .x j k e -i(ξ-ξ ′ ).x k = 1 + 1 N i =k e iξ.(x i -x k ) + 1 N j =k e -iξ ′ .(x j -x k ) + 1 N i =j =k e iξ.(x i -x k ) e -iξ ′ .(x j -x k ) , (1.58)
or equivalently

S (3) (-ξ, ξ ′ ) = S(ξ) + S(ξ ′ ) -1 + 1 N i =j =k e iξ.(x i -x k ) e -iξ ′ .(x j -x k ) (1.59)
Here, x i , x j and x k are respectively the locations of the particles i, j and k. We shall now continue with the last term in (1.59) which is equal to

1 N e iξ.(x-x ′′ ) e -iξ ′ .(x ′ -x ′′ ) i =j =k δ(x -x i )δ(x ′ -x j )δ(x ′′ -x k ) dxdx ′ dx ′′ = ρ3 N e iξx e -iξ ′ x ′ g (3) (x + x ′′ , x ′ + x ′′ , x ′′ )dxdx ′ dx ′′ = ρ2 e iξx e -iξ ′ x ′ g (3) (x, x ′ )dxdx ′ ≃ ρ2 e iξ.x e -iξ ′ .x ′ g(x -x ′ )g(x)g(x ′ )dxdx ′ , (1.60)
with ρ = ρ(x) = N/V being the average particle density. To derive the last equation of (1.60), we have applied the KSA scheme to g (3) . This is the exact expression of S (3) (-ξ, ξ ′ ) based on KSA. However, this expression is still difficult to use since it is a double integral in 3D space. For simplification, we approximately replace g(xx ′ ) with its average 1 and obtain

ρ2 e iξ.x e -iξ ′ .x ′ g(x -x ′ )g(x)g(x ′ )dxdx ′ ≃ ρ2 e iξ.x e -iξ ′ .x ′ g(x)g(x ′ )dxdx ′ = (S(ξ) -1)(S(ξ ′ ) -1). (1.61)
Finally, the triplet structure factor S (3) (-ξ, ξ ′ ) can be recast in the approximate simple form (1.42).

Chapitre 2

Overall elastic properties of composites from optimal strong contrast expansion

Quy-Dong To

Abstract

In this paper, we propose a new systematic procedure of estimating elastic properties of composites constituted of two phases, matrix and inclusions. A class of integral equations based on eigenstrain (or eigenstress) with the matrix as reference material is constructed with an explicit form in Fourier space. Each integral equation belonging to this class can yield estimates of the overall elastic tensor via Neumann series expansion. The best estimates and series are selected based on the convergence rate criteria of the series, i.e the spectral radius must be minimized. The optimized series is convergent for any finite contrast between inclusions and matrix. Applying the optimized series and the associated estimates to different microstructures yields very satisfying results when compared with the related full solution. For the case of a random distribution of spherical inclusions, exact relations between the elastic tensor and n-th order structure factors are demonstrated.

Introduction

We consider the problem of finding the effective stiffness tensor C e of periodic heterogeneous matrix-inclusion materials. The existence of C e and homogenization procedure have been rigorously founded in the literature [START_REF] Sanchez-Palencia | Non-homogeneous media and vibration theory[END_REF]. Given the distribution of the constituents, the cell problem must be solved first and the linear relation between average stress and strain is then established. Estimates can be obtained by making relevant approximation to the ingredients constituting the effective tensor [see e.g] [START_REF] Jd Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF][START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF]. Although the present contribution concerns the theory of optimally estimating C e from the microstructure, it is closely related to FFT numerical homogenization methods.

By introducing a reference material C 0 , the heterogeneity effect can be viewed as a distribution of eigenstrains within an homogeneous material. Using the related Green tensor, our problem can be formulated as a Lippmann-Schwinger equation for eigenstrain [START_REF] Nemat-Nasser | On composites with periodic structure[END_REF]. The integral equation is the origin of resolution methods based on iteration and Fast Fourier Transform (FFT) techniques [START_REF] Michel | Effective properties of composite materials with periodic microstructure : a computational approach[END_REF][START_REF] Bhattacharya | A model problem concerning recoverable strains of shape-memory polycrystals[END_REF]. Significant progresses have been made regarding the improvement of convergence rate [START_REF] Michel | Effective properties of composite materials with periodic microstructure : a computational approach[END_REF][START_REF] Dj Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF][START_REF] Gw Milton | The theory of composites[END_REF][START_REF] Monchiet | A polarization-based fft iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF][START_REF] Brisard | Fft-based methods for the mechanics of composites : A general variational framework[END_REF]. The study of convergence rate in those works will be extended in the present contribution in the case of new integral equations.

We can remark that the iteration scheme used to solve the Lippmann-Schwinger equation corresponds to the Neumann series summation. The latter can be used to derive exact theoretical relations and estimates [see e.g and the references therein] [74, 42], using for example the weak and strong contrast expansions. As an extension to previous works on conductivity, [START_REF] Torquato | Effective stiffness tensor of composite media : II. Applications to isotropic dispersions[END_REF][START_REF] Torquato | Effective stiffness tensor of composite media. I. Exact series expansions[END_REF] introduced the cavity strain field and derived the strong contrast estimates. Connections with statistical information related to the distribution of two phases, the n-point correlation function, were also demonstrated in this context. At the second order, the estimates produce satisfying results for systems composed of dispersed non overlapping spheres. However, the series of Torquato is increasingly complicated at higher order without guarantee of convergence. As already shown in [START_REF] Gw Milton | The theory of composites[END_REF], the starting point of the strong contrast expansion is a conditionally convergent series. In a similar way, bounds of the effective elasticity tensor are also depending on correlation functions [see pps 313, 554] [START_REF] Gw Milton | The theory of composites[END_REF] and in this context, the use of the matrix as a reference material simplifies significantly the expression of these bounds In this paper, we propose a new estimate based on series expansion that works at all finite contrast ratio, while using the matrix as a reference material. Additionally, we can control and optimize the convergence rate so that the series converges in the quickest way, and therefore produces the best estimates when using a finite sum in the series expansion. A class of integral equations for eigenstrain depending on two parameters α, β is first derived. The spectral radius and norm of the corresponding operators are bounded by analytical expressions. Different optimization methods are proposed to find the fastest series convergence and the associated estimates. It is noteworthy that a convergence will be obtained for any finite contrast between matrix and inclusions, while using always the matrix as a reference material, which simplifies significantly the series expansion.

Similarly to the estimations of the effective elasticity tensor using correlation functions [START_REF] Gw Milton | The theory of composites[END_REF][START_REF] Torquato | Random heterogeneous materials : microstructure and macroscopic properties[END_REF], the new method presented in this paper allows to estimate the effective elasticity tensor using the n-order structure factors, which represent the counterpart in Fourier space of correlation functions. As an example, a direct connection of the effective elasticity tensor to n-th order structure factors is given in the case of randomly distributed spheres. Numerical applications for cubic arrays and random distribution of spheres yield very good results in comparison with FFT based methods and other results from the literature. The details of those contributions are presented in the following.

Mathematical preliminaries

As it will be seen in the following sections, the problem of the convergence of the Neumann series associated to a given integral equation is related to the spectrum of a fourth order iteration tensor which appears in the integral equation. In this section, the notations and mathematical results related to these fourth order tensors will be presented.

Notations and definitions

Before proceeding to the considerations on the spectrum of fourth order operators, we introduce first the system of abstracted notations and some important definitions which will be used throughout the paper. Most of our calculations involve symmetric second order tensors and fourth order tensors with minor symmetries. Unless specified, two tensors standing next to each other implies their double contraction product. Given two fourth order tensors A and B and two second order tensors u and v, the double contraction products AB, Au and uv read

(AB) ijkl = A ijmn B nmkl , (Au) ij = A ijkl u lk , uv = u ij v ji , i, j, k, l = 1, 2, 3 (2.1)
The symbol ⊗ denotes the usual tensorial product and ⊗ a special symmetrized tensorial product between two second order tensors u and v

(u⊗v) ijkl = 1 2 (u ik v jl + u il v jk ), (u ⊗ v) ijkl = u ij v kl , i, j, k, l = 1, 2, 3. (2.2)
We are also dealing with periodic functions using Fourier analysis. Any V -periodic tensor field u, function of coordinate x(x 1 , x 2 , x 3 ) can be expressed as an infinite Fourier series

u(x) = ξ u(ξ)e iξ.x , (2.3)
with u(ξ) being the Fourier transform of u(x)

u(ξ) = 1 V V u(x)e -iξ.
x dx, ∀ξ.

(2.4)

In (2.3), the sum is taken over all the wavevectors ξ whose components ξ 1 , ξ 2 and ξ 3 are defined from the relation

ξ i = 2πn i l i , i = 1, 2, 3, n 1 , n 2 , n 3 ∈ Z, (2.5)
where l 1 , l 2 , l 3 are the periods along the three directions, i.e the dimensions of the unit cell V . According to the convolution theorem, products in Fourier space correspond to convolution products, (notation * , in physical space) by

[A * u](ξ) = A(ξ)u(ξ), (2.6)
In the paper, we will encounter frequently equations in the form

u = U + A * u (2.7)
for a given second order tensorial function U (x) and fourth order tensorial operator A.

The solution u of the above equation is the following Neumann series

u = ∞ n=0 (A * ) n U . (2.8)
Numerically, the above summation is equivalent to repeating the iterative scheme

u 0 = U , u n+1 = U + A * u n . (2.9)
until convergence. The convergence rate of the Neumann series (or the iterative scheme) can be estimated from the spectral radius or the norm of the associated operator A. The spectral radius of an operator A is the absolute maximum of all eigenvalues λ

ρ(A) = max |λ| : ∃u, A * u = λu, or A(ξ)u(ξ) = λu(ξ) (2.10)
The set of symmetric second order tensors has the structure of an Hilbert space related to the scalar product between two second order tensor fields, notation ., .

u, v = 1 V V u(x)v(x)dx = ξ u(ξ)v(ξ), (2.11)
where the overline notation • stands for the complex conjugate value. Consequently, the related norm, notation . , of any tensor field u(x) can be defined by the expression u = u, u .

(2.12)

The norm of operator A is then defined as

A = inf{c ≥ 0 : A * u ≤ c u } (2.13)
and in the Hilbert space, it can be evaluated as

A = ρ(A † A), A ijkl (ξ) = A † klij (ξ) (2.14)
with A † being the adjoint operator of A. The norm and the spectral radius meet the inequality

ρ(A) ≤ A , (2.15)
where the equality is achieved if A is self-adjoint or at least normal, i.e A † A = AA † .

Walpole base

We consider operators which are transversely isotropic fourth order tensors in Fourier space and have minor symmetries, the plane of isotropy being normal to the wave vector ξ. The related direction vector along ξ is denoted as ξ. It has been shown by [START_REF] Walpole | Elastic behavior of composite materials : theoretical foundations[END_REF] that these tensors are a linear combination of the following six tensors which constitute the Walpole base :

E 1 = 1 2 k ⊥ ⊗ k ⊥ , E 2 = k ⊗ k, E 3 = k ⊥ ⊗k ⊥ -E 1 , E 4 = k ⊥ ⊗k + k⊗k ⊥ , E 5 = k ⊗ k ⊥ , E 6 = k ⊥ ⊗ k, (2.16)
where k and k ⊥ are the second order projection tensors

k = ξ ⊗ ξ, k ⊥ = I -k, (2.17)
and I the second order identity tensor. As shown later, most of our calculations involve transversely isotropic tensors (Green tensors in Fourier space) and constant isotropic tensors (elasticity tensors of isotropic constituents). Those calculations can be considerably simplified using Walpole base elements. For example, the fourth order identity tensor I and the spherical and deviatoric projection tensors J and K can be expressed in the Walpole base in the form

I = E 1 + E 2 + E 3 + E 4 , J = 1 3 I ⊗ I = 1 3 (2E 1 + E 2 + E 5 + E 6 ), K = I -J = 1 3 (E 1 + 2E 2 -E 5 -E 6 ) + E 3 + E 4 . (2.18)
Given any two transversely isotropic tensors A and B in the Walpole base

A = a 1 E 1 + a 2 E 2 + a 3 E 3 + a 4 E 4 + a 5 E 5 + a 6 E 6 , B = b 1 E 1 + b 2 E 2 + b 3 E 3 + b 4 E 4 + b 5 E 5 + b 6 E 6 , (2.19)
their double inner product can be computed via the formula

AB = (a 1 b 1 + 2a 6 b 5 )E 1 + (a 2 b 2 + 2a 5 b 6 )E 2 + a 3 b 3 E 3 + a 4 b 4 E 4 +(a 2 b 5 + a 5 b 1 )E 5 + (a 1 b 6 + a 6 b 2 )E 6 .
(2.20)

These properties are important for the calculation of the spectral radius or norm, which will be presented in the later sections. From the definition of the adjoint operator in the previous section, we can also deduce that

E † i = E i , ∀i = 1, 2, 3, 4, E † 6 = E 5 , E † 5 = E 6 . (2.21)
In other words, the operators E 1 , E 2 , E 3 and E 4 are self-adjoint and the operators E 5 and E 6 are adjoint operators of each other. The fourth order operator A is self-adjoint if its components along E 5 and E 6 are identical.

Spectral radius and norm of operators via Walpole base

Always working in Fourier space, we shall determine the eigenvalues of A(ξ) in the form (2.19). Multiplying both sides of the eigen equation (2.10) with E i i = 1, 2, .., 6, and then with v, we have a system of linear equations

               a 1 (vE 1 v) + a 6 (vE 6 v) = λ(vE 1 v) a 2 (vE 2 v) + a 5 (vE 5 v) = λ(vE 2 v) a 3 (vE 3 v) = λ(vE 3 v) a 4 (vE 4 v) = λ(vE 4 v) a 1 (vE 5 v) + 2a 6 (vE 2 v) = λ(vE 5 v) a 2 (vE 6 v) + 2a 5 (vE 1 v) = λ(vE 6 v) (2.22)
It turns out that the vector of components vE i v with i = 1, 2, ..6 is an eigenvector of the matrix built with the elements a i appearing in the left side of (2.22) and λ is an eigenvalue of this matrix. As a result, λ must be the solution of the eigen equation

a 1 -λ 0 0 0 0 a 6 0 a 2 -λ 0 0 a 5 0 0 0 a 3 -λ 0 0 0 0 0 0 a 4 -λ 0 0 0 2a 6 0 0 a 1 -λ 0 2a 5 0 0 0 0 a 2 -λ = 0 (2.23) or (a 3 -λ)(a 4 -λ)[(a 1 -λ)(a 2 -λ) -2a 5 a 6 ] 2 = 0. (2.24)
Finding the roots of Eq. (2.24) yields the values for λ and the spectral radius of A can now be computed with the formula

ρ(A) = max |a 3 |, |a 4 |, 1 2 |(a 1 + a 2 ) ± (a 1 -a 2 ) 2 + 8a 5 a 6 | . (2.25)
Applying the same procedure as before to find the eigenvalues of A † A, we can determine

A A = ρ(A † A) = max |a 3 |, |a 4 |, 1 2 [ (a 1 -a 2 ) 2 + 2(a 5 + a 6 ) 2 + + (a 1 + a 2 ) 2 + 2(a 5 -a 6 ) 2 ] . (2.26)

Homogenization of elastic periodic composites

The computation of effective properties of composites can be expressed in the real space. However, one important ingredient of the homogenization process is the use of Green's tensors which are highly singular in the real space, while their Fourier transforms can be expressed analytically in a simple way. In addition, it will be shown that the use of the formulation in Fourier space introduces naturally the structure factors which are related to the geometrical distribution of the component materials. As a result, the formulation will be done jointly in real space and in Fourier space. In a first step, the homogenization problem will be expressed by using integral equations. Next, the optimization of the solution of these integral equations by using Neumann series will be studied.

Governing integral equations

We consider a heterogeneous material where the local isotropic stiffness C(x) (compliance S(x)) is a V -periodic function of the coordinates x. The local bulk and shear moduli are denoted respectively as κ(x) and µ(x). To determine the overall elastic properties of the material, we need to solve the following periodic boundary value problem for stress σ and strain ǫ, where the average of one of these fields over V is given. Then the effective elasticity tensor C e can be computed from the linear relation between the average strain E and stress Σ

Σ = C e E, Σ = σ V , E = ǫ V .
(2.27)

Here, we adopt the notation V to refer to the average over volume V of the quantity inside the brackets. Introducing a reference material with an isotropic stiffness C 0 (compliance S 0 ), we have

σ = C 0 ǫ + τ , ǫ = S 0 σ + e, (2.28)
where the eigenstress τ and the eigenstrain e have been introduced. By usual ways [START_REF] Gw Milton | The theory of composites[END_REF][START_REF] Nemat-Nasser | On composites with periodic structure[END_REF], the integral equation for the eigenstress τ may be expressed as :

τ = δ ′ C(E -Γ 0 * τ ), δ ′ C(x) = C(x) -C 0 (2.29)
and the dual integral equation for the eigenstrain e e = δ ′ S(Σ -∆ 0 * e), δ ′ S(x) = S(x) -S 0 .

(2.30)

The Green operators Γ 0 and ∆ 0 for strain and stress are defined in Fourier space by :

Γ 0 (ξ) = 3 3κ 0 + 4µ 0 E 2 + 1 2µ 0 E 4 ∀ξ = 0, Γ 0 (0) = 0, ∆ 0 (ξ) = 18µ 0 κ 0 3κ 0 + 4µ 0 E 1 + 2µ 0 E 3 ∀ξ = 0, ∆ 0 (0) = 0, (2.31)
where ξ is the wave vector and the tensors E i are defined in Fourier space from the Walpole base by using ξ as the unit vector in the direction of ξ. The elastic constants κ 0 and µ 0 appearing in (2.31) are respectively the bulk modulus and the shear modulus associated to the reference tensor C 0 . It is clear that the two equations (2.29,2.30) are equivalent to strain and stress integral equations which Fast Fourier Transform (FFT) numerical homogenization methods are essentially based on (see Appendix A). They also have the same form as (2.7) and can be solved using Neumann series expansion (or iterative scheme techniques) based on (2.8,2.9).

In general, the solutions for strain and stress are unique and independent from the choice of the reference medium. However, the reference material can affect the convergence if Neumann series expansion is used, especially when one of the constituents is chosen as a reference material. From another point of view, the choice of one of the constituents as a reference material simplifies drastically the form of the Neumann series associated to the integral equation. Connections to important statistical information like structure factors appear directly and naturally using this formulation, as it will be shown thereafter. Therefore, instead of using stress or strain formulation, our strategy is to derive a unified formulation by exploiting the relations between stress σ, strain ǫ and eigen-stress τ (or eigen-strain e). The two latter quantities are connected via the relation τ = -C 0 e. It will be shown thereafter that this strategy is compatible with the choice of one of the constituents as a reference material.

From the two elementary integral equations, we can construct a family of integral equations for τ (or e equivalently) by linear combination. Using two tensors L(x) and I-L(x), we can obtain

τ = A ′ E + B ′ * τ , (2.32)
in which

A ′ = ((I -L)(δ ′ C) -LC 0 (δ ′ S)C e ) B ′ = -((I -L)(δ ′ C)S 0 P + LC 0 (δ ′ S)Q). (2.33)
Here, the two projection operators P and Q are defined again in Fourier space by :

P = C 0 Γ 0 = E 2 + E 4 + 3κ 0 -2µ 0 3κ 0 + 4µ 0 E 6 , Q = ∆ 0 S 0 = E 1 + E 3 - 3κ 0 -2µ 0 3κ 0 + 4µ 0 E 6 , (2.34)
Equation (2.32) can be used to solve τ with an iteration scheme corresponding to the associated Neumann series as long as ρ(B ′ ) < 1. This condition guarantees that the Neumann series, which is obtained by repeating the recurrence (2.32)

τ = A ′ E + B ′ * τ = A ′ E + B ′ * A ′ E + B ′ * B ′ * τ = ... = ∞ n=0 (B ′ * ) n A ′ E (2.35)
will surely converge. The convergence rate depends on the magnitude of ρ(B ′ ). The smaller it is, the faster the series converge. For convenience, let us call the series associated to (2.29) as E-series (ES) and (2.30) as Σ-series (SS). Those are particular cases of (2.32) where L = 0 and L = I. The series in the form (2.32) with suitable L that yields the fastest convergence rate will be called as optimized series and studied in the following.

We consider a two-phase material composed of inclusions embedded in the matrix. At this stage, it can be remarked that the choice of the reference material is crucial in the expression of the power of B ′ in (2.3.1). Indeed, if the reference material is not chosen properly, the operator B ′ will contain simultaneously the characteristic functions of both constituents defined thereafter and the expression of (B ′ * ) n will become increasingly complicate. Therefore, in the following, the reference material is chosen to be the matrix C 0 and the stiffness of the inclusion is denoted as C 1 . The tensor L is assumed to be constant and isotropic with the representation

L = 2αK + 3βJ (2.36)
in which α and β are two constants. Like L, the tensor C i with i = 0, 1 can be expressed in a similar way using the bulk stiffness κ i and the shear stiffness µ i as parameters.

Optimization based on the Green operators

In this section, we employ the method developed in [START_REF] Mt Nguyen | Conductivity estimates of sphericalparticle suspensions based on triplet structure factors[END_REF][START_REF] Gw Milton | The theory of composites[END_REF] for heat conduction to estimate the convergence rate of the general series (2.32) in elasticity. This method will require the mathematical preliminaries that we presented previously. Adopting the notation χ for the characteristic function,

χ(x) = 0 in matrix, χ(x) = 1 in inclusion (2.37)
and recalling that the reference tensor is the elasticity tensor of the matrix, we can rewrite the Neumann series in the following form

τ = ∞ j=0 (B ′ ) * j A ′ E (2.

38)

where

A ′ = χA, B ′ = χB and A = ((I -L)(δC) -LC 0 (δS)C e ) B = -((I -L)(δC)S 0 P + LC 0 (δS)Q) δC = C 1 -C 0 , δS = S 1 -S 0 (2.39)
It is noteworthy that the choice of the matrix as a reference material allows us to obtain an expression involving only χ in its expression and power series, while another choice would contain two characteristic functions of the two materials at the same time.

Another strategy behind this expression is also that it is easier to quantify separately the spectral radius and norm of B and χ, which in turn provides information on B ′ . The function χ takes only 0 and 1 values and depends on the distribution of the inclusions in the cell. Its spectral radius and norm can be obtained straightforwardly, for example

ρ(χ) = χ = 1. (2.40)
Intuitively, to optimize the convergence rate of B ′ , we need to optimize ρ(B) or the norm B . Since B is generally not self-adjoint, optimizing B and ρ(B) with L may lead to different results (see Figs. 2.1,2.2). For the sake of clarity, we shall use the two abbreviated names : OR for optimization by spectral radius and ON for optimization by norm.

As an example, let us consider the case where the two materials are isotropic with the same Poisson ratio ν, or

C 1 = εC 0 , (2.41)
where ε is the contrast ratio. We are also limited to the case where

L = 2αI, (2.42)
and minimize ρ(B). First, B can be written explicitly as

B = -(ε -1)(1 -2α)(E 2 + E 4 ) + 2α (ε -1) ε (E 1 + E 3 ) - - 3κ 0 -2µ 0 3κ 0 + 4µ 0 E 6 (1 -2α)(ε -1) + 2α (ε -1) ε (2.43)
It is clear that the spectral radius of B can be evaluated as

ρ(B) = max |(ε -1)(1 -2α)|, 2α (ε -1) ε . (2.44)
Since both eigenvalues are linear function of α, optimizing ρ(B) yields the following results

2α = ε ε + 1 , ρ(B) = ε -1 ε + 1 , B = ε -1 ε + 1 E 1 -E 2 + E 3 -E 4 - 2(3κ 0 -2µ 0 ) 3κ 0 + 4µ 0 E 6 . (2.45)
Comparing with two original series, for example the ES series

B = -(ε -1) E 2 + E 4 + 3κ 0 -2µ 0 3κ 0 + 4µ 0 E 6 , ρ(B) = |ε -1|, (2.46)
and the SS series

B = - 1 -ε ε E 1 + E 3 - 3κ 0 -2µ 0 3κ 0 + 4µ 0 E 6 , ρ(B) = ε -1 ε , (2.47)
the series after optimization shows a clear advantage. The spectral radius of the latter is the smallest and less than 1. The ES series can diverge if ε > 2 and the SS series can diverge if ε < 1/2.

Minimizing ρ(B) and B for the general case is not analytically simple, but it can be done numerically without difficulties. The main problem is that we need to ensure the convergence related to operator ρ(B ′ ). However, optimizing ρ(B) only guarantees the convergence of the series associated to B not B ′ while optimizing B may overestimate ρ(B ′ ). This is due to the inequalities As a result, in the case where B is not self-adjoint, those estimations may not be strict enough. To overcome those issues, we shall present now a more direct optimization method (OD) to estimate ρ(B ′ ), detailed in the next subsection.

ρ(B ′ ) ≤ B ′ , ρ(B) ≤ B , B ′ ≤ χ B = B . ( 2

Direct estimation of ρ(B ′ ) and optimization

We introduce now an approach which is closely related to the one presented by [START_REF] Michel | Effective properties of composite materials with periodic microstructure : a computational approach[END_REF] to estimate the convergence rate of the series. We begin first with the ES series whose eigen-equation reads

-(C -C 0 )S 0 P * τ = λτ .
(2.49)

Equivalently, we can rewrite the above equation as

-C ′ -C 0 S 0 P * τ = τ , C ′ = 1 λ C - (1 -λ) λ C 0 . (2.50)
Comparing with (2.29), we conclude that this is the equation for the boundary value problem with zero strain E = 0 and local stiffness C ′ . Since there exists a non-trivial solution τ , the local stiffness C ′ of the imaginary material can not be either negative definite (notation ≺ 0) or positive definite (notation ≻ 0) for the whole physical space, or equivalently

∃x | C(x) -(1 -λ)C 0 ⊀ 0, ∃x | C(x) -(1 -λ)C 0 ⊁ 0. (2.51)
Equation (2.51) is identical to the conditions obtained for the basic scheme in [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF]. Hence, similar results can be deduced

λ ≥ min min x 1 - κ(x) κ 0 , min x 1 - µ(x) µ 0 , λ ≤ max max x 1 - κ(x) κ 0 , max x 1 - µ(x) µ 0 . (2.52)
For a two phase material with the same Poisson ratio in both phases C 1 = εC 0 , we can deduce that

|λ| ≤ |1 -ε| . (2.53)
It is clear that for high contrast ratios, for example ε > 2, we can not guarantee the convergence of the series. In this case, the convergence of the series can only be achieved by using another reference material than the one of the matrix.

Analogously, considering the SS series with the eigen equation

-C 0 (S -S 0 )Q * τ = λτ , (2.54)
we can show that τ is the solution of the new boundary value problem with local compliance S ′ and zero stress Σ = 0, namely

-C 0 S ′ -S 0 Q * τ = τ , S ′ = 1 λ S - (1 -λ) λ S 0 . (2.55)
As a result, we can deduce the similar condition for λ

|λ| ≤ 1 -ε -1 , (2.56)
which means that the SS series can diverge at small contrast ratio, say ε < 1/2.

Considering the general series, we shall bound the eigenvalue λ issued from the equation

-(I -L)(C -C 0 )S 0 P * τ -LC 0 (S -S 0 )Q * τ = λτ . (2.57)
Like for the two previous cases, our strategy is to find the imaginary material that admits non-trivial solution τ with zero average strain E = 0 and stress Σ = 0. Our technique is to assume that the relations

- 1 λ (I -L)(C -C 0 )S 0 P * τ = Mτ , - 1 λ LC 0 (S -S 0 )Q * τ = (I -M)τ , (2.58)
hold true at the same time and compute the corresponding isotropic tensor M(x). In order to do that, we have to rearrange (2.58) in the following way

-(C ′ -C 0 )S 0 P * τ = τ , -C 0 (S ′ -S 0 )Q * τ = τ (2.59)
where C ′ and S ′ are the stiffness and compliance tensors

C ′ = 1 λ M -1 (I -L)(C -C 0 ) + C 0 S ′ = 1 λ S 0 (I -M) -1 LC 0 (S -S 0 ) + S 0 . (2.60)
The next step is to choose suitably M so that C ′ and S ′ are compatible, meaning that they are stiffness and compliance of the same imaginary material

C ′ (x)S ′ (x) = I, ∀x. (2.61)
This choice will guarantee that τ corresponds to the same boundary value problem. As a consequence, τ satisfies automatically both constraints (2.59) at the same time. It is clear that if x belong to the matrix, this condition is automatically verified because

C ′ (x) = C 0 ≻ 0, S ′ (x) = S 0 ≻ 0. (2.62)
The compatibility condition can be recast in the form

1 λ M -1 (I -L)(δC)S 0 + I 1 λ (I -M) -1 LC 0 (δS) + I = I. (2.63)
After finding M from (2.63), it is sufficient to bound λ via the positive-definite and negative definite properties of C ′ and S ′ as before. Due to the fact that C ′ and S ′ are already positive definite in the matrix (2.62), those quantites can not be positive definite in the inclusion

1 λ M -1 (I -L)(δC) + C 0 ⊁ 0. (2.64)
Since (2.63) involves isotropic tensors, it is better to use a decomposition based on orthogonal tensors J and K. Assuming that

M = 3β ′ J + 2α ′ K (2.65)
and posing ε m = µ 1 /µ 0 and ε k = κ 1 /κ 0 , the corresponding equation for α ′ , β ′ reads

(1 -2α) (ε m -1) 2λα ′ + 1 1 -2α (ε m -1) ε m λ(1 -2α ′ ) = 1, (1 -3β) (ε k -1) 3λβ ′ + 1 1 -3β (ε k -1) ε k λ(1 -3β ′ ) = 1.
(2.66) Solving (2.66) for α ′ , β ′ yields the solution

2α ′ = ε m λ(1 -2α) -2α(1 -2α)(ε m -1) [(1 -2α)ε m + 2α]λ , 3β ′ = ε k λ(1 -3β) -3β(1 -3β)(ε m -1) [(1 -3β)ε m + 3β]λ . (2.67)
Substituting α ′ , β ′ back to (2.64), we obtain the condition

2α(ε m -1)/ε m -λ (1 -2α)(ε m -1) + λ ≥ 0, or 3β(ε k -1)/ε k -λ (1 -3β)(ε k -1) + λ ≥ 0. (2.68)
In any case, we must have the inequalities

|λ| ≤ max {|(1 -2α)(ε m -1)|, |2α(ε m -1)/ε m |} or |λ| ≤ max {|(1 -3β)(ε k -1)|, |3β(ε k -1)/ε k |} (2.69)
which bound the spectral radius of the general operators. To optimize the acceptable value of |λ|, one must have

2α = ε m /(ε m + 1), 3β = ε k /(ε k + 1), (2.70)
so that we can obtain the inequalities

|λ| ≤ max ε m -1 ε m + 1 , ε k -1 ε k + 1 , (2.71)
which guarantees the convergence of the series for any finite contrast. Having derived the optimal values of α and β, we can write explicitly the expressions for A and B as follows

A = 3 ε k -1 ε k + 1 κ 0 J + 2 ε m -1 ε m + 1 µ 0 K + ε k -1 ε k + 1 J + ε m -1 ε m + 1 K C e , B = ε k -1 ε k + 1 J + ε m -1 ε m + 1 K E 1 -E 2 + E 3 -E 4 - 2(3κ 0 -2µ 0 ) 3κ 0 + 4µ 0 E 6 .
(2.72)

To differentiate from the ON and OR series obtained by different optimization methods, the series issued from (2.72) will be associated to the notation OD (optimization by direct method). Table 2.1 shows numerical examples corresponding to the three methods of optimization. Taking the properties of the constituents as input, one can compute the values α and β in order to minimize the quantities B (ON) or ρ(B) (OR) or ρ(B ′ ) (OD). One can remark that the three methods can yield different results. Those differences can have a impact on the estimates presented in the next section. However, it is interesting to note that at the infinite contrast limits, numerical results show that the coefficients α, β tend to the same value (0, 0) for void inclusion and (1/2, 1/3) for rigid inclusion.

Estimation of the overall elastic properties

To determine the effective stiffness tensor C e , we need to find the average τ over the inclusion domain τ Ω . The latter is connected to C e via the relation

f τ Ω = (C e -C 0 )E.
(2.73)

On the other hand, τ Ω can be obtained by averaging (2.38), for example

τ Ω = ∞ j=0 D j AE, (2.74)
in which the tensors D 0 , D 1 , .. are determined with the formulas

D 0 = (Bχ) 0 Ω = I, D 1 = (Bχ) 1 Ω = f -1 ξ χ(-ξ)B(ξ)χ(ξ), D 2 = (Bχ) 2 Ω = f -1 ξ χ(-ξ)B(ξ) ξ ′ χ(ξ -ξ ′ )B(ξ ′ )χ(ξ ′ ),
...

D j = (Bχ) j Ω = f -1 ξ 1 ,ξ 2 ,...,ξ n χ(-ξ 1 )χ(ξ 1 -ξ 2 )...χ(ξ j-1 -ξ j )χ(ξ 1 ) B(ξ 1 )...B(ξ j ). ( 2
C e -C 0 = f ∞ j=0 D j [(I -L)(δC) -LC 0 (δS)C e ].
(2.76)

Generally, in numerical applications, we need to truncate the series up-to a sufficiently high value n

C e ≃ I + f n j=0 D j LC 0 (δS) -1 I + f n j=0 D j (I -L)(δC)S 0 C 0 . (2.77)
The estimate (2.77), based on truncation, works well if the remainder of the series is negligible. It can be improved by a better treatment of the series [START_REF] Qd To | A numerical-analytical coupling computational method for homogenization of effective thermal conductivity of periodic composites[END_REF][START_REF] Mt Nguyen | Conductivity estimates of sphericalparticle suspensions based on triplet structure factors[END_REF]. Indeed, repeating the recurrence at step n we obtain an equation for τ

τ = χ n-1 j=0 (Bχ) j A + (χB) n τ . (2.78)
Averaging both sides over the inclusion volume and making the approximation

τ ≃ χ τ Ω , (2.79)
we obtain the following equation for τ

Ω τ Ω ≃ n-1 j=0 D j AE + D n τ Ω . ( 2
C e ≃ (I -D n ) + f n-1 j=0 D j LC 0 (δS) -1 (I -D n ) + f n-1 j=0 D j (I -L)(δC)S 0 C 0 . (2.81)
In the case where the effective material is isotropic or at least cubic, we can extract the main shear modulus µ e and the bulk modulus κ e using the expressions

µ e µ 0 ≃ 1 + f (1 -2α)δµ/µ 0 + 2αδµ/µ 1 1-2αn n-1 j=0 2α j -2αf δµ/µ 1 κ e κ 0 ≃ 1 + f (1 -3β)δκ/κ 0 + 3βδκ/κ 1 1-3βn n-1 j=0 3β j -3βf δκ/κ 1 , µ e = C e 1212 , 3κ e = C e 1111 + 2C e 1122 , α j = D j 1212 , 3β j = D j 1111 + 2D j 1122 , (2.82)
From (2.82), it is interesting to remark that all the microstructure information is contained in the parameters 1-2αn n-1 j=0 2α j and 1-3βn

n-1 j=0 3β j .

Distributions of non overlapping spheres

We assume now that the unit cell V of dimensions l × l × l contains N identical non overlapping spheres of radius R. The shape functions χ(ξ) become

χ(ξ) = V s V F (ξ)ρ(ξ), ρ(ξ) = N i=1 e -iξ.x i , (2.83)
where x i is the center location of the inclusion numbered i in the cell. The form factor F (ξ) and inclusion volumes V s admit the following form

V s = 4π 3 R 3 , F (ξ) = 3 sin η -η cos η η 3 , η = |η|, η i = Rξ i = 2πn i R/l, i = 1, 2, 3 (2.84)
For cubic arrangements of spheres (see Fig. 2.3), the formulas of χ(ξ) are known explicitly [START_REF] Nemat-Nasser | On composites with periodic structure[END_REF][START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF].

-Simple Cubic (SC)

χ(ξ) = 3f sin η -η cos η η 3 (2.85) -Body Centered Cubic (BC) χ(ξ) = 3f 2 sin η -η cos η η 3 [1 + (-1) n 1 +n 2 +n 3 ] (2.86) -Face Centered Cubic (FCC) χ(ξ) = 3f 4 sin η -η cos η η 3 [(-1) n 1 + (-1) n 2 + (-1) n 3 + (-1) n 1 +n 2 +n 3 ]
(2.87)

Using these expressions, we can compute numerically the intermediate tensors D n at any order n (see Appendix B). The first order tensor D 1 can be expressed more explicitly and related to Hashin-Shtrikman estimates, as shown in the next subsection. It is possible to derive analytically expressions for the first two tensors D 0 = I (or 2α 0 = 1 and 3β 0 = 1) and D 1 [START_REF] Qd To | Explicit effective elasticity tensors of two-phase periodic composites with spherical or ellipsoidal inclusions[END_REF] which is a lattice sum. In the case where D 1 is a cubic tensor, we can compute explicitly the coefficients β 1 and α 1 associated to D 1

2α 1 = -(1 -2α) δµ µ 0 2 3 (1 -f ) -S 3 4(3κ 0 + µ 0 ) 3κ 0 + 4µ 0 - -2α δµ µ 1 1 3 (1 -f ) + S 3 4(3κ 0 + µ 0 ) 3κ 0 + 4µ 0 , 3β 1 = -(1 -f ) (1 -3β) δκ κ 0 3κ 0 3κ 0 + 4µ 0 -3β δκ κ 1 4µ 0 3κ 0 + 4µ 0 (2.88)
via the lattice sums S 1 , S 2 and S 3 defined as

S 1 = f -1 ξ ξ2 i χ(ξ)χ(-ξ), i = 1, 2, 3 S 2 = f -1 ξ ξ4 i χ(ξ)χ(-ξ), i = 1, 2, 3 S 3 = f -1 ξ ξ2 i ξ2 j χ(ξ)χ(-ξ), i, j = 1, 2, 3, i = j (2.89)
Expression (2.88) has been obtained by making use of the additional identities

S 1 = (1 -f )/3, S 2 + 2S 3 = S 1 (2.
κ 0 δκ + (1 -f ) 3κ 0 3κ 0 +4µ 0 (2.91)
Analogously, all the estimates of µ e are reduced at order 1 to

µ e µ 0 = 1 + f µ 0 δµ + 2 3 (1 -f ) -4S 3 3κ 0 +µ 0 3κ 0 +4µ 0 (2.92)
where the lattice sum S 3 can be evaluated by semi-analytical expressions [START_REF] Qd To | Explicit effective elasticity tensors of two-phase periodic composites with spherical or ellipsoidal inclusions[END_REF]. We will show later that in the case of an isotropic random distribution, this first order estimate of µ e is identical to HS estimate.

For higher orders n ≥ 2, numerical summation can be used to compute D n which in turn provides better estimates of κ e and µ e .

Random distributions of non overlapping spheres and relation with structure factors

Regarding random distribution (RD) in Fig. 2.3, we adopt the ergodicity hypothesis implying that at the infinite volume limit (both N and V → ∞), the tensors D j are identical to their ensemble average, both being denoted by ... . As a result, we obtain now statistical relations

D 1 = V s V ξ F (ξ)F (-ξ)B(ξ)S (2) (ξ) D 2 = V 2 s V 2 ξ,ξ ′ F (ξ)F (ξ ′ )F (ξ -ξ ′ )B(ξ)B(ξ ′ )S (3) (-ξ, ξ ′ ), etc. (2.93)
with S (2) , S (3) being the structure factor and the triplet structure factor given by

S (2) (ξ) = 1 N ρ(ξ)ρ(-ξ) , S (3) (ξ, ξ ′ ) = 1 N ρ(ξ)ρ(ξ ′ )ρ(-ξ -ξ ′ ) , etc. (2.94)
We note that these results are written using the structure factors which are important quantities in condensed matter physics and crystallography since they are directly related to statistical information, i.e the relative arrangement of the inclusions. For example, S (2) (ξ) are connected to radial distribution g (2) (r) via the expression

S (2) (ξ) = 1 + ρ V g (2) (r)e -iξ.r dr (2.95)
where ρ is the average inclusion density. For a system of non overlapping spheres in equilibrium, analytical expressions for both g (2) and S (2) exist from the solution of Ornstein-Zernike equation [START_REF] Ls Ornstein | Accidental deviations of density and opalescence at the critical point of a single substance[END_REF] and Percus-Yevick closure approximation [START_REF] Percus | Analysis of classical statistical mechanics by means of collective coordinates[END_REF][START_REF] Ms Wertheim | Exact solution of the Percus-Yevick integral equation for hard spheres[END_REF]. More importantly, S (2) (ξ) can be determined experimentally via scattering techniques. Higher order structure factors like S (3) , S (4) ,... are more difficult to obtain by experiment [START_REF] Wochner | X-ray cross correlation analysis uncovers hidden local symmetries in disordered matter[END_REF] but they can be obtained by approximations [START_REF] Hansen | Theory of simple liquids[END_REF][START_REF] Mt Nguyen | Conductivity estimates of sphericalparticle suspensions based on triplet structure factors[END_REF][START_REF] Jl Barrat | Factorization of the triplet direct correlation function in dense fluids[END_REF][START_REF] Kirkwood | Statistical mechanics of fluid mixtures[END_REF] or atomistic based computer simulations.

In this work, we are concerned with estimates based on the first two structure factors S (2) (ξ) and S (3) (ξ) which already perform well. The sums (2.93) can be rewritten as integrals in the Fourier space as follows

D 1 = 1 6π 2 F (ξ)F (-ξ)B(ξ)S (2) (ξ)dη D 2 = 1 36π 4 F (ξ)F (ξ ′ )F (ξ -ξ ′ )B(ξ)B(ξ ′ )S (3) (-ξ, ξ ′ )dηdη ′ . (2.96)
For isotropic distributions of non overlapping spheres, we have the following property

1 6π 2 F 2 (ξ)S (2) (ξ)dη = (1 -f ) (2.97)
regardless the local pair distribution S (2) . As in the discrete case, the expression of the equivalent main shear modulus is given by expression 2.92, but using the continuous equivalent of S 3 given by :

S 3 = 1 6π 2 F (ξ)F (-ξ) ξ2 i ξ2 j S (2) (ξ)dη (2.98)
Using the general relation 2.97, we can finally evaluate S 3

S 3 = 1 -f 15 (2.99)
and finally (93), giving the first order evaluation of the main shear modulus µ e coincides with the HS bound for µ e [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF]. As seen thereafter, the evaluation at the second order departs from HS bound .

Numerical applications and analysis

Let us consider the application of our theory to a cubic array of spheres in comparison with FFT and some literature results. From Figs. 2.4, 2.5 and Tabs. 2.2, 2.3, we find that the second order estimates have improved significantly the HS bound which coincides with the first order estimates. The degree of improvement depends on the properties considered, estimation scheme, microstructure and the elastic properties of constituents. For BCC array with contrast ratio as high as 10 (see Fig 2 .4, 2.5)), the second order estimates of the bulk modulus κ e /κ 0 and the shear modulus µ e /µ 0 issued from the three schemes are close to the FFT results at convergence. The agreement is good upto a very high volume fraction near the percolation limit. Detailed results on FCC array also have the same trend as those for BCC cases. Tables 2.2 and 2.3 show that the two series OR and OD yield very good results while the series ON works less well. It is interesting to note that at volume fraction as large as 0.5 and the rigidity contrast ratio as high as 100, our three estimates perform well.

Next, we study microstructures constituted of randomly isotropic distribution of spheres. Two extreme cases of rigid spheres and voids will be considered. Fifty sample composed of 500 non overlapping spheres are prepared by standard Event Driven Molecular Dynamics [START_REF] Dc Rapaport | The Art of Molecular Dynamics Simulation[END_REF]. To compute the effective properties of the material, we shall limit to OD based estimates and final results are obtained by averaging over the 50 samples. It may be noticed that, the iterative scheme is not theoretically convergent for fields in void or rigid inclusions but the effective properties exist for the considered microstructures. In this case, the expressions for µ e and κ e at first and second order can be used for infinite contrast. Simulations on the systems show that the OD-2 estimate is close to the estimate of [START_REF] Torquato | Effective stiffness tensor of composite media : II. Applications to isotropic dispersions[END_REF][START_REF] Torquato | Effective stiffness tensor of composite media. I. Exact series expansions[END_REF]. For spherical voids (see Fig. 2.9), our second order estimate again shows a significant improvement with respect to the first order (HS bound). The estimate is also close to [START_REF] Torquato | Effective stiffness tensor of composite media : II. Applications to isotropic dispersions[END_REF][START_REF] Torquato | Effective stiffness tensor of composite media. I. Exact series expansions[END_REF] using three point parameters. Those results again confirm the robustness of our estimation scheme at high rigidity contrast and high volume fraction. We note that the good performance comes from the benefit of the fast convergence series and the high order correlation information.

Concluding remarks

In this paper, we have presented a new estimate of the overall stiffness tensor of elastic composites. Starting from a class of Lippmann-Schwinger integral equations for eigens- tress (or eigenstrain), the optimization procedure is then carried out to find the best Neumann series, i.e those with the fastest convergence rate. To this end, we have introduced tools to bound the spectral radius and norm of fourth order operators in Fourier space and methods to obtain the optimal series. The series are then used to derive estimates at different order n. We have also shown that n-order statistical information on the microstructure, in this case corresponding to the structure factors, also appear in the estimates. Numerical applications of the procedure on some test cases show that our estimates perform very well in comparison with FFT results and those from the literature.

Not only concerning the theoretical estimates, the contribution of the present paper is closely related to the FFT resolution method. Results issued from the paper correspond to new computation methods based on the eigenstress whose theory is rigorously founded and convergence rate is controllable. Indeed, by not taking the reference material as one of the constituents while keeping all the remaining ingredients, we obtain a very general and promising class of series which can be developed into new FFT schemes. The next step is to optimize and examine the convergence performance in comparison with the existing FFT schemes. Those perspectives will be investigated in a future work.

ment (NAFOSTED) under grant number 107.02-2016.05.

Appendix

Appendix A : Lippmann-Schwinger equations in elasticity and associated FFT resolution method

The boundary value problem can be solved using either one of the two Lippmann-Schwinger integral equations

ǫ = E -Γ 0 δ ′ Cǫ, σ = Σ -∆ 0 δ ′ Sσ (2.100)
It is clear that these two integral equations are equivalent to (2.29) and (2.30) via the relations τ = δ ′ Cǫ and e = δ ′ Sσ. Taking the first equation (2.100) as an example, ǫ can be computed by repeating the iterative scheme

ǫ 0 = E, ǫ n+1 = E -Γ 0 δ ′ Cǫ n (2.101)
until convergence. Since we know the explicit expression of Γ 0 in the Fourier space and the material distribution δ ′ C in the physical space, it is preferable to use the FFT algorithm to convert results back and forth between the two spaces. In most of our applications, the resolution 256 × 256 × 256 is adopted for both spaces. Finally, improving the resolution and increasing the number of iterations up to convergence, we obtain the exact solution of the problem within any given accuracy. We note that the adopted FFT method is related to ES series and different from series used as the background of the estimation presented in the paper. Using this classic scheme, the bulk stiffness and the shear stiffness of the reference material are chosen as averages of the related moduli of inclusion and matrix. Details of the scheme can be found in previous works [START_REF] Michel | Effective properties of composite materials with periodic microstructure : a computational approach[END_REF][START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF].

Appendix B : Method of computing D n and explicit expressions for D 1

To compute numerically tensor D n (j ≥ 1), we need to base on (2.75) which can be recast as

D j = f -1 ξ χ(ξ)F j (ξ)χ(-ξ), F 1 (ξ) = B(ξ), F j (ξ) = B(ξ)(χ * F j-1 )(ξ), j ≥ 2, ( 2.102)
where * stands for convolution in Fourier space. Given the explicit expressions of B and χ in Fourier space, the convolution can be done efficiently using FFT techniques.

To compute the cubic tensor D 1 , we need to calculate first the following elementary H i to Walpole base elements E i

H i = f -1 ξ E i (ξ)χ(ξ)χ(-ξ), i = 1, 2, .., 6 (2.103)
The explicit Mandel matrix representation of those tensors in terms of S 1 , S 2 and S 3 are given as follows

[H 1 ] =         1 2 (S 1 + S 2 ) 1 2 (S 1 + S 3 ) 1 2 (S 1 + S 3 ) 0 0 0 1 2 (S 1 + S 3 ) 1 2 (S 1 + S 2 ) 1 2 (S 1 + S 3 ) 0 0 0 1 2 (S 1 + S 3 ) 1 2 (S 1 + S 3 ) 1 2 (S 1 + S 2 ) 0 0 0 0 0 0 S 3 0 0 0 0 0 0 S 3 0 0 0 0 0 0 S 3         , (2.104) [H 2 ] =         S 2 S 3 S 3 0 0 0 S 3 S 2 S 3 0 0 0 S 3 S 3 S 2 0 0 0 0 0 0 2S 3 0 0 0 0 0 0 2S 3 0 0 0 0 0 0 2S 3         , (2.105) [H] 3 =         1 2 (S 1 + S 2 ) 1 2 (S 3 -S 1 ) 1 2 (S 3 -S 1 ) 0 0 0 1 2 (S 3 -S 1 ) 1 2 (S 1 + S 2 ) 1 2 (S 3 -S 1 ) 0 0 0 1 2 (S 3 -S 1 ) 1 2 (S 3 -S 1 ) 1 2 (S 1 + S 2 ) 0 0 0 0 0 0 S 1 + S 3 0 0 0 0 0 0 S 1 + S 3 0 0 0 0 0 0 S 1 + S 3         , (2.106) [H] 4 =         2(S 1 -S 2 ) -2S 3 -2S 3 0 0 0 -2S 3 2(S 1 -S 2 ) -2S 3 0 0 0 -2S 3 -2S 3 2(S 1 -S 2 ) 0 0 0 0 0 0 2(S 1 -2S 3 ) 0 0 0 0 0 0 2(S 1 -2S 3 ) 0 0 0 0 0 0 2(S 1 -2S 3 )         , (2.107) [H] 5 = [H] 6 =         S 1 -S 2 S 1 -S 3 S 1 -S 3 0 0 0 S 1 -S 3 S 1 -S 2 S 1 -S 3 0 0 0 S 1 -S 3 S 1 -S 3 S 1 -S 2 0 0 0 0 0 0 -2S 3 0 0 0 0 0 0 -2S 3 0 0 0 0 0 0 -2S 3         , (2.108)
From the above results, given any tensor B(ξ) with representation B(ξ) = 6 i=1 b i E i (ξ), we can derive D 1 = 6 i=1 b i H i and the parameters α 1 and β 1 in (2.88).

Introduction

The Brinkman equation [START_REF] Hc Brinkman | On the permeability of media consisting of closely packed porous particles[END_REF][START_REF] Hc Brinkman | A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles[END_REF] has been adopted to simulate complex fluid flows in porous solids with various challenging applications. The Laplacian has been introduced in the Darcy equation in order to describe non-slip boundary conditions along the solid wall which confine the porous solid or the conditions at the frontier with a free fluid (note that alternative methods use the Beaver-Joseph-Saffman interfacial model [START_REF] Beavers | Boundary conditions at a naturally permeable wall[END_REF][START_REF] Pg Saffman | On the boundary condition at the surface of a porous medium[END_REF]). The domain of application of the Brinkman equation has concerned fluid flow in porous media featuring open channels or fractures [START_REF] Lesinigo | A multiscale darcy-brinkman model for fluid flow in fractured porous media[END_REF], the lubrication theory [START_REF] Silva | Stokes-brinkman-darcy solutions of bimodal porous flow across periodic array of permeable cylindrical inclusions : cell model, lubrication theory and lbm/fem numerical simulations[END_REF], the mass transfer in biological system or in modern technologies (see for instance the books of Truskey et al. 2004 [START_REF] Ga Truskey | Transport phenomena in biological systems[END_REF] and Bejan et al. 2013 [START_REF] Bejan | Porous and complex flow structures in modern technologies[END_REF]).

The Brinkman equation introduces two coefficients, the permeability (or resistivity) and the effective dynamic viscosity µ e that is different of the dynamic viscosity µ which enters in the Stokes equations. The determination of the effective dynamic viscosity has been the subject of many works, experimental ones [START_REF] Givler | A determination of the effective viscosity for the brinkman-forchheimer flow model[END_REF] or numerical approaches [START_REF] Martys | Computer simulation study of the effective viscosity in brinkmans equation[END_REF][START_REF] Vm Starov | Effective viscosity and permeability of porous media[END_REF]. More recently Auriault et al. 2005 [3] derived the coefficients of the Brinkman law in the framework of the asymptotic homogenization. Starting from the Stokes equation at the local scale and expanding the solution with the scale parameter, Auriault and Sanchez-Palencia 1977 [START_REF] Jl Auriault | Study of the macroscopic behavior of a saturated porous elastic medium[END_REF], Sanchez Palencia 1980 [START_REF] Sanchez-Palencia | Non-homogeneous media and vibration theory[END_REF], Levy 1983 [34], etc. have derived the famous Darcy's Law. By exploiting the higher order terms of the series, Auriault et al. 2005 [3] found that the correction to the Darcy law introduce the dependence with the double, triple etc gradient of pressure. Moreover, substituting the higher order derivatives of the pressure gradient by those of the macroscopic velocity, and considering the case of a macroscopic isotropic porous material they retrieve the Brinkman equation.

It is interesting to note that the results obtained with the Stokes equations has some similarities with those obtained for elastic composites ten years before by Boutin 1996 [START_REF] Boutin | Microstructural effects in elastic composites[END_REF] with the difference that, in the elasticity problem, the correctors depends on the double, tripe, etc gradient of macroscopic displacement. The consideration of these higher order solutions led to regularized elasticity models (see [START_REF] Vp Smyshlyaev | On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media[END_REF][START_REF] Th Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF])

In this paper we provide a FFT based numerical method for the resolution of the hierarchy of unit cell problem at any order. There is various reason to develop the FFT method instead of conventional numerical methods based on finite elements or volume elements. First, the method based on FFT do not require any special treatment due to the incompressibility of the flow while the FEM need the introduction of a bubble function to make the solution stable. Next, the FFT method is less memory saving and open interesting possibilities for the application to complex 3D microstructure (see for instance Ly et al. 2016 [START_REF] Hb Ly | Computation of permeability with fast fourier transform from 3-d digital images of porous microstructures[END_REF]).

The paper is organized as follows. In section 3.2 we recall the main results obtained by Auriault 2005 that we reformulate in a form more adapted for the implementation with the FFT method. In section 3.3, we derive the "Brinkman type" macroscopic filtration law written for an anisotropic centrosymetric porous medium. The method of resolution of the unit cell problems is presented in section 3.4. In order to validate the accuracy of the method we compare the numerical solution with analytic ones in the case of the flow along parallel pores in section 3.5. In the last section, we derive the anisotropic Brinkman coefficients in the case of the 2D transverse flow in a regular array of cylinders.

Hierarchy of local problems

Consider a periodic porous medium saturated by a homogeneous Newtonian viscous fluid with the dynamic viscosity µ. By V , we denote the total volume of the unit cell, by V f and V s the volume occupied by the fluid and the solid respectively, where V f and V s are both interconnected domains. The frontier between the fluid and the solid is denoted Γ. Due to the fluid motion, the mass, momentum and energy can be transported along with the stream currents, which results in the following set of local equations :

   µ∆v -∇p = 0 in V f div v = 0 in V f v = 0 on Γ (3.1)
in which v and p represent the local velocity and pressure fields.

Classically in fluid flow in porous materials, we assume that the Reynolds number is low, the inertia term is then neglected in the first equation in (3.1). Following Auriault and Sanchez-Palencia 1977, Sanchez Palencia 1980, Levy 1983 and more recently, Auriault 2005 [3], the solution of the Stokes equations are determined by using the methods based on asymptotic series. Each term of the expansion series for the velocity field are determined by solving a hierarchy of local problems which are fully computed on a unit cell of the porous solid. Here, the main results of Auriault 2005 are recalled but in a different form which introduces the strain rate tensors and the stress tensors at different order. The main reason is inherent with the use of the FFT method for solving the local problems. In The FFT method, a continuation by continuity of the velocity and the stress field within the solid phase is needed that requires to introduce explicitly the stress fields at each order. The strain rate and stress tensors are given by :

d = ∇ s v, σ = 2µd -pI (3.2)
where ∇ s denote the symmetric part of the gradient and I is the two order identity tensor.

The strain rate tensor is traceless due to incompressibility while the stress comply with the local equilibrium : div(σ) = 0 (3.

3)

The cornerstone of the asymptotic homogenization is to define a scale factor :

ǫ = h L (3.4)
in which h and L are the characteristic lengths of the micro-and macrostructures respectively. In order to formulate the equations of the macroscopic law, we need to introduce non-dimensional quantities and to evaluate the magnitude of each physical quantities at the local scale with respect to the scale factor ǫ. To this end, the following change of variable is considered :

v * = 1 v c v, p * = p p c , y = 1 h X, x = 1 L X, µ = µ * µ c (3.5)
where v c and p c are the characteristic velocity and pressure while x and y are the low and rapid coordinates. The solution is expanded along power series in ǫ :

v(x, y) = v c +∞ i=0 ǫ i v i (x, y), p(x, y) = p c +∞ i=0 ǫ i p i (x, y) (3.6)
The differential operators are split with respect to coordinates x and y respectively :

∇ = ∇ x + 1 ǫ ∇ y .
The strain rate tensor reads :

d(x, y) = v c h +∞ i=0 ǫ i d i (x, y) (3.7)
in which the strain rate tensors at each order are given by :

d 0 = d y (v 0 ), d 1 = d y (v 1 ) + d x (v 0 ), d 2 = d y (v 2 ) + d x (v 1 ), ... (3.8)
and where d y and d x represents the symmetric part of the gradient with respect to y and x respectively. Assuming now that [3] :

µ c v c hp c = ǫ (3.9)
We deduce that :

p c = µ c v c hǫ = µ c v c Lǫ 2 (3.10)
Using this result in the computation of the stress, we can read :

σ(x, y) = µ c v c h +∞ i=-1 ǫ i σ i (x, y) (3.11)
in which :

σ -1 = -p 0 I, σ 0 = 2µ * d 0 -p 1 I, σ 1 = 2µ * d 1 -p 2 I, ... (3.12)
Finally, the local equilibrium for the stress reads :

1 ǫ div y (σ) + div x (σ) = 0 (3.13)
All the term with the same power in ǫ are afterward regrouped to get a hierarchy of local problems. The method has been well documented in the literature and the reader could refer to Auriault 2005 and Auriault, et al. 2010 [2]. The problem in ǫ -2 is ∇ y p 0 = 0 which suggests that p 0 is only function of the coordinates x and ∇ x p 0 is the macroscopic pressure gradient and is denoted by thereafter J . The problem in ǫ -1 (the Darcy problem) reads :

       div y (σ 0 ) = J in V f σ 0 = -p 1 I + 2µ * d y (v 0 ) in V f tr[d y (v 0 )] = div y (v 0 ) = 0 in V f v 0 = 0 on Γ (3.14)
The solution linearly depends on the applied macroscopic pressure gradient J . Considering the average of the velocity field v 0 , one recover the Darcy equation :

V 0 = c f < v 0 > V f = - 1 µ * K 0 .J (3.15)
where c f is the volume fraction of the fluid and <> V f represent the average over the volume V f . Eq. (3.15) gives the first approximation for the macroscopic flux. Investigating the higher order terms of the expansion series, we observe that the equilibrium equation for the stress fields

σ i with i ≥ 1 is div y (σ i ) + div x (σ i-1 ) = 0 in V f (3.16)
The volume balance equation for the velocity field v i with i ≥ 1 is :

div y (v i ) + div x (v i-1 ) = 0 in V f (3.17)
As shown in [3], the integration of the above balance equation over V f leads to :

div x (< v i-1 > V f ) = 0 (3.18)
where the term < div y (v i ) > V f vanishes since the velocity v i is null on Γ. The above equation represents the volume balance equation at the macroscopic level. Particularly, for i = 1, this equation becomes div x (V 0 ) = 0 where V 0 has been introduced in (3.15). Doing the same with the higher order term in the expansion series of v, we denote by

V i = c f < v i > V f for i = 1, 2
, .... Since the bracket <> V represents the integration over the volume of the cell (with respect to coordinate y), the macroscopic velocity V i is only function of the coordinate x. These vectors represent the correction to the Darcy law.

Owing to relation (3.18), all these vectors satisfy to div x (V i ) = 0. Following [3], equation (3.18) is removed from the mass balance equation (3.16), this leading to :

div y (v i ) + div x v i-1 - 1 c f V i-1 = 0 in V f (3.19)
The relations (3.12) must be also modified to be compatible with the above relation. Indeed, in (3.12), the strain fields d i for i = 0, 1, 2, ... are traceless and -p i for i = 1, 2, ... are the trace of the stress tensor σ 0 , σ 1 , etc. However, with the modification (3.19) made in the mass balance equation, the strain fields d i are not traceless anymore. Consequently, the following definition must be now used for the local strain fields :

d i = d y (v i ) + d x (v i-1 ) - 1 3 tr(D i-1 ) c f I, i ≥ 1 (3.20)
with the following notation :

D i = d x (V i ), i ≥ 1 (3.21)
Note that this has no incidence on the equilibrium equation since the extra term in (3.20) is proportional to I and the solution of each cell problem for the pressure fields p i are defined up to an added constant pressure.

Considering the last results, all higher order cell problems can be recast in the following form :

               div y (σ i ) + div x (σ i-1 ) = 0 in V f σ i = -p i+1 I + 2µ * d y (v i ) + d x (v i-1 ) - 1 3 tr(D i-1 ) c f I in V f tr d y (v i ) + d x (v i-1 ) - 1 3 tr(D i-1 ) c f I = 0 in V f v i = 0 on ∂V f (3.22)
with i = 1, 2, 3.... These higher order cell problems introduce the dependence of the local solution with the higher order derivatives of the macroscopic pressure, i.e. with the double, triple,... gradient of the macroscopic pressure. For instance, in the second order problem, the enforcing terms depends on the derivative of the first order solution (v 0 , σ 0 ) with respect to x. The solution at first order linearly depends on the pressure gradient J i , then the solution at second order (v 1 , σ 1 ) depends on J ij which represents the components of the double gradient of the macroscopic pressure. In the third order problem, the enforcing terms are given by the derivative of (v 1 , σ 1 ) with respect to x, the solution (v 2 , σ 2 ) then linearly depends on the triple gradient of the macroscopic pressure whose components are J ijk .

Macroscopic filtration law 3.3.1 Derivation of the macroscopic velocity field

The macroscopic filtration law is derived by computing the flux over each surface of the unit cell :

V i = 1 Σ i Σ f i v i dy (3.23)
where Σ i represent the surface of the cubical unit cell of normal unit vector e i , Σ f i is the surface of Σ i crossed by the fluid. As shown in [3], the above formula can be transformed, using the divergence theorem, into :

V =< v > V + 1 V V y div y (v)dy (3.24)
The macroscopic velocity is not equal to the average of the local velocity. However, when only the first term of the series v 0 is kept in the expansion series of the velocity field, we have V =< v 0 > V = V 0 since div y (v 0 ) = 0. When the higher order terms of the series are considered, div y (v i ) = 0 and the equality V =< v > V is no valid anymore. The solution at first order linearly depends on the applied pressure gradient J i . The solution at second, third, etc order linearly depends on double, triple, etc gradient of pressure. As a consequence, the macroscopic velocity computed from eq. (3.24) reads :

V i = - 1 µ K 0 ij J j + ǫK 1 ijk J jk + ǫ 2 K 2 ijkl J jkl + ... (3.25)
which involves higher order permeabilities K 1 ijk , K 2 ijkl ,... in addition to the classic Darcy permeability K 0 ij . The double, triple, etc, gradient of macroscopic pressure J ij , J ijk , etc, are invariants by any permutation of theirs indices. Consequently, K 1 ijk is invariant by the permutation of j and k, K 2 ijkl is invariant by any permutation of the indices j, k and l,...

Implementation of unit cell problems with the FFT

In this section, we propose a unified iterative scheme to compute the solution of cell problems at any order. The method use the Fourier series and the FFT algorithm to compute the Fourier transform and its inverse. The method has been already developed in the case of fluid flow in porous solids (see for instance [START_REF] Monchiet | A fft-based method to compute the permeability induced by a stokes slip flow through a porous medium[END_REF][START_REF] Tk Nguyen | A fourier based numerical method for computing the dynamic permeability of periodic porous media[END_REF], Ly et al. 2016 [45, 55, 35]) but has not been considered for the derivation of the higher order terms of the series. The computation of higher order terms requires particular development presented in the next of the section.

Continuation of local fields in the solid phase

A first step of the application of the FFT method in the case of fluid flow in porous solid consists to make a continuation of the local fields within the solid phase. This necessary when the Fourier series are used to discretized the cell problem since they must be considered in a rectangular (or cubic for 3D problems) domain. For example, in the case of fluid flow around a rigid cylinder, the fluid phase is not a rectangular region. Since all the velocities v i for i = 1, 2, 3, ... are null on the boundary Γ, the continuation by continuity is made with a null velocity field within the solid region. A continuation of the stress fields σ i for i = 0, 1, 2, ... is considered with the following hypotheses :

• the traction σ i .n is antiperiodic on the boundary of the unit cell,

• the traction σ i .n is continuous at the solid/fluid interface,

• the stress σ i comply with the equilibrium equation :

div y (σ i ) + f i (y) = 0 in V s (3.34)
for i = 0, 1, 2, 3, ... and where the value of the source terms f i is determined in the next of the section. The first two hypotheses are automatically verified when the Fourier series are used. The stress fields σ i are periodic on two opposite side of the cell, so, the traction is antiperiodic on the opposite sides. By taking the average of the divergence of the stress σ i over the volume of the cell V and making use of the divergence theorem we deduce that :

< div y (σ i ) > V = 1 V V σ i .ndx = 0 (3.35)
The divergence of σ i being equal to f i (y), we also deduce that :

< f i (y) > V = 0 (3.36)
Consider the first order cell problem, the source term f 0 (y) is equal to J , the macroscopic pressure gradient, in the fluid region. In the solid phase, the value of f 0 (y) could be chosen arbitrarily as soon as the condition (3.36) is satisfied. The simplest choice consists to consider a constant value for f 0 (y) in the solid phase. Denoting by f 0 s this constant, the verification of (3.36) leads to :

f 0 s = - c f c s J (3.37)
where c f and c s denotes the volume fraction of the fluid and solid phase respectively. The source term f 0 s physically represents the drag force due to the fluid motion. This result was obtained in [START_REF] Monchiet | A fft-based method to compute the permeability induced by a stokes slip flow through a porous medium[END_REF] and the first order cell problem could be rewritten in the form :

   div y (σ 0 ) = f 0 (y) in V d y (v 0 ) = S(y) : σ 0 in V v 0 = 0 in V s (3.38)
where we have introduced :

S = 1 2µ(y) K, µ(y) = µ * in V f +∞ in V s (3.39)
and

f 0 (y) =    J in V f - c f c s J in V s (3.40)
Since µ(y) = +∞ in the solid phase, the strain rate is null in V s and the velocity is given by a rigid motion. Obviously, a large but finite value of the dynamic viscosity could be used in the solid phase. This dynamics viscosity could then be interpreted as a penalty coefficient. However, and as shown later, the method of resolution based on the FFT use the inverse of the dynamic viscosity that is taken equal to zero in the solid region. In this problem, the velocity field being defined up to an added rigid motion, the condition v 0 = 0 in V s must be added. Consider now the higher order cell problems that is the central subject of this work. In eq. (3.22), the solution depends on the local fields v i-1 and σ i-1 computed at lower order. For example when i = 1 the velocity v 0 and the stress σ 0 has been computed from (3.14) or equivalently from (3.38) with the difference that the velocity fields v 0 and the stress field σ 0 are defined anywhere in the cell while the solution from (3.14) is defined only in the fluid region. By making a continuation of the velocity v 1 and the stress σ 1 in V s we obtained the new following equations for the first order cell problem :

   div y (σ i ) = f i (y) in V d y (v i ) = S(y) : σ i -g i in V v i = 0 in V s (3.41)
where :

g i = d x (v i-1 ) - 1 3 tr(D i ) c f I f (y) I (3.42)
and where I f (y) is the characteristic function of the fluid phase which takes the value 1 if y ∈ V f and 0 if not. With the introduction of I f (y), tensor g i is null in V s and d y (v i ) is also null in the V s that is needed to recover v = 0 in V s .

The following expression is used for the source term f i (y) :

f i (y) = -div x (σ i-1 ) in V f f i s in V s (3.43)
where, again, a constant expression is used for f i (y) in the solid phase. Introducing this expression in (3.36) leads to :

f i s = c f c s div x (< σ i-1 > V f ) (3.44)
The case i = 1 and i ≥ 2 must be distinguished. If i = 1, we have :

f 1 s = c f c s div x (< σ 0 > V f ) (3.45)
in which

< σ 0 > V f = 2µ * < d y (v 0 ) > V f -< p 1 > V f I (3.46)
The term < d y (v 0 ) > V f is null since the velocity v 0 is null at the boundary Γ and v 0 .n is antiperiodic on the boundary of the cell. In each local problem, the pressure fields is defined up to an added constant pressure, then, we can put < p 1 > V f = 0 without any modification of the solution in the fluid region. Consequently, f 1 s is null in the solid region and :

f 1 (y) = -div x (σ 0 ) in V f 0 in V s (3.47)
In consequence, there is no drag force for the problem associated with the first corrector. Consider now the case i ≥ 2. The stress σ i-1 is given by (3.22) and we have :

< σ i-1 > V f = 2µ * < d y (v i-1 ) > V f +d x (< v i-2 > V f ) - 1 3 tr(D i-2 )I (3.48)
-< p i > V f I (3.49)
For the same reasons, we put < p i > V f = 0 and < d y (v i-1 ) > V f is null. Since the velocity field is null within the solid phase, we have :

< v i > V f = 1 c f V i (3.50)
It follows that :

< σ i-1 > V f = 2µ * c f D i-2 - 1 3 tr(D i-2 )I = 2µ * c f D i-2 , i ≥ 2 (3.51)
where

D i-2 is the deviatoric part of D i-2 .
The expression of f i (y) is :

f i (y) =    -div x (σ i-1 ) in V f 2µ * c s div x D i-2 in V s , i ≥ 2 (3.52)
The unit cell problems are hierarchized since the the source term at a given order depends on the solution at lower order. These unit cell problems are given by (3.38) and (3.41) in which the source term f i is given by (3.40) for i = 0, by (3.47) for i = 1 and (3.52) for i ≥ 2.

The integral equation and the iterative scheme

The cornerstone of the FFT numerical method is the reformulation of the system of partial differential equations in an integral equation which uses the Green operator and the convolution product. Such formulation of the problem is possible for the Stokes problem only after making the continuation by continuity of the local fields within the porous medium. This has been the subject of the latter section. It has been shown that the cell problems could be put into the following general form :

   div y (σ i ) = f i (y) in V d y (v i ) = S(y) : σ i -g i (y) in V v i = 0 in V s (3.53)
where the value of f i (y) has been given by (3.40), (3.47) and (3.52), and the value of g i (y) is given by (3.42) for i ≥ 1 and is null for i = 0. For simplicity, the subscript "y" in the differential operators is removed. Classically, the integral equation is determined by introducing a homogeneous medium of dynamic viscosity µ 0 and the polarization tensor or equivalently the eigenstrain tensor given by :

d * = (S(y) -S 0 ) : σ i -g i (y) (3.54)
and finally by solving the problem in Fourier space. All the details of this resolution are not provided here and could be found in [START_REF] Monchiet | A fft-based method to compute the permeability induced by a stokes slip flow through a porous medium[END_REF] for instance. The solution for the stress field is :

σ i =< σ i > V +Ω * f i -∆ 0 * d * (3.55)
where " * " represents the convolution product, ∆ 0 is the Green tensor for the stress defined by :

∆ 0 ijkl = 2µ 0 Q ij Q kl + 1 2 (Q ik Q jl + Q il Q jk ) (3.56)
and the components of Ω are :

Ω ijk = i ξ 2 [δ ij ξ k + Q ik ξ j + Q jk ξ i ] (3.57)
with :

Q = I - ξ ⊗ ξ ξ 2 (3.58)
Replacing the eigenstrain d * by its expression given by (3.54), we obtain an integral equation for the stress field :

σ i =< σ i > V +Ω * f i -∆ 0 * (S -S 0 ) : σ i -g i (3.59)
This equation is very similar to that obtain in [START_REF] Monchiet | A fft-based method to compute the permeability induced by a stokes slip flow through a porous medium[END_REF] but contains two additional terms. The first is "g i " which only appears in higher-order cell problems. The second term is

< σ i > V
that is null for the first order cell problem but is different of zero for higher order cell problems. This has an important consequence on the final form of the iterative scheme. Considering eq. (3.51), we have :

< σ i > V = 2µ * c f D i-1 -< Ω * f i > V f + < ∆ 0 * (S -S 0 ) : σ i -g i > V f (3.60)
Introducing

T i = 2µ * c f D i-1 + Ω * f i -< Ω * f i > V f , (3.61)
the integral equation reads :

σ i = T i + < L(σ i ) > V f -L(σ i ), L(σ i ) = ∆ 0 * (S -S 0 ) : σ i -g i (3.62)
The solution of this integral equation can be computed with the following iterative process :

u (k+1) = T i + < L(u (k) > V f -L(u (k) ) (3.63)
which is initialized by :

u (k=1) = T i (3.64)
When the convergence is achieved after k = N iterations, the stress is given by :

σ i = u (N ) (3.65)
and

d y (v) i = S(y) : u (N ) -g i (3.66)
The velocity field can be computed by the integration of the term d y (v) i which reads :

v i = V i + R * [d y (v i )] (3.67)
where the third order tensor R has the following components in the Fourier space :

∀ξ = 0 : R ijk = i ξ 2 [δ jk ξ i + Q ik ξ j + Q ij ξ k ] (3.68)
and R ijk = 0 for ξ = 0. Consequently, the field R * [d y (v i )] has a null volume average. The macroscopic velocity V i in (3.67) is determined by the condition v = 0 in V s .

Elimination of a rotation

In some particular cases, a rotation of the solid phase may occurs. This is possible when the problem is dissymmetric due to the geometry of the solid (flow around a triangle for instance) or when the enforcing terms in the local problem are dissymmetric. This rotation is eliminated by superposing the solution of the elementary problem :

   µ * ∆v -∇p = 0 in V f div v = 0 in V f v = W ∧ x on Γ (3.69)
where W represents a prescribed rotation of the solid. The solution of this elementary problem is computed with the FFT by considering the following source term :

f (y) = 0 in V f C ∧ y in V s (3.70)
where the couple of force C is prescribed. This couple generates the rotation W of the solid. The velocity field in V f is determined with the iterative scheme (3.62) in which g i = 0 and the source term f is given by equation (3.70). Since the problem is linear, the rotation is proportional to the applied couple of force, W = P .C, the rotation of the solid being computed with :

W =< w > Vs , w = 1 2 rot(v) ⇔ w = 1 2 ξ ∧ v (3.71)
As soon as the components of the matrix P are determined, the value of the couple of force needed to obtained the desired value of the rotation W is C = P -1 .W . The solution of this elementary problem is computed one time, is stored, and used thereafter to eliminate the rotations of the solids.

As an example, consider the case of a fluid flow trough a periodic array of cylinders with the radius R. The problem is two dimensional, the couple of force has a unique component 3.5 Validation : the fluid along parallel pores

Statement of the problem

The case of fluid flow between parallel pores is considered as a validation of the approach. The porous solid is constituted of parallel layers which are periodically distributed along the direction y 2 . A fluid flow is generated in the direction y 1 and the problem is invariant along y 3 direction. The unit cell is defined on figure 3.2. It consist of a porous layer of thickness 2h in a squared unit cell with the dimensions 2H. The problem is solved in the plane Oy 1 y 2 . Since the problem is invariant along the direction y 1 , all the local variables depends on variable y 2 which is just denoted y for simplicity. Each unit cell problems have been examined by Auriault [3]. Here the solutions are derived (various corrections for the higher order velocity fields are given from the paper of Auriault [3]). Particular attention is also paid to the stress field inside the solid phase which are given in order to make the comparison with the FFT solution. The derivation of analytic solution also allows to check the accuracy of the numerical solution with the grid refinement.

General equations

The local stress fields only depends on the coordinate y 2 ≡ y. The equilibrium of the stress fields at each order then reads :

∂σ i 12 ∂y = f i 1 (y), (3.72)
∂σ i 22 ∂y = f i 2 (y) (3.73)
The stress components σ i 12 and σ i 22 are determined by integration of the two above equations and by considering the continuity of the traction at y = ±h. The resulting stress is defined up to a constant which is determined by considering equation (3.51). In the second equation in (3.41), the local strain rate-stress relation reads in the fluid region (in the case of plane strain condition) :

d i 11 = 1 2µ * σ i 11 -σ i 22 2 -g i 11 , d i 22 = 1 2µ * σ i 22 -σ i 11 2 -g i 22 , d i 12 = σ i 12 2µ * -g i 12 (3.74)
Since all the derivatives with respect to y 1 are null, we have d i 11 = 0 and

σ i 11 = σ i 22 + 4µ * g i 11 (3.75)
and the remaining non null components of the strain rate are :

d i 22 = -(g i 11 + g i 22 ) = -tr(g i ), d i 12 = σ i 12 2µ * -g i 12 (3.76)
Note that the stress is periodic since σ 12 (y = ±H) = 0. The exact solution for the displacement v 0 1 and the shear stress σ 0 12 are represented on figure 3.3 and 3.4 respectively for h/H = 0.5 and for an applied pressure gradient J 1 = 1. The solution in the fluid phase lies in the range [-0.5, 0.5], the solution for |y| > 0.5 corresponds to the continuation of the fields within the solid region. It is observed that the numerical data greatly reproduce the exact solution. The FFT solution has been computed with the grid 64 × 64.A more complete analysis of the convergence with the resolution is provided in the next.

First corrector

In equation (3.41), the stress field σ 1 is in equilibrium with the source term f 1 whose general expression has been given in equation (3.47). In the fluid region, this source term depends on the stress field σ 0 which has only one non null component σ 0 12 and whose expression is given by (3.81). Replacing this expression in (3.47), we obtain the following expression for the components of the source term f 1 :

f 1 1 (y) = -yJ 12 in V f 0 in V s , f 1 2 (y) = -yJ 11 in V f 0 in V s (3.86)
where J 11 and J 12 represents the derivatives of J 1 with respect to x 1 and x 2 respectively. The solutions of equations (3.72) and (3.73) for i = 1 are :

σ 1 12 =      1 2 (B -h 2 )J 12 in V s 1 2 (B -y 2 )J 12 in V f (3.87) σ 1 22 =      1 2 (C -h 2 )J 11 in V s 1 2 (C -y 2 )J 11 in V f (3.88)
where the constants B and C are determined with (3.51). Particularly, considering the case i = 2 and the components 12 and 22 in equation (3.51), it gives :

< σ 1 12 > V f = 2µ * H h D 0 12 = - h 2 3 J 12 (3.89) < σ 1 22 > V f = - µ * H h D 0 11 = h 2 3 J 11 (3.90)
Introducing expressions (3.87) and (3.88) in the two above equations leads to the following expressions for the two constants B and C :

B = - h 2 3 , C = h 2 (3.91)
The stress σ 1 11 is computed with equation (3.75) :

σ 1 11 =      1 2 3y 2 - 5h 2 3 J 11 in V f ? in V s (3.92)
It must be noted that the stress σ 1 11 is undetermined in the solid region. The velocity field is now computed from eq. (3.77) and (3.78). Considering the solution for the shear stress (3.85) with (3.91), that of g 1 12 :

g 1 12 = 1 2 ∂v 0 1 ∂x 2 = 1 4µ * (y 2 -h 2 )J 12 (3.93)
and the adherence condition v 1 1 (y = ±h) = 0 in equation (3.78), we deduce that :

v 1 1 = 1 3µ * (h 2 -y 2 )yJ 12 (3.94)
Considering now equation (3.77) with :

tr(d x (v 0 )) = 1 2µ * (y 2 -h 2 )J 11 , tr(D 0 ) = - h 3 3µ * H J 11 , (3.95)
we obtain :

v 1 2 = 1 6µ * (h 2 -y 2 )yJ 11 (3.96)
The macroscopic velocities V 1 1 and V 1 2 are both equal to zero since the v 1 1 and v 1 2 are odd function with respect to the coordinate y. On figure 3.5 and 3.6 we represent the velocity v 1 2 (y) and the stress σ 1 22 (y) due to the applied double gradient of pressure J 11 = 1. Figures 3.7 and 3.8 give the solutions v 1 1 (y) and σ 1 12 (y) for J 12 = 1.

-0.5 0

Second corrector

In the equilibrium equations (3.72) and (3.73), the source terms have the following expressions :

f 1 1 (y) =        1 2 5h 2 3 -3y 2 J 111 + 1 2 h 2 3 + y 2 J 122 in V f - h 2 3 (J 111 + J 122 ) in V s , (3.97)
f 1 2 (y) =      h 2 3 -y 2 J 112 in V f 0 in V s (3.98)
which, by integration gives the following expressions for the stress field :

σ 2 12 =                  1 2 5h 2 3 y -y 3 J 111 + 1 6 h 2 y + y 3 J 122 + D |y| ≤ h - h 2 3 (J 111 + J 122 )y + E y ≥ h - h 2 3 (J 111 + J 122 )y + E ′ y ≤ -h , (3.99)
σ 2 22 =          1 3 h 2 y -y 3 J 112 + F |y| ≤ h G y ≥ h G ′ y ≤ -h (3.100)
where D, E, E ′ , F, G, G ′ are six constants. They are determined by the continuity of the traction at y ± h and with condition (3.51) in which D 1 = 0 since the velocity V 1 is null. We obtain :

E = 2h 3 3 (J 111 + J 122 ), E ′ = - 2h 3 3 (J 111 + J 122 ), D = 0, F = 0, G = 0, G ′ = 0 (3.101)
The stress component σ 2 11 is computed from equation (3.75) :

σ 2 11 =    5 3 h 2 y -y 3 J 112 in V f ? in V s (3.102)
The velocity components are determined by equations (3.77) and (3.78) :

v 2 2 = 1 8µ * (h 2 -y 2 ) 2 J 112 v 2 1 = - 1 12µ * (h 2 -y 2 )(3h 2 -y 2 )J 111 - 1 24µ * (h 2 -y 2 )(h 2 + 3y 2 )J 122 (3.103)
The components of the macroscopic velocity V 2 are :

V 2 1 = - 1 45
h 5 µ * H (7J 111 + 2J 122 ), V 2 2 = 1 15
h 5 µ * H J 112 (3.104)
The solutions are represented on figures 3.9 to 3.14.

Macroscopic filtration law

Introducing the velocities at first, second and third order in eq. (3.24) gives the following expression for the components of the macroscopic velocity :

V 1 = - h 3 3µ * H J 1 -ǫ 2 h 5 45µ * H (7J 111 + 2J 122 ), V 2 = 0 (3.105)
Note that the second corrector exhibits a velocity along direction Ox 2 but when it is introduced in equation (3.24) it gives a null component for the flux in direction Ox 2 that is obviously expected for the aligned parallel pores. Note also that the contribution of the first corrector in the computation of the flux is null because the local velocity is an odd function with the coordinate y. The correction to the Darcy law is then proportional to J 111 and J 122 . Note that term proportional to J 122 was missing in [3].

Accuracy of the FFT solution

In order to evaluate the accuracy of the FFT numerical solution we compute the relative error :

error = 1 - i=N -1 i=-N |v F F T (y n )| |v exact (y n )| (3.106)
Where 2N represent the number of wave vectors taken along each space direction in the FFT algorithm, then 2N × 2N is the resolution. The errors for the first order solution and the correctors are provided on table (3.1) for increasing values of the resolution. It is observed that a good accuracy is achieved for 128 × 128. This resolution is used for the problem of flow in a regular array of cylinders.

J 1 = 1 J 11 = 1 J 12 = 1 J 111 = 1 J 112 = 1 J 122 = 1 Resolution : v 0 1 v 1 2 v 1 1 v 2 1 v 2 2 v 2

Fluid flow through a periodic array of cylinders

In this section we computed the correctors of the Darcy law in the case of a porous solid made up of aligned cylinders periodically distributed in the Ox 1 x 2 plane. The unit cell is given on figure 3.15. The radius of the cylinder is denoted by R and the size of the unit cell is H. The problem is two dimensional. The unit cell problems are computed at the third order (second corrector to the Darcy law) with the FFT iterative scheme. No analytic solution, either numeric one, have been provided for the two first correctors in the literature and by so on, no comparison of our solutions is possible. The first order solution (Darcy problem) has been well documented using semi-analytic methods (Sangani and Acrivos 1982 [63], Drummond and Tahir 1984 [22], Wang 2001 [79]). As a consequence it is not discussed anymore. Just as for comparison with the correctors, the variation of the Darcy transverse permeability k 0 with the radius of the sphere is given on figure 3.16. The results for the first and second correctors are discussed in the next.

First Corrector

Let us recall that the enforcing terms of the first corrector depends on the solution at first order (not detailed in the present paper). The solution at first order is computed by considering the applied pressure gradient component J 1 , the case corresponding to J 2 is not computed since it is equivalent to J 1 due to the symmetry of the unit cell. The enforcing terms of the second order cell problem then depend the derivative of J 1 with respect to x 1 and x 2 as in the case of the flow along parallel pores. The first corrector has then two terms proportional to J 11 and J 12 . Again, the term proportional to J 22 would be determined from the solution related to J 11 . It must be also noted that the contribution of the first correction to the Darcy law is also null dur to the symmetry of the unit cell. However it is necessary to compute it since these fields enters into the definition of the enforcing terms at third order.

Second corrector

The solution for the second corrector are determined for the components of the macrosocpic pressure gradient J 111 , J 112 . It is not necessary to compute the solution corresponding to J 122 and J 222 due to the geometry of the unit cell (that is invariant by any rotation of an angle θ = π/2 around the direction Ox 3 ). Since the macroscopic filtration law remains unchanged if we invert the axes Ox1 and Ox2, the macroscopic flux is :

V 1 = - 1 µ * k 0 J 1 + ǫ 2 k 2 1 J 111 + ǫ 2 k 2 2 J 122 , (3.107)
V 2 = - 1 µ * k 0 J 1 + ǫ 2 k 2 1 J 222 + ǫ 2 k 2 2 J 112 (3.108)
which involves three coefficients k 0 , k 2 1 and k 2 2 . The variations of coefficients k 2 1 and k 2 2 with the radius of the cylinder are represented on figures 3.17 and 3.18. It is observed that k 2 1 is positive but k 2 2 is positive for small values of the radius R and becomes negative when R is superior to a value close to 0.07H. It is also noted that coefficients k 2 1 and k 2 2 are smaller than k 0 with a factor ∼ 20.

Conclusion

In this paper we propose a FFT-based numerical algorithm to compute the unit cell problem coming from the asymptotic homogenization analysis of Stokes flow problem in a porous solid. The resolution of the first unit cell problem provide the famous Darcy's permeability. Considering the higher order cell problems we compute the correction to the Darcy permeability which depends on the double, triple, etc gradient of the macroscopic pressure. Substituting the Gradient of macroscopic velocity to the Gradient of macroscopic pressure, the reformulation of the macroscopic law leads to a generalization of the Brinkman equation for an anisotropic porous material (the Brinkman equation being established for an isotropic material). The Correctors have been computed for two particular microstructure. First, the flow along parallel pores is considered as a Benchmark problem. The good accuracy of the FFT solution with analytic ones has been shown in the paper. Next the transverse flow through a regular array of cylinders is considered. It is shown that one coefficient is either positive or negative depending on the value of the cylinder. The method could be now applied to more realistic 3D porous microstructures such as metal foams are fibrous materials, in order to investigate the importance of the correctors to Darcy's law.

Chapitre 4

Derivation of FFT-based bounds for the homogenized elastic coefficients of polycristals

Introduction

An alternative approach to Finite Elements Methods (FEM) has been proposed in the middle of the nineties by Moulinec ans Suquet [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] for the computation of the effective properties of linear elastic periodic composites. The unit cell problem is solved by means of an iterative scheme which uses the periodic Green's tensor for the strain and could be directly applied to digital images which come from microtomography. The main advantages of this approach over the FEM is that it does not requires the introduction of a "high dimension" stiffness matrix. The FFT methods only need the storage of tables whose dimension are the number of DoF (the FEM requires the storage of the rigidity matrix whose dimension is the square of the DoF). The memory needed for solving the unit cell problem with FFT methods is then significantly reduced. The FFT based method [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] allows to expand the solution of the Lippman-Schwinger equation into Neumann series, along the lines of a method which was first introduced for composite conductors by Brown [START_REF] Wf Brown | Solid mixture permittivities[END_REF] and later by Kröner [START_REF] Kröner | Statistical continuum mechanics[END_REF]. The convergence of the method has been largely studied in the literature particularly for the problems with high contrasts between the phase elastic properties. Other formulation have been then provided to improve the convergence : the (dual) stress based formulation [START_REF] Bhattacharya | A model problem concerning recoverable strains of shape-memory polycrystals[END_REF][START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF] and the accelerated schemes [START_REF] Dj Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF][START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF][START_REF] Monchiet | A polarization-based fft iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF]. Accelerated schemes are not considered in the present study since the contrast in polycrystal is sufficiently low to use the basic strain and stress based iterative schemes. An alternative approach based on the shape functions has been developed by Bonnet [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF] to improve the FFT solutions. The shape functions are the exact expressions of the Fourier transform of the characteristic functions of the phase. Such close-form expressions are available for inclusions having ellipsoidal shapes in the book of Nemat-Nasser [START_REF] Nemat-Nasser | Micromechanics : overall properties of heterogeneous materials[END_REF], for instance. In the paper of Bonnet [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF], the method has been illustrated in the case of a composite with aligned cylindrical fibers. Later, in [START_REF] Monchiet | Combining fft methods and standard variational principles to compute bounds and estimates for the properties of elastic composites[END_REF], it has been proved that the use of the shape functions in the FFT-based iterative schemes provides rigorous bounds of the effective elastic properties of the composites. More precisely, the strain and stress based FFT solutions deliver an upper and lower bound respectively. In this chapter, we propose to apply the method to polycrystals. A lower and an upper bounds of the homogenized elastic coefficients are computed by means of FFT method combined with the use of the shape functions. The microstructure of the polycrystal is generated by Voronoi tesselations. The polycristal is constituted of single crystal which are polygons (for 2d problems) and polyhedrons (for 3d ones). The shape functions of the polygons and polyhedrons have been recently provided in the paper of Wuttke [START_REF] Wuttke | Form factor (Fourier shape transform) of polygon and polyhedron[END_REF]. We consider a heterogeneous elastic material defined by a parallelepipedic unit cell and three (two for plane strain or plane stress problems) vectors of translation invariance. The unit cell is made up of M phases whose elastic tensor (resp. compliance) is denoted by C α (resp. S α = (C α ) -1 ) for α = 1..M . Classically, the local problem involves compatibility equations, linear elastic constitutive equations, equilibrium and periodic conditions at the boundary of the considered unit cell :

           ε(x) = 1 2 (∇u(x) + ∇ t u(x)), ∀x ∈ V σ(x) = C(x) : ε(x), ∀x ∈ V div(σ(x)) = 0, ∀x ∈ V u(x) -E.x periodic σ(x).n antiperiodic (4.1)
in which the stiffness tensor C(x) (resp. the compliance) of the heterogeneous medium is given by :

C(x) = α I α (x)C α , S(x) = α I α (x)S α with : I α (x) = 1 if x ∈ V α 0 if x ∈ V -V α (4.2)
I α (x) for α = 1..M are the characteristic functions describing volumes V α which comply with α I α (x) = 1. Prescribed macroscopic strain E =< ε > V or macroscopic stress Σ =< σ > V are classically considered (the brackets denote the volume average over V ).

Variational principles based on the macroscopic elastic energy are now considered to determine a lower and an upper bound for the effective elastic coefficients of the composites. Let us introduce the strain elastic energy U(ε(u)) and the complementary elastic potential W(σ) by :

U(ε(u)) = 1 2 < ε(u) : C : ε(u) > V , W(σ) = 1 2 < σ : S : σ > V (4.3)
Since the elastic tensors C and S are positive definite, the following energy minimization principles hold for any kinematically admissible displacement u * and and statically admissible stress field σ * (see for instance [START_REF] Castaneda | Nonlinear composites[END_REF]) :

1 2 E : C hom : E ≤ U(ε(u * )) (4.4) 1 2 Σ : S hom : Σ ≤ W(σ * ) (4.5)

Discretization with Fourier series

Due to the periodicity, the strain and stress fields are expanded into Fourier series. The expansion in Fourier series of any real quantity f is denoted by :

f (x) = n=N -1 n=-N f (ξ n ) exp(iξ n .x), f (ξ) =< f (x) exp(-iξ.x) > V , (4.6)
where ξ n denote the discrete wave vectors given by

ξ n = 2πnζ, n = -N..N -1, ζ i = 1 L i (4.7)
and L 1 , L 2 , L 3 are the half of the cell along the three space directions x 1 , x 2 and x 3 .

Obviously, the problem is discretized along each space direction, this would involve the use of two index n 1 and n 2 for 2-d problems and three for 3-d one. However, only one index n is used for simplicity. Still for simplicity, we shall use the notation f n , the Fourier transform of f associated with the wave vector ξ n . Particularly, the Fourier component corresponding to n = 0 of the quantity f represents its average over the volume of the cell, f 0 =< f (x) > V . Then, the quantity ε 0 = E is known when the macroscopic strain E is prescribed to the unit cell. In this case, the macroscopic stress is Σ = σ 0 that is determined by the resolution of the unit cell problem. The homogenized elasticity tensor is determined such that Σ = C hom : E and we also denote by S hom = (C hom ) -1 the homogenized compliance tensor. Alternatively, when the macroscopic stress is applied, the macroscopic strain have to be determined. In the next section, classic variational principles are considered to obtain elementary bounds for the macroscopic elastic energy. These variational principles are based on the consideration of admissible subsets made up of compatible fields for the strain, and equilibrated fields for the stress. These admissibility conditions can be formulated by means of two suitable projection operators, P n and Q n , defined for every ξ n = 0 by :

Q n = E 1 (ξ n ) + E 3 (ξ n ), P n = E 2 (ξ n ) + E 4 (ξ n )
in which the E i (ξ) for i = 1..4 are the first four tensors of the Walpole basis [START_REF] Walpole | Elastic behavior of composite materials : theoretical foundations[END_REF]. Let us recall the six tensors of this basis :

     E 1 (ξ) = 1 2 Q(ξ) ⊗ Q(ξ), E 2 (ξ) = P (ξ) ⊗ P (ξ) E 3 (ξ) = Q(ξ)⊗Q(ξ) -E 1 (ξ), E 4 (ξ) = Q(ξ)⊗P (ξ) + P (ξ)⊗Q(ξ) E 5 (ξ) = P (ξ) ⊗ Q(ξ), E 6 (ξ) = Q(ξ) ⊗ P (ξ) (4.8)
where P and Q are given by :

P (ξ) = 1 |ξ| 2 ξ ⊗ ξ, Q(ξ) = I -P (ξ) (4.9)
and I is the second order identity tensor. As already shown in [START_REF] Monchiet | A polarization-based fft iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF], the strain field is compatible if its fourier coefficients ε n comply with ∀n = 0, Q n : ε n = 0 (4.10)

Alternatively, the equilibrium for the stress field leads to the following condition for its fourier coefficients σ n :

∀n = 0, P n : σ n = 0 (4.11)
With these notations, the subset of kinematical and statical admissible fields, denoted E N and S N respectively, are defined by :

E N = ε = n=N -1 n=-N ε n exp(iξ n .x), ε 0 = E, Q n : ε n = 0 for n = 0 S N = σ = n=N -1 n=-N σ n exp(iξ n .x), σ 0 = Σ, P n : σ n = 0 for n = 0 (4.12)
Using the variational principles (4.4) and (4.5), we deduce that two estimates for the effective elastic properties can be deduced from :

1 2 E : C hom : E ≤ 1 2 E : C ε N : E = min ε ∈ E N U(ε) 1 2 Σ : S hom : Σ ≤ 1 2 Σ : S σ N : Σ = min σ ∈ S N W(σ) (4.13)
where C ε N and S σ N are two estimates of the homogenized elastic tensor C hom = (S hom ) -1 in which the index N refers to the dimension of the subset of admissible fields for the strain and the stress. Classically, the following inequality is deduced for the elasticity tensors :

C σ N ≤ (S hom ) -1 = C hom ≤ C ε N (4.14)
with the notation C σ N = (S σ N ) -1 .

Derivation of bounds with FFT

The stationary point of U with respect to ε n can be determined by means of a FFT based iterative scheme which uses the shape functions. The terminology "shape function" will more explained in the next. The details of the minimization procedure can be found in [START_REF] Monchiet | Combining fft methods and standard variational principles to compute bounds and estimates for the properties of elastic composites[END_REF] but are note detailed here. The following recurrence relation is then used to computed the Fourier components of the strain field :

( ε n ) i+1 = ( ε n ) i -Γ 0 n : α=p α=1 C α : ( ε α n ) i (4.15)
which starts with ε n = 0 for n = 0 and ε 0 = E. In (4.15), Γ 0 n are the Fourier coefficients of the Green tensor associated to the reference material of rigidity C 0 . The fourth order tensor Γ 0 is the Green operator for the strain. Its components depends on the wave vectors and also on the elastic coefficients of a reference material C 0 . The reference material is chosen in order to provide the better rate of convergence of the iterative scheme. The method for its determination has been well documented in the literature and is not detailed here (the reader could refer to [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF][START_REF] Gw Milton | The theory of composites[END_REF][START_REF] Monchiet | A polarization-based fft iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF] for more details on these issues). Considering an isotropic elastic reference material, the Green operator can be decomposed along the Walpole basis :

∀n = 0 : Γ 0 n = 1 λ 0 + 2µ 0 E 2 (ξ n ) + 1 2µ 0 E 4 (ξ n ) (4.16)
and Γ 0 n = 0 for n = 0. In (4.16), λ 0 and µ 0 are the Lamé coefficients of the reference medium. In equation (4.15), ε α represents the product between the strain and the characteristic function I α (x) of the the phase V α :

ε α (x) = I α (x)ε(x) (4.17)
and ε α n are its Fourier coefficients. The exact expressions of ε α n for α = 1..M are then computed by the discrete convolution product between the Fourier components of the strain and the Fourier components of the characteristic functions I α (x) :

ε α n = m=N -1 m=-N I α (ξ n -ξ m ) ε m (4.18) σ α n = m=N -1 m=-N I α (ξ n -ξ m ) σ m (4.19)
Since n and m vary from -N to N -1, then n-m vary from -2N to 2N -2. It follows that in equation (4.18), I α (ξ n -ξ m ) must be computed on a double grid (dimension is 4N ×4N ) while the coefficients ε m are computed on the simple grid (dimension is 2N × 2N ). The components I α n = I α (ξ n ) must be computed before the iteration process and stored. The procedure is computationally more expensive than in [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] because the convolution product is made on the double grid while the former method uses a representation of the elasticity tensor on the simple grid. However, the representation on the double grid is necessary to compute exactly the discrete convolution product and then to ensure the bound character of the solutions. Moreover, the method accounts for the real geometry of the cell when exact Fourier coefficients of the characteristic function are used while, in [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF], the product between the elastic tensor with the strain is made with the approximation :

ε α n ∼ F F T I α (x n )F F T -1 ( ε n ) (4.20)
in which I α (x n ) are the values of I α (x) computed at the nodes of a regular grid in the real space. The shape functions account for the real geometry of the unit cell when the exact expressions of these functions could be derived. For example, for the case of a composite containing fibers with circular cross section, the Fourier transform of the characteristic function can be analytically computed as function of the Bessel function. In this case, there is no approximations on the geometry, while, in the former method [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF], the real image was replaced by a pixelized image which approximates the inclusion geometry by means of rectangular elements. Note also that finite element methods also approximate the real geometry since they use triangular elements (or other kinds of elements) in order to mesh the inclusion. Truly curved objects are then approximated by polygon meshes. The lower bound for the elastic tensor is computed with the stress based iterative scheme. This scheme has been formulated by [START_REF] Bhattacharya | A model problem concerning recoverable strains of shape-memory polycrystals[END_REF][START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF], it uses the compliance and the dual, stress based, Green operator :

( σ n ) i+1 = ( σ n ) i -∆ 0 n : α=p α=1 S α : ( σ α n ) i (4.21)
which is initialized with σ n = 0 for any n = 0 and σ 0 = Σ. In (4.21), ∆ 0 n are the Fourier coefficients of the Green tensor for the stress, defined by :

∆ 0 n = C 0 -C 0 : Γ 0 n : C 0 = 2µ 0 (3λ 0 + 2µ 0 ) λ 0 + 2µ 0 E 1 (ξ n ) + 2µ 0 E 3 (ξ n ) (4.22)
for n = 0 and ∆ 0 n = 0 for n = 0. At each step of the stress based iterative scheme, we need to compute σ α n obtained as the convolution product between the shape functions and the stress. σ α n is computed from (4.18) in which the strain ε is replaced by the stress σ at both sides of the equality. The convergence test used for the strain based iterative scheme is based on the equilibrium for stress. The dual scheme use a convergence condition based on the compatibility of the strain : strain scheme : P n :

α=p α=1 C α : ( ε α n ) i < ǫ stress scheme : Q n : α=p α=1 S α : ( σ α n ) i < ǫ (4.23)
where the precision ǫ = 10 -4 is considered in the applications.

It is worth noted that the iterative schemes (4.15) and (4.21) are formally equivalent to that introduced in [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF][START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] and [START_REF] Bhattacharya | A model problem concerning recoverable strains of shape-memory polycrystals[END_REF][START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF]. They use the Green operator that is explicitly determined by the wave vectors when it is represented in the Fourier space. However, the difference lies in the computation of the local elastic law. The original scheme uses the characteristic functions I α (x) evaluated on a regular grid in the real space while the modified method used the shape functions I α (ξ) evaluated along a regular grid in the Fourier space.

Properties of the shape function

Let us recall that the shape functions are the exact expressions of the Fourier transform of the characteristic functions :

I α (ξ) = 1 V V α exp(-iξ.x)dx (4.24)
Particularly, when ξ = 0, we have :

I α (0) = V α V = c α (4.25)
where c α represents the volume fraction of the phase α.

The characteristic functions of the phase comply to the equality :

α I α (x) = 1 (4.26)
As a consequence, the shape functions verify the equality :

α I α (ξ) = 0 (4.27)
for any ξ = 0. Indeed, the Fourier transform of 1 is null whatever ξ excepted for ξ = 0 for which it takes the value 1. For ξ = 0, the shape functions take the value c α and α c α = 1. Consider a two phase composite made-up of inclusions surrounded by a matrix. The phase α = 1 corresponds to the matrix and α = 2 to the inclusion. Due to the properties of the shape function, we deduce that :

I 1 (ξ) = -I 2 (ξ) (4.28)
for any ξ = 0. So, if the shape function of the inclusion is known, that of the matrix is also known. Only the storage of I 2 (ξ) is required for the implementation of the FFT method. Consider for instance a composite with cuboidal inclusions. The size of the inclusion is 2a, 2b, 2c along the directions x 1 , x 2 and x 3 . The unit cell is a cube and its size is L.

It contains one cuboidal inclusion centered at the origin. The computation of the shape function I 2 is straightforward since the triple integral can be decomposed has the product of three simple integrals :

I α (ξ) = 1 V a -a exp(-iξ 1 x 1 )dx 1 × b -b exp(-iξ 2 x 2 )dx 2 × c -c exp(-iξ 3 x 3 )dx 3 (4.29)
It leads to :

I 2 (ξ) = 8abc L 3 sinc(aξ 1 )sinc(bξ 2 )sinc(cξ 3 ) (4.30)
for all ξ = 0 and where sinc is the cardinal sinus. For ξ = 0, the shape function I 2 (ξ) is equal to the volume fraction of inclusions, that is 8abc/L 3 . Note that this value, can be also recovered by taking the limits ξ 1 → 0, ξ 2 → 0 and ξ 3 → 0, in the above expression.

In the 2d case, the rectangular inclusion in a squared unit cell, the shape function is :

I 2 (ξ) = 4ab L 2 sinc(aξ 1 )sinc(bξ 2 ) (4.31)

Application to 2d-polycristals 4.3.1 Local elastic law

We consider a 2d-polycristal made-up of M single cubic-crystals. The elastic tensor of the crystal α is denoted C α . For convenience, it is more suitable to read the local elasticity rule with a matrice representation. For instance, in the 2d case, the elastic law written in the basis B α oriented along the axis of the crystal is :

  σ 1 σ 2 σ 6   Bα =   C α 11 C α 12 C α 16 C α 12 C α 22 C α 26 C α 16 C α 26 C α 66   Bα   ε 1 ε 2 ε 6   Bα (4.32)
Where the following notation has been used for the components of the strain and the stress :

  σ 1 σ 2 σ 6   =   σ 11 σ 22 √ 2σ 12   ,   ε 1 ε 2 ε 6   =   ε 11 ε 22 √ 2ε 12   (4.33)
and the components C α ij are related to that of C α ijkl by : 

 C α 11 C α 12 C α 16 C α 12 C α 22 C α 26 C α 16 C α 26 C α 66   Bα =   C α 1111 C α 1122 √ 2C α 1112 C α 1122 C α 2222 √ 2C α 2212 √ 2C α 1112 √ 2C α 2212 2C α 1212   Bα (4.34)
The elasticity law in the crystal is assumed to be cubic. As a consequence the material has three independent elastic coefficients λ, µ 1 , µ 2 and :

  C α 11 C α 12 C α 16 C α 12 C α 22 C α 26 C α 16 C α 26 C α 66   Bα =   λ + 2µ 1 λ 0 λ λ + 2µ 1 0 0 0 2µ 2   Bα (4.35)
In the FFT method, the strain and the stress are computed in a global frame B. The components of the elastic tensor must be computed in this global frame. Let us denote by ε ′ i and σ ′ i the components of the strain and stress written in the global frame. The orientation of the crystal in the global frame is given by the angle θ. The base change relations giving the components of the strain ε i as function of ε ′ i are :

ε 1 = c 2 ε ′ 1 + s 2 ε ′ 2 - √ 2csε ′ 6 (4.36) ε 2 = s 2 ε ′ 1 + c 2 ε ′ 2 + √ 2csε ′ 6 (4.37) ε 6 = - √ 2cs(ε ′ 2 -ε ′ 1 ) + (c 2 -s 2 )ε ′ 6 (4.38)
in which c = cos(θ) and s = sin(θ). The relations giving the components σ ′ i as function of σ i are :

σ ′ 1 = c 2 σ 1 + s 2 σ 2 + √ 2csσ 6 (4.39) σ ′ 2 = s 2 σ 1 + c 2 σ 2 - √ 2csσ
 σ ′ 1 σ ′ 2 σ ′ 6   B =   C α 11 ′ C α 12 ′ C α 16 ′ C α 12 ′ C α 22 ′ C α 26 ′ C α 16 ′ C α 26 ′ C α 66 ′   B   ε ′ 1 ε ′ 2 ε ′ 6   B (4.42)
In which the components C α ij ′ are given by :

C α 11 ′ = λ + 2µ 1 + 4(µ 2 -µ 1 )c 2 s 2 (4.43) C α 12 ′ = λ + 4(µ 1 -µ 2 )c 2 s 2 (4.44) C α 16 ′ = 2 √ 2(µ 1 -µ 2 )cs(s 2 -c 2 ) (4.45) C α 22 ′ = C α 11 ′ (4.46) C α 26 ′ = -C α 16 ′ (4.47) C α 66 ′ = 2µ 2 + 8(µ 1 -µ 2 )c 2 s 2 (4.48)
It is readily observed that the components of C α ij ′ are the same that C α ij as soon as µ 1 = µ 2 . This corresponding to the particular case of an isotropic elastic medium.

The shape function of a polygon

The representative cell of the polycristal is generated by Voronoi tessellations. Each single crystal is then represented by a polygon. The number of edges of the polygon is arbitrary. In a given microstructure, the polygons must contains 3,4,5,... edges. Consider a polygon and let us denote the positions of the corners by r 1 , r 2 , r 3 , ..., r J , the corners being numbered in counter-clockwise direction. The shape function I α and the area S α of the α-polygon are given by the expressions in [START_REF] Hl Nguyen | Conductivity of composites with multiple polygonal aggregates, theoretical estimates and numerical solutions from polarization series[END_REF] :

I α (ξ) = - i S e 3 ∧ ξ |ξ| 2 . J j=1 (r j -r j-1 )sinc ξ. (r j -r j-1 ) 2 exp -iξ. r j + r j-1 2 (4.49) S α = 1 2 e 3 . J j=1 r j-1 ∧ r j (4.50)
with the convention r 0 = r J , |ξ| is the norm of ξ, e 3 is the normal unit vector to the working plane (e 1 , e 2 ). In the above equation S = L 1 L 2 is the area of the squared unit cell.

It is noted that this expression is preferred over another equivalent expression existing in literature for the regularity of the sinc function. In the latter analytical expression, the denominator can vanish at some values of the wave vectors, requiring attention when implementing these formula. Indeed, sinc(x) tends to zero when x also tends to zero but, numerically, when x = 0, sinc(x) leads to singularity. The limit must be correctly computed when the term ξ.

(r j -r j-1 ) 2 is null. Consider as an illustration the case of a rectangular polygon. The rectangle is centered at the origin and the position of the four corners are : That is exactly the expression obtained in section (4.2.4). The implementation of the shape functions is now checked by computing the inverse FFT of I(ξ). For the polygon, the exact inverse Fourier transform of I(ξ) is the characteristic function of the polygone that is equal to 1 inside the polygon and is null outside. Obviously, I(ξ) is computed along the discrete wave vectors ξ n . The inverse FFT is used to compute F F T -1 I(ξ), that is closed to the characteristic function of the polygon when N is sufficiently large but it is not rigorously equal at finite value of the resolution. The equality I(x n ) = F F T -1 ( I(ξ n )) is valid only when the resolution tends to infinity. As a first example, we consider a triangle, its three corners have the coordinates r 1 = (-0.25, 0.

Illustration

A representative cell of the polysristal is obtained with Voronoï tesselation and is represented on figure 4.5. The orientation of each crystal, defined by the angle θ, is randomly chosen in the interval [0, π]. The elastic moduli of the cubic crystal are µ 1 = 1, µ 2 = 2 and λ = 1. The calculations are performed on 40 reconstructed unit cells. The average value of the effective elastic moduli are λ hom and µ hom . The variations of the homogenized elastic shear modulus with the number of wave vectors are provided on Fig. 4.6. Three solutions are provided on this figure, the strain and the stress based FFT bounds computed with the shape functions, the solution obtained with the original FFT iterative scheme of Moulinec and Suquet [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] without using the shape function. It is observed that the bounds have a uniform convergence and the solution of obtained with the original scheme of Moulinec and Suquet is comprised between the two bounds.

Extension to 3d-polycristals

In this section we propose to extend the method based on the shape function to the case of 3D polycrystals.

Local elastic law

We still use the modified Voigt notations to put the strain-stress relation in a matricial form : R ij = p -q 0 q p 0 0 0 1 (4.58)

        ε 1 ε 2 ε 3 ε 4 ε 5 ε 6         Bα =         C α
with p = cos(θ) and q = sin(θ). The relations between the components of the strain in basis B α and B is classically given by :

ε ′ ij = R ip R jq ε pq (4.59)
When arranging the components of the strain in a vector which use the modified Voigt notation, we obtain :

        ε ′ 1 ε ′ 2 ε ′ 3 ε ′ 4 ε ′ 5 ε ′ 6         = R 3 (θ)         ε 1 ε 2 ε 3 ε 4 ε 5 ε 6         (4.60)
where the components of the matrix R 3 (θ) are :

R 3 (θ) =         p 2 q 2 0 0 0 - √ 2pq q 2
p 2 0 0 0 √ 2pq 0 0 1 0 0 0 0 0 0 p -q 0 0 0 0 q p 0 √ 2pq -√ 2pq 0 0 0

p 2 -q 2         (4.61)
The 6 × 6-matrix corresponding to a rotation of an angle θ around the axis x 2 is :

R 2 (θ) =         p 2 0 q 2 0 √ 2pq 0 0 1 0 0 0 0 q 2 0 p 2 0 - √ 2pq 0 0 0 0 p 0 -q - √ 2pq 0 √ 2pq 0 p 2 -q 2 0 0 0 0 q 0 p         (4.62)
and for a rotation of an angle θ around the axis x 1 the 6 × 6-matrix is :

R 1 (θ) =         1 0 0 0 0 0 0 p 2 q 2 √ 2pq 0 0 0 q 2 p 2 - √ 2pq 0 0 0 - √ 2pq √ 2pq p 2 -q 2 0 0 0 0 0 0 p -q 0 0 0 0 q p         (4.63)
Now, the components of the elastic 6 × 6-matrices C α′ and C α is :

C α′ = R T C α R T (4.64)
with : R = R 3 (θ 3 )R 2 (θ 2 )R 1 (θ 1 ) (4.65)

The shape function of a polyhedron

Each crystal of the 3d-polycristal is represented by a polyhedron. It is defined by its K faces denoted Γ k for k = 1...K. Each face Γ k is a J k -polygon given by the simple polygonal vertex chain r k,1 ; ...; r k,J , and by the normal unit vector n k pointing towards the outside of the polyhedron. The shape function of the α-polyhedron is ∀ξ = 0 :

I α (ξ) = - i V K k=1 ξ. n k |ξ| 2 χ k (ξ) (4.66)
where V denotes the volume of the unit cell and :

χ k (ξ) = -i n k ∧ ξ ξ.q k .ξ . J j=1
(r k,jr k,j-1 )sinc ξ.

(r k,jr k,j-1 ) 2 exp -iξ. r k,j + r k,j-1 2 (4.67)

where the convention r k,0 = r k,J must be used. The two order tensor q k is the projector onto the plane normal to the unit vector n k :

q k = I -n k ⊗ n k (4.68)
Note that χ k (ξ) is the 2d-shape function of the polygon k without the area S of the 2dunit cell. Indeed, when the polygon if embedded in a squared unit cell, the shape function is χ k (ξ)/S where S denotes the area of the unit cell. The volume of a polyhedron can be conveniently computed from :

V α = 1 6 K k=1 . J j=1 n k .(r k,j-1 ∧ r k,j )( n k .r k,j ) (4.69)
For instance, consider a cuboidal polyhedron centered at x = 0. Its dimensions are a, b, c along the axes x 1 , x 2 and x 3 respectively. The cuboidal is constituted of 6 faces having the outward normal unit vectors ±e 1 , ±e 2 , ±e 3 . Consider the face k = 1, defined by its outward normal unit vector e 1 .

Conclusion

In this chapter we provide rigorous bounds for 2d and 3d elastic polycrystals. The effective elastic coefficients are computed from reconstructed unit cell using standard methods based one Vornoi tesselations. In each single crystal the elastic law is assumed to be cubic. The effective elastic coefficients are computed by using the FFT methods combined with the use of the shape functions. The method has been introduced by [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF] for composites with ellipsoidal inclusions. In this chapter, the approach has been extended to the case of polycrystal by taking advantages of the exact expressions of the shape function of polygon and polyhedron recently provided in the paper of Wuttke [START_REF] Wuttke | Form factor (Fourier shape transform) of polygon and polyhedron[END_REF]. The approach could be easily considered for hexagonal and tetragonal local elasticity rule.

Conclusion générale et perspectives

In this work, we propose new results about the determination of the macroscopic properties of heterogeneous materials. We determine effective elastic properties and also the heat and mass transfer properties. Despite their differences, those problems are formulated by a Lippmann-Schwinger (LS) equation that is solved with Neumann series. New methods have been developed to estimate or to compute numerically the series.

The first two chapters of the thesis concerns the effective conductivity and elastic stiffness of composites. We have formulated a general class of LS equations of polarization (or eigenstress for elasticity problem) and established bounds for spectral radius and norm of the associated operators. These results allow the derivation of optimal equations and best estimates for the effective properties. The estimates are based on truncated Neumann series and approximation of its residual. The relations with other statistical quantities in Fourier space, the structure factors of order n ≥ 1 are also explicitly derived. Comparisons with literature examples on hard sphere inclusions and FFT numerical solutions show that estimates based on the first two structure factors already yield satisfactory results. Given the generality of the approach, it is possible to deal with more complex microstructures, for example those involving polygonal and cracks geometry. Such microstructures have numerous applications in practice including concrete, soil and polymer materials. Those directions will be explored in the future works.

In chapter 3, we compute numerically the corrective terms to the Darcy equation. At the local scale, the fluid flow is described by the Stokes equation and, by homogenization, by the Darcy equation at the macroscopic scale when the hypothesis of strict separation of scales is used (null value of the scale factor). At finite but non null value of the scale factor, the corrective terms to the Darcy equation depend on the double, triple, etc, gradient of pressure. A FFT iterative scheme has been proposed to compute the higher order cell problems coming from the asymptotic series expansion method. The results are in good agreement with the exact solution for the flow along parallel pores. The method has been next applied in the case of a periodic network of aligned rigid cylinders, the results show the relative importance of the corrective terms comparatively to the Darcy one. The FFT code has been developed in the 2d case and would be extended in the 3d case to provide more realistic values of the corrective terms. The method could then be applied to fibrous media or metal foams.

In the last chapter, we provide new FFT solutions for the effective properties of poly-cristals. The computations are performed on 2d and 3d voronoi unit cells. The method based on the shape function is used in the FFT method to deliver exact bounds for the effective elastic coefficients. The shape functions, which are the exact expressions of the Fourier transform of the phase, account for the exact geometry of the microstructure. Particularly, for a polycristal, the geometry of the microstructure is accounted by computing exact the Fourier transform of the polygons (for 2d microstructure) and of the polyhedrons (for the 3d polycristal). An upper and a lower bounds have been determined in the case of a polycristal made up of cubic single crystals. The method could be extended to account for another material symmetries such has hexagonal or tetragonal crystals. Also the case of non linear behavior could be considered.

Chapitre 1 Conductivity

 1 estimates of spherical-particle suspensions based on triplet structure factors 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Homogenization of the periodic thermal conduction problem and approximation method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Class of integral equations based on polarization . . . . . . . . . . . 1.2.2 Estimation based on Neumann series expansion and integral equation 1.2.3 Cubic and random distribution of spherical particles . . . . . . . . . 1.2.4 Integral evaluation of C 2 for random distribution with S (3) approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.5 Sample generation and direct computation of C 2 for random distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Numerical applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Cubic lattice arrangements . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Randomly distributed spheres . . . . . . . . . . . . . . . . . . . . . 1.4 Final remarks and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Statement of the problem . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 General equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.3 First order solution . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.4 First corrector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.5 Second corrector . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.6 Macroscopic filtration law . . . . . . . . . . . . . . . . . . . . . . . 3.5.7 Accuracy of the FFT solution . . . . . . . . . . . . . . . . . . . . . 3.6 Fluid flow through a periodic array of cylinders . . . . . . . . . . . . . . . 3.6.1 First Corrector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.2 Second corrector . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapitre 4 Derivation of FFT-based bounds for the homogenized elastic coefficients of polycristals 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Derivation of bounds with FFT schemes . . . . . . . . . . . . . . . . . . . 4.2.1 The cell problem for linear elastic composites with prescribed macroscopic strain or stress . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Discretization with Fourier series . . . . . . . . . . . . . . . . . . . 4.2.3 Derivation of bounds with FFT . . . . . . . . . . . . . . . . . . . . 4.2.4 Properties of the shape function . . . . . . . . . . . . . . . . . . . . 4.3 Application to 2d-polycristals . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Local elastic law . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 The shape function of a polygon . . . . . . . . . . . . . . . . . . . . 4.3.3 Illustration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Extension to 3d-polycristals . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Local elastic law . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2 The shape function of a polyhedron . . . . . . . . . . . . . . . . . . 4.4.3 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1. 1

 1 Integral over a sphere surface . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Radial distribution function of hard sphere in equilibrium at f = 0.4. The inset is a snapshot of a typical sample . . . . . . . . . . . . . . . . . . . . 1.3 Comparison the effective conductivity (k e /k m ) of FCC lattice structure between the approximation OS and FFT solution with ratio k i /k m = 100. 1.4 Comparison the effective conductivity (k e /k m ) of randomly distributed spheres with ratio k i /k m = ∞. Comparison of the second order OS-2 estimations and Torquato's expression (TO), Kim and Torquato (KT), the simulation data of Bonnecaze and Brady (BB) and the Clausius-Mossoti (CM) approximation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Isolines of spectral radius ρ(B) (left) and norm B (right) as functions of α and β. The results are obtained for the case where the two materials are of the same Poisson ratio ν = 0.3 and stiffness ration ε = 3. The optimal values by the two methods are respectively (α, β) = (0.375, 0.250) and (α, β) = (0.338, 0.349) . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Isolines of spectral radius ρ(B) (left) and norm B (right) as functions

3. 1

 1 Variations of P 33 = W 3 /C 3 as function of the radius of the cylinder. . . . . Liste des tableaux 1.1 Comparison the effective conductivity (k e /k m ) of SC lattice structure between the approximation schemes based on series OS, JS and ES and the FFT solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 The effective conductivity (k e /k m ) of randomly distributed spheres. Com-

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Comparison of the results α, β issued from three methods of optimization. Notations : OR for optimization based on spectral radius of B, ON for optimization based on norm of B and OD for optimization based on the direct estimation of spectral radius of B ′ . . . . . . . . . . . . . . . . . . . . 2.2 Normalized effective bulk modulus κ e /κ 0 vs inclusion volume fraction f of FCC. Elastic properties of the constituents are ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 100. The results are computed by first order estimates which coincide with HS estimates, second order estimates of the three methods (OR,ON and OD) and the numerical method FFT at convergence. . . . . . . . . . . . . 2.3 Effective shear moduli µ e /µ 0 vs inclusion volume fraction f of FCC. Elastic properties of the constituents are ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 100. . . . . . . 3.1 Relative error between the FFT and the exact solutions. . . . . . . . . . . 4.1 Average values of the effective shear modulus and effective compressibility computed for the 3d-polycristal as function of the resolution. Comparison between the lower bound (LB) and the upper bound (UB) and the solutions obtained with the orifinal scheme of Moulinec and Suquet (M &S). . . . . .

Figure 1 . 1 -

 11 Figure 1.1 -Integral over a sphere surface

Figure 1 . 2 -

 12 Figure 1.2 -Radial distribution function of hard sphere in equilibrium at f = 0.4. The inset is a snapshot of a typical sample

Figure 1 . 3 -

 13 Figure 1.3 -Comparison the effective conductivity (k e /k m ) of FCC lattice structure between the approximation OS and FFT solution with ratio k i /k m = 100..

Figure 1 . 4 -

 14 Figure 1.4 -Comparison the effective conductivity (k e /k m ) of randomly distributed spheres with ratio k i /k m = ∞. Comparison of the second order OS-2 estimations and Torquato's expression (TO), Kim and Torquato (KT), the simulation data of Bonnecaze and Brady (BB) and the Clausius-Mossoti (CM) approximation.

Figure 2 . 1 -Figure 2 . 2 -

 2122 Figure 2.1 -Isolines of spectral radius ρ(B) (left) and norm B (right) as functions of α and β. The results are obtained for the case where the two materials are of the same Poisson ratio ν = 0.3 and stiffness ration ε = 3. The optimal values by the two methods are respectively (α, β) = (0.375, 0.250) and (α, β) = (0.338, 0.349)

ν 1 =Table 2 . 1 -

 121 ν 0 = 0.3 0.0040 -0.0031 0.0050 0.0027 -0.0012 0.0033 µ 1 /µ 0 = 0.01 ν 1 = 0.1, ν 0 = 0.3 0.0525 -0.0215 0.0528 0.0430 0.0037 0.0159µ 1 /µ 0 = 0.1 ν 1 = 0.4, ν 0 = 0.3 0.923 -1.000 0.250 0.259 0.353 0.228 µ 1 /µ 0 = 1 ν 1 = 0.4, ν 0 = 0.3 Comparison of the results α, β issued from three methods of optimization. Notations : OR for optimization based on spectral radius of B, ON for optimization based on norm of B and OD for optimization based on the direct estimation of spectral radius of B ′ .

Figure 2 . 3 -

 23 Figure 2.3 -Different arrangements of non-overalapping spheres in cubic unit cell. From left to right : simple cubic (SC), body centered cubic (BCC), face centered cubic (FCC) and random distribution (RD). The random distribution is generated by Event Driven Molecular Dynamics method.

  90) and some intermediate results in Appendix B. Substituting n = 1 and β 1 in (2.88) into (2.82) and making use of 2α 0 = 1 and 3β 0 = 1, we can recover the Hashin-Shtrikman (HS) estimate for κ e regardless the chosen values of α and β κ e κ 0 = 1 + f

Figure 2 . 4 -

 24 Figure 2.4 -Normalized effective bulk modulus κ e /κ 0 vs inclusion volume fraction f of BCC array. Elastic properties of the constituents are ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 10. The results are computed by first order estimates which all coincide with HS estimates, second order estimates of the three methods (OR,ON and OD) and the numerical method FFT at convergence.

Figure 2 . 5 -

 25 Figure 2.5 -Normalized effective shear moduli µ e /µ 0 vs inclusion volume fraction f of BCC array. Elastic properties of the constituents are ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 10. The results are computed by first order estimates, second order estimates of the three methods (OR,ON and OD) and the numerical method FFT at convergence.

Figure 2 . 6 -

 26 Figure 2.6 -Normalized effective bulk modulus κ e /κ 0 vs inclusion volume fraction f of FCC array. Elastic properties of the constituents are ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 10. The results are computed by first order estimates, second order estimates of the three methods (OR,ON and OD) and the numerical method FFT at convergence.

Figure 2 . 7 -

 27 Figure 2.7 -Normalized effective shear moduli µ e /µ 0 vs inclusion volume fraction f of FCC array. Elastic properties of the constituents are ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 10.

Figure 2 . 8 -

 28 Figure2.8 -Normalized effective bulk modulus κ e /κ 0 vs inclusion volume fraction f for random distribution of rigid spheres (µ 1 /µ 0 = κ 1 /κ 0 = ∞). The solutions of the present work (OD-1, OD-2) are compared with the results of[START_REF] Torquato | Effective stiffness tensor of composite media. I. Exact series expansions[END_REF][START_REF] Torquato | Effective stiffness tensor of composite media : II. Applications to isotropic dispersions[END_REF]. The first order estimates coincide with HS estimates.

Figure 2 . 9 -

 29 Figure 2.9 -Normalized effective bulk modulus κ e /κ 0 vs inclusion volume fraction f with random distributions of spherical voids (µ 1 /µ 0 = κ 1 /κ 0 = 0).

C 3 1 0Figure 3 . 1 -

 3131 Figure 3.1 -Variations of P 33 = W 3 /C 3 as function of the radius of the cylinder.

2 Figure 3 . 2 -

 232 Figure 3.2 -Unit cell of the periodic porous solid with parallel pores.

Figure 3 . 3 -Figure 3 . 4 -

 3334 Figure 3.3 -Variations of the velocity v 0 1 (y) with the coordinate y for J 1 = 1.

Figure 3 . 5 -Figure 3 . 6 -Figure 3 . 7 -Figure 3 . 8 -

 35363738 Figure 3.5 -Variations of the velocity v 1 2 (y) with the coordinate y for J 11 = 1.

Figure 3 . 9 -Figure 3 . 10 -Figure 3 . 11 -Figure 3 . 12 -Figure 3 . 13 -Figure 3 . 14 -

 39310311312313314 Figure 3.9 -Variations of the velocity v 2 1 (y) with the coordinate y for J 111 = 1.

3. 6 .2HFigure 3 . 15 -

 6315 Figure 3.15 -Definition of the unit cell.

Figure 3 . 16 -

 316 Figure 3.16 -Variation of Darcy's permeability k 0 with the radius of the cylinder.

Figure 3 . 17 -

 317 Figure 3.17 -Variations of permeability k 2 1 with the radius of the cylinder.

Figure 3 . 18 -

 318 Figure 3.18 -Variations of permeability k 2 2 with the radius of the cylinder

4. 2

 2 Derivation of bounds with FFT schemes 4.2.1 The cell problem for linear elastic composites with prescribed macroscopic strain or stress

Figure 4 . 1 -

 41 Figure 4.1 -Orientation θ of the single crystal in the global frame

Figure 4 . 2 -

 42 Figure 4.2 -Numbering of the four corners of the rectangle.

Figure 4 . 3 -Figure 4 . 4 -

 4344 Figure 4.3 -Inverse FFT of the shape function of the triangle for N = 128 wave vectors.

Figure 4 . 5 -

 45 Figure 4.5 -unit cell of the Voronoi periodic structure 2D

Figure 4 . 6 -

 46 Figure 4.6 -Variation of the effective elastic shear modulus µ hom as function of the resolution.

  by 21 independent coefficients. Now, the local elasticity rule is expressed in a global frame B. To this end, the three Euler angles θ 1 , θ 2 and θ 3 are introduced. Let us first consider the rotation of an angle θ around the axis x 3 . The change of basis matrix is :

Figure 4 . 7 - , r 1 , 4 = r 1

 47141 Figure 4.7 -face k = 1 of the cuboidal polyhedron.

Figure 4 . 9 -

 49 Figure 4.9 -Variations of the effective shear modulus µ hom as function of the resolution.

Table 1 .

 1 2 -The effective conductivity (k e /k m ) of randomly distributed spheres. Comparison of the second order OS-2 estimations (with and without S (

	1.4. Final remarks and discussion

  .80) Solving (2.80) for τ Ω and substituting back into (2.73), we obtain a new expression for C e

Table 2 .

 2 2 -Normalized effective bulk modulus κ e /κ 0 vs inclusion volume fraction f of FCC. Elastic properties of the constituents are ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 100. The results are computed by first order estimates which coincide with HS estimates, second order estimates of the three methods (OR,ON and OD) and the numerical method FFT at convergence.

		2 OR-2 OD-2 FFT
	0.10 1.1784	1.1754 1.1791 1.1795 1.1907
	0.20 1.4009	1.3950 1.4024 1.4031 1.4254
	0.30 1.6863	1.6769 1.6898 1.6911 1.7293
	0.40 2.0659	2.0536 2.0758 2.0785 2.1424
	0.50 2.5952	2.5872 2.6287 2.6357 2.7533
	0.60 3.3850	3.4182 3.5059 3.5276 3.8286
	0.70 4.6900	4.9483 5.1744 5.2562 7.2886

Table 2 .

 2 3 -Effective shear moduli µ e /µ 0 vs inclusion volume fraction f of FCC. Elastic properties of the constituents are ν 1 = 0.4, ν 0 = 0.3, µ 1 /µ 0 = 100.

	f	Order-1 ON-2 OR-2 OD-2 FFT
	0.10 1.2338	1.2348 1.2350 1.2355 1.2523
	0.20 1.5365	1.5433 1.5437 1.5457 1.5819
	0.30 1.9347	1.9667 1.9675 1.9745 2.0417
	0.40 2.4719	2.5742 2.5770 2.5970 2.7232
	0.50 3.2230	3.4793 3.4896 3.5381 3.8078
	0.60 4.3281	4.8812 4.9222 5.0277 5.5844
	0.70 6.0830	7.2978 7.4748 7.7229 10.7332

Table 3 .

 3 1 -Relative error between the FFT and the exact solutions.

	1

  × 8 1.526 1.4347 1.489 1.7178 1.7186 1.7168 16 × 16 × 16 1.5121 1.453 1.4795 1.7195 1.7189 1.712 32 × 32 × 32 1.5054 1.463 1.478 1.7202 1.7205 1.7148 64 × 64 × 64 1.5018 1.4682 1.4811 1.7205 1.7207 1.7178 96 × 96 × 96 1.5007 1.47 1.482 1.7207 1.7207 1.7182 128 × 128 × 128 1.5001 1.4709 1.4823 1.7208 1.7203 1.7173 192 × 192 × 192 1.4995 1.4715 1.4826 1.7208 1.7203 1.7177 256 × 256 × 256 1.4992 1.472 1.4828 1.7207 1.7208 1.7175 Table 4.1 -Average values of the effective shear modulus and effective compressibility computed for the 3d-polycristal as function of the resolution. Comparison between the lower bound (LB) and the upper bound (UB) and the solutions obtained with the orifinal scheme of Moulinec and Suquet (M &S). .

	Resolution	UB	µ hom LB	M&S	UB	k hom LB	M&S
	8 × 8

Remerciements Acknowledgement

This project is funded by the National Foundation for Science and Technology Develop-

Isotropic case

In the case of isotropy, all tensors of odd order are null and the tensors of even number can be expressed in term of the Kronecker symbol δ ij :

[δ ij δ kl + δ ik δ jl + δ il δ jk ], ... (3.26) and the relation (3.25) becomes :

It can be noted that J ikk is the Laplacian of J i where the derivatives of the Laplacian are made with respect to the coordinates x. We note J ikk = ∆ x J i . Considering that :

the Brinkman equation can be recovered by eliminating J ikk in eq. (3.27) using eq. (3.28) and by neglecting all the terms of order ǫ n with n ≥ 4 :

Anisotropic and centrosymetric porous microstructure

Assume that the microstructure of the porous solid has the centrosymmetry, i.e. the macroscopic law is invariant if we replace x by -x. In eq. (3.25), V i is then replaced by -V i , J i by -J i , J ij remains unchanged and J ijk is replaced by -J ijk . Eq. (3.25) remains unchanged if the term proportional to "J ij " in eq. (3.25) vanishes. Consequently the components of the third order tensor K 1 are null. Moreover, from eq. (3.25) we have also :

which, by invertion of K 0 ij , leads to :

where H 0 ij is the inverse of K 0 ij . Introducing this result in eq. (3.25) and truncating at the second order, we obtain :

in which K 2 ijkl is defined by :

Eq. (3.32) is the generalization of the Brinkman equation written for an anisotropic centrosymmetric material.

Also the velocity depends only on the coordinate y, consequently :

)

It is observed that the velocity v i 2 is in fact fully determined by the mass balance equation.

First order solution

Only the component J 1 of the macroscopic pressure gradient is applied along the direction x 1 (it is obviously not physically possible to apply the component J 2 ). The case of fluid flow in direction x 3 is equivalent to that in x 1 direction that's why the problem in the plane Ox 1 x 2 is only considered. The velocity is given by the well-known Poiseuille solution :

And the average values are :

The shear stress within the fluid region is :

In the solid region, the drag force f 0 1s is :

By integration of :

we obtain the expression of the shear stress within the solid region :

where A is a constant. It is determined by the continuity of the shear stress at y = ±h, this leads to :

The corresponding 2d-shape function is :

Let k = 2 be the face having for outward normal unit vector -e 1 . The corresponding 2d-shape function is :

The sum of the faces k = 1 and k = 2 in (4.66) gives :

By adding the faces ±e 2 and ±e 3 in the sum, we finally obtain the expression :

which has been obtained in section (4.2.4) by a direct computation of the Fourier transform.

Application

The homogenized elastic coefficients of the polycristal are now computed with the method based on the shape functions. A unit cell of the polycrystal is shown on figure

Résumé

Ce travail propose de nouvelles contributions aux méthodes d'homogénéisation avec des applications aux composite, aux polycristaux et aux milieuxporeux. Les propriétés effectives sont déterminées en résolvant un problème élémentaire sur la cellule unitaire que l'on peut reformuler avec l'équation de Lippmann-Schwinger (LS). Celle-ci est résolue en utilisant des développements en série de Neumann. Plusieurs approches sont alors proposées pour calculer les différents termes de la série, en utilisant des approches analytiques ou numériques. Ainsi, dans les deux premiers chapitres, on établit une famille d'équation LS pour la polarisation dans le contexte de la conductivité thermique et de l'élasticité. L'opérateur de cette équation est optimisé afin d'obtenir la meilleur convergence de la série de Neumann et par conséquent la meilleur estimation des propriétés effectives du composite. L'estimation proposée est basée à la fois sur une série tronquée et une estimation du résidu de la série de Neumann. Le travail présenté au chapitre 3 concerne le calcul des propriétés de transport de masse en milieu poreux. De manière classique, la loi de filtration est donnée par la loi de Darcy à l'échelle macroscopique. Dans ce travail, on calcule les termes correctifs à l'équation de Darcy lorsque la condition de stricte séparation des échelles n'est pas vérifiée. Ces termes correctifs sont calculés numériquement en résolvant une équation LS et en utilisant un schéma itératif basé sur la transformée de Fourier Rapide (TFR). Finalement, au chapitre 4, on détermine numériquement des bornes pour les propriétés élastiques des polycristaux en utilisant toujours les approches basées sur la TFR. L'approche proposée permet de tenir compte de la géométrie exacte de la cellule de Voronoi en utilisant les expressions exactes des fonctions formes pour des polygones et des polyèdres. La méthode est appliquée à des polycristaux constitués de monocristaux cubiques.

Abstract

This work proposed some contributions to the homogenization methods with applications to composites materials, polycristals and porous media. The effective properties are determined by solving the unit cell problem and the corresponding Lippmann-Schwinger (LS) equation. The latter is solved by means of Neumann series. Different approaches are considered to evaluate each terms of the series using analytic or numerical approaches. In the first two chapters, we formulate a general class of LS equations for the polarization in the case of conductivity and then elasticity. The operator of the latter is optimized to obtain the best convergence of the associated Neumann series and then of the better estimate of the effective of the composite. The estimate is based on both a truncated Neumann series and an approximation of its residual. In chapter 3, we deal with the mass transport properties of porous media. Classically, the filtration law is given by the Darcy equation at the macroscopic scale. In the present work we compute the corrective terms of the Darcy equation in the situation of no strict scale separation. These corrective terms are determined numerically by solving a LS equation with a fast Fourier Transform (FFT) based iterative scheme. Finally, in chapter 4, we derivative numerically some bounds for the elastic properties of polycristalsstill by means of an FFT iterative scheme. The approach uses an exact description of the voronoi-unit cell geometry by using the shape functions of polygons and polyhedrons. The method is applied to polycristals constituted of cubic single crystals.