Viet-Thanh To

Vincent Monchiet
email: vincent.monchiet@u-pem.fr

Quy Dong

A FFT method for the Computation of thermal diffusivity of porous periodic media

Keywords: Porous media, Thermal diffusivity, Permeability, Fast Fourier Transform Viet-

In this paper, we provide Fast Fourier Transform iterative schemes to compute the thermal diffusivity of periodic porous medium. We consider the fluid flow through a porous rigid solid due to a prescribed macroscopic gradient of pressure and a macroscopic gradient of temperature. As already proved in the literature, the asymptotic homogenization procedure is reduced to the resolution of two separated problems for the unit cell: (i) the fluid flow governed by the Stokes equations with an applied gradient of pressure, (ii) the heat transfer by both convection and conduction due to an applied macroscopic gradient of temperature. We develop new numerical approaches based on fast Fourier transform for the implementation of the cell problems. In a first approach, a simple iterative based on the primal variable (gradient of temperature) is provided to solve the heat transfer problem. In order to improve the convergence in the range of high values of the prescribed gradient of pressure, we propose a more sophisticated iterative scheme based on the polarization. In order to evaluate their capacities, these FFT algorithms are applied to some specific microstructures of interest including flows past parallel pores (Poiseuille flows) and periodically or randomly distributed cylinders.

Introduction

Studying transport processes in porous media and determining their effective behavior are problems of fundamental interest. Using the volume average [START_REF] Carbonell | Heat and mass transfer in porous media[END_REF][START_REF] Hsu | Thermal dispersion in a porous medium[END_REF][START_REF] Moyne | Two-equation model for a diffusive process in porous media using the volume averaging method with an unsteady-state closure[END_REF][START_REF] Moyne | Thermal dispersion in porous media: one-equation model[END_REF][START_REF] Whitaker | Diffusion and dispersion in porous media[END_REF] or asymptotic homogenization [START_REF] Auriault | Upscaling heterogeneous media by asymptotic expansions[END_REF][START_REF] Auriault | Homogenization of coupled phenomena in heterogenous media[END_REF][START_REF] Bloch | Heat transfer in nonsaturated porous media. modelling by homogenisation[END_REF][START_REF] Sanchez-Palencia | Homogenization Techniques for Composite Media[END_REF], a large number of macroscopic models have been derived to capture successfully the phenomena. For periodic media, effective properties can generally be proved to exist uniquely and determined by solving the periodic boundary value problem on the unit cell. Regarding the coupled fluid flow/heat transfer in porous media, there are numerous works dedicated to the resolution of flow/thermal diffusion in porous media. Standard numerical methods based on Finite Element Method (FEM), Finite Volume Method (FVM) or Boundary element Methods (BEM) have been often considered for the problem of Stokes flow in porous media. For instance, [START_REF] Alcocer | Permeability of periodic porous media[END_REF][START_REF] Alcocer | Permeability of periodic arrays of cylinders for viscoelastic flows[END_REF][START_REF] Berdichevsky | Perform permeability predictions by self consistent method and finite element simulation[END_REF][START_REF] Edwards | The influence of reynolds number upon the apparent permeability of spatially periodic arrays of cylinders[END_REF][START_REF] Idris | Microstructural effects on the flow law of power law fluids through fibrous media[END_REF] determined the effective permeability for fluid flow past a regular array of cylinders. The computation of permeability for 3D microstructures has been performed by [START_REF] Barrere | Détermination numérique de la perméabilité en milieu poreux périodique tridimensionnel[END_REF] in the case of an array of sphere and more recently by [START_REF] Malinouskaya | Wave propagation through saturated porous media[END_REF] on a real media imaged by microtomography. The resolution with the FEM of the Stokes flow problems requires special care, since classical interpolation based on piecewise linear function fails to perform due to the flow incompressibility. In such situation, more sophisticated interpolation known as "MINI element" must be used [START_REF] Arbogast | A computatonal method for approximating a darcystokes system governing a vuggy porous medium[END_REF][START_REF] Arnold | A stable finite element for the stokes equations[END_REF][START_REF] Fortin | Old and new finite element for icompressible flows[END_REF]. Thermal dispersion in convective flow in porous media has been numerically investigated for 2D problems constituted of regular arrays of cylinders, squared or elliptic rods [START_REF] Alshare | Modeling of unsteady and steady fluid flow, heat transfer and dispersion in porous media using unit cell scale[END_REF][START_REF] Kuwahara | A numerical study of thermal dispersion in porous media[END_REF][START_REF] Moghari | A numerical study of non-equilibrium convective heat transfer in porous media[END_REF][START_REF] Pedras | Thermal dispersion in porous media as a function of the solid-fluid conductivity ratio[END_REF][START_REF] Vu | Forced convection of air through networks of square rods or cylinders embedded in microchannels[END_REF]. More recently, a new numerical method based on the fast Fourier Transform (FFT) has been proposed to numerically evaluate the effective permeability of porous media [START_REF] Ly | Computation of permeability with fast fourier transform from 3d digital images of microstructures[END_REF][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] Nguyen | A fourier based numerical method for computing the dynamic permeability of porous media[END_REF] on the basis of the method first introduced by Moulinec et al. [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] for elastic composites. The local cell problem for the fluid flow is solved by means of an iterative scheme which uses the periodic Green's tensor and the images of the microstructures (which come from imagerie technics for instance). The method has significantly increased performance since the memory needed for solving the problem is greatly reduced compared to other methods. In this paper we provide new FFT based iterative schemes to handle the problem of thermal diffusion in porous media. In a first approach, we propose an iterative formulated with the gradient of temperature. However, due to the presence of the convective term, the rate of convergence is greatly reduced at large values of the applied pressure gradient. This motivates an alternative iterative scheme based on the polarization. The paper is organized as follows. In section 2, we present the periodic boundary value problems and we recall the definitions for effective permeability and diffusivity. Section 3 is devoted to the implementation of cell problems with FFT algorithms. After a brief recall on the resolution of the stokes problem and the implementation of the local velocity (which appears in the convection term), we focus our attention on the heat transfer problem by both conduction and convection. In section 4, we provide some illustrations. First, we consider the simple case of flows past parallel pores which has the advantage to provide analytic expression of the conductivity which serves to evaluate the accuracy of our numerical solutions. Illustrations for 2D and 3D microstructures are also provided. Particularly, the resolution of the heat transfer problem past randomly oriented fibers is performed to show the capacity of the method.

Effective transport properties for coupled heat transfer and fluid flow problems

We consider a periodic porous medium saturated by a homogeneous Newtonian viscous fluid with the dynamic viscosity η f . Both the fluid and the skeleton are thermally conductive media with the diffusivity k f and k s respectively. By V , we denote the total volume of the cell, by V f and V s the volume occupied by the fluid and the solid respectively, where V s and V f 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:

                 η f ∆v -∇p = 0 in V f div v = 0 in V f v = 0 on Γ ∂T ∂t + v.∇T = ∇.(k∇T ) in V (1)
where v, p, T and k are respectively the velocity, pressure, temperature and thermal diffusivity. In our problem, the Reynolds number is assumed to be small Re ≪ 1 and the fluid is incompressible and stationary. In writing (1), the first three equations are valid for the fluid while the last equation is valid for both the solid and the fluid by admitting a null velocity field in the solid phase. The local diffusivity k is a piecewise constant function:

k = k s in solid phase, k = k f in fluid phase. (2)
The periodic homogenization procedure is based on the existence of the two following length scales: L the size of the porous medium and l the size of the periodic cell. The asymptotic expansion method in terms of the ratio ǫ = l/L ≪ 1 has been used to determine the equations corresponding to every orders of ǫ k . The detailed derivation of the effective permeability, diffusivity, etc... has been well documented in the literature [START_REF] Allaire | Homogenization of a convection-diffusion model with reaction in a porous medium[END_REF][START_REF] Auriault | Upscaling heterogeneous media by asymptotic expansions[END_REF][START_REF] Auriault | Homogenization of coupled phenomena in heterogenous media[END_REF][START_REF] Bloch | Heat transfer in nonsaturated porous media. modelling by homogenisation[END_REF].

The macroscopic description for the fluid flow in the porous medium is given by the Darcy law:

V = - 1 η f B.J (3)
where V is the macroscopic velocity, J the macroscopic pressure gradient and B the permeability tensor. The macroscopic equivalent description for the heat transfer is:

∂T ∂t + V .E = div (K.E) ( 4
)
where E is the macroscopic gradient of temperature and K is the effective diffusivity tensor. Finally, to find the homogenized permeability and diffusivity tensors of the porous medium, we must account for the periodicity conditions prescribed for the variables v, T and p when solving (1). These conditions will be given in the next. For the sake of clarity, two periodic boundary value problems will be treated separately, as described in the following:

-Stokes flow problem:

         η f ∆v -∇p -J = 0 in V f div v = 0 in V f v = 0 on ∂V f (5)
where v and p are V -periodic, i.e. they take the same values on two opposite sides of the unit cell. J is constant in the cell and denotes the applied macroscopic pressure gradient. Due to the linearity of the equations, the local velocity is on the form:

v = A.J (6)
where A is the two-order localization tensor. It components A ij are determined by solving Eq. ( 5) and by considering the component v i due to the applied component of the pressure gradient J j .

The permeability tensor B is given by:

B = -η f < A > V (7)
where < • > V denotes the volume average of the quantity "•" over the volume of the cell.

-Heat transfer problem:

∇. [k∇T ] = v.∇T -V .E (8)
where T -E.x, in which x denotes the vector position, is V -periodic and v is determined from the Stokes-flow problem. In other hand, the flux j.n, with j = k∇T and n is the outward normal unit vector taken on the boundary of the cell, is antiperiodic that is the consequence of the thermal flux equilibrium between the cell and its neighbors. Another approximations could be derived for the macroscopic filtration law considering the Forchheimer Law [START_REF] Forchheimer | Wasserbewegung durch boden[END_REF], Darcy-Brinkman Law [START_REF] Brinkman | A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles[END_REF][START_REF] Brinkman | On the permeability of media consisting of closely packed porous particles[END_REF] or Forchheimer-Brinkman law [START_REF] Hsu | Thermal dispersion in a porous medium[END_REF]. The determination of such macroscopic law in the framework of the asymptotic homogenization theory has been the subject of many studies. The derivation of the Forchheimer or more complex non linear filtration law has been proposed by Mei and Auriault [START_REF] Mei | The effect of weak inertia on flow through a porous medium[END_REF], Wodie and Levy [START_REF] Wodie | Correction non lineaire de la loi de darcy[END_REF], Giorgi [START_REF] Giorgi | Derivation of the forchheimer law via matched asymptotic expansions[END_REF], Skjetne and Auriault [START_REF] Auriault | New insights on steady, non-linear flow in porous media[END_REF] or Bahloff et al. [START_REF] Bahloff | Polynomial filtration law for low reynolds number flows through porous media[END_REF] by accounting for the inertial term which appears in the Navier-Stokes. Also the Brinkman equation can be derived from the asymptotic homogenization theory by accounting for the higher order terms of the series [1]. Note that the present works focus on the resolution of the Heat transfer problem and the algorithms developed in the present paper could be also used by considering a local velocity in the convection term which accounts for such effects.

Denoting by e the temperature gradient, its average over the volume is the macroscopic gradient of temperature E:

e = ∇T, E = e V . (9)
The solution of the heat transfer problem then depends on two forcing variables, the macroscopic gradient of temperature E and the gradient of pressure J due to the presence of the convection term (see relation [START_REF] Arbogast | A computatonal method for approximating a darcystokes system governing a vuggy porous medium[END_REF] giving the velocity as function of J).

Due to the linearity of the problem in E, there exists a uniquely zero mean function g satisfying to:

T = g.E + T , g V = 0 ( 10
)
where T is the average temperature. Note however that the local temperature T has a non linear dependence with the applied pressure gradient J.

The effective diffusivity K is given by:

K = k I + k∇ T g -(v -V ) ⊗ g ( 11
)
where I is the identity for two order tensors. Note that the above definition for the diffusivity tensor is non symmetric due to the presence of the dispersive term (v -V ) ⊗ g . However, the permeability tensor K is defined merely by its inner product with the symmetric tensor ∇E in relation ( 4), so that only the symmetric part K s of K has a contribution at the macroscopic scale. Note also that, still due to the presence of this convection term, the diffusion tensor K introduces a dependence with the macroscopic pressure gradient J.

The determination of the diffusivity then requires successively the resolution of the Stokes flow problem and then the heat transfer problem. Since the two problems are periodic, it is convenient to work with Fourier series and the associated FFT technique. Considering a parallelepipedic unit cell V , a V -periodic function φ admits the following representation:

φ(x) = n φ(ξ n ) exp(iξ n .x) (12)
where i 2 = -1 and φ(ξ) denotes the Fourier transform of φ(x) given by:

φ(ξ) = F (φ(x)) =< φ(x) exp(-iξ.x) > V ( 13
)
where ξ n are the wave vectors given by:

ξ n = 2πnζ, n = 0, ±1, ... ± ∞, ζ i = 1 h i (14)
and where h 1 , h 2 , h 3 are the dimensions of the cell. In Eq. ( 13) and in the next, we shall use the notation φn , a quantity symbol φ associated with wave vector ξ n , to refer to the discrete Fourier transform of φ(x) which is associated with the real variable x. When dealing with Fourier series representation, it is convenient to use the FFT based iterative schemes which use exact expressions of the Green operator in the Fourier space. These methods, initially introduced by Moulinec and Suquet [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] for composite materials, are particularly adapted when dealing with unit cell elementary problems with periodic conditions.

3 Resolution of cell problems with FFT

Resolution of Stokes problem

The description of the FFT method dealing with this problem have been presented in [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] Nguyen | A fourier based numerical method for computing the dynamic permeability of periodic porous media[END_REF] and will be recalled briefly for completeness. Generally, we can make use of the analogy between the Stokes equations and the elasticity problem of composites to formulate the integral equation in the sense of Lippmann-Schwinger-Dyson [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF][START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF][START_REF] Milton | The Theory of Composites[END_REF]. However, there are some details of viscous flow problem that render the problem different from the original approach developed for elastic composites. The unit cell of the porous medium can be assimilated to an elastic composite reinforced by perfectly rigid inclusions. For such class of problems a stress formulation is more suitable since it uses the compliance that is null within the inclusion while a strain based approach uses the rigidity that is singular in the inclusion. The compliance reads:

S(x) = 1 2η f I f (x)K (15)
in which the fourth order tensor K is defined by:

K = I -J, J = 1 3 I ⊗ I. (16)
In the above equation, I is the two order identity tensor, I is the fourth order identity tensor. J and K are fourth order projectors. The action of K on a two order tensor provides its deviatoric part while the action of J provides the spherical part. Moreover, I f (x) is the characteristic function of the volume V f given by:

I f (x) =    1 in V f 0 in V s (17)
The main difference with the problem of elastic composite lies in the forcing term, which could be represented here by the body force f which takes the value:

f = -J in V f , f = F + C ∧ x in V s . (18)
The body force in the solid phase represents the drag force due to the flow through the rigid skeleton. F is prescribed in order to makes null the translation while C is added to makes null the rotation of the solid. F and C physically represents the drag force due to the fluid motion. Introducing a reference medium of compliance S 0 , it has been shown (see [START_REF] Monchiet | A fft-based method to compute the permeability induced by a stokes slip flow through a porous medium[END_REF]) that the stress field is solution of the integral equation (written in Fourier space):

σ n = Ω n . f n -∆ n : S n * σ n -S 0 : σ n ( 19
)
where ∆ n = ∆(ξ n ) is the Green tensor for the stress whose expression can 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]. The symbol " * " represents the discrete convolution product defined by:

S n * σ n = m S(ξ n -ξ m ) : σ(ξ m ). ( 20
)
In Eq. ( 19), Ω n = Ω(ξ n ) where the components of the third order tensor Ω(ξ) are:

Ω ijk (ξ) = i |ξ| [2ξ i ξ j ξ k -δ ij ξ k -δ ik ξ j -δ jk ξ i ], ξ = ξ/ ξ (21)
if ξ = 0 but Ω ijk (ξ) = 0 if ξ = 0. In Eq. ( 21), ξ denotes the norm of ξ. The solution of the integral equation is computed with the following recurrence relation:

σ i+1 n = σ i n -∆ n : S n * σ i n (22)
which starts from the following first term:

σ 1 n = Ω n . f n . ( 23
)
At convergence of the iterative scheme, we compute the strain rate d n :

d n = S n * σ i n ( 24
)
and then the velocity field v n , obtained by the integration of the strain rate in the Fourier space:

v n = - 2i |ξ n | 2 d n .ξ n ∀ξ = 0. ( 25
)
The velocity field is defined by its Fourier coefficients for all values of ξ n except for n = 0. It means that the velocity field is defined up to an added constant which represents the mean velocity field V . The latter is identified by the condition v = 0 in the solid phase.

Resolution of heat transfer problem

Preliminaries

Consider a field j in real space which is in equilibrium with a source term q:

∇. j + q = 0 ( 26
)
and a curl free field e:

∇ ∧ e = 0. (27)
We denote by e n , j n and q n the Fourier coefficient of e, j and q respectively. To read the conditions ( 26) and [START_REF] Koutsourelakis | Stochastic upscaling in solid mechanics: an exercice in machine learning[END_REF] in Fourier space, it is convenient to introduce the projectors P n , Q n and the vector ω n defined, for n = 0, by:

P n = ξ n ⊗ ξ n ξ n 2 , Q n = I -P n , ω n = -i ξ n ξ n 2 (28)
and:

P 0 = Q 0 = 0, ω 0 = 0. ( 29
)
P n and Q n are two projectors introduced by Milton [START_REF] Milton | The Theory of Composites[END_REF]. Tensor ω n is also considered due to the presence of the source term in [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF]. In the Fourier space, equations ( 26) and ( 27) read:

P n . j n + ω n q n = 0, Q n . e n = 0. ( 30
)
When the last condition is verified, e derives from the temperature T , whose Fourier coefficients are obtained from:

T n = ω n . e n . (31)
In the next, two FFT based iterative schemes are proposed to solve the heat transfer problem.

A primal iterative scheme

In this section we provide a simple iterative scheme to compute the solution of the heat transfer problem. It is formulated with the primal variable: the gradient of temperature e (the associated dual variable being the flux j). Introducing in Eq. ( 8) a reference medium of diffusivity k 0 , it leads to:

∇.(k 0 ∇T + τ ) + q = 0 ( 32
)
where we have introduced:

τ = (k -k 0 )∇T, q = V .E -v.∇T. (33)
Eq. ( 32) can be also read in Fourier space (using Eq. ( 30)):

P n .(k 0 e n + τ n ) + ω n q n = 0, Q n . e n = 0. ( 34
)
The first relation expresses the equilibrium of the flux with the applied polarization τ and the source term q. The second relation is the compatibility for the gradient of temperature e.

Considering that τ n and q n are applied to the system, the solution for e n is:

e n = -G 0 n .( τ n + ω n q n ) for n = 0, e 0 = E (35)
where:

G 0 n = 1 k 0 P n ( 36
)
is the Green function of the reference medium of diffusivity k 0 . Replacing τ n and q n by their expressions as function of e n we obtain:

e n = -G 0 n .( k n * e n -k 0 e n -ω n ( v n * e n )) for n = 0, e 0 = E (37)
where the definition for the discrete convolution product " * " has been already given in Eq. [START_REF] Giorgi | Derivation of the forchheimer law via matched asymptotic expansions[END_REF]. The solution of this equation is expanded along Neumann series, each term of this series being computed from the following iterative scheme:

e i+1 n = -G 0 n .( k n * e i n -k 0 e i n -ω n ( v n * e i n )) for n = 0, e i+1 0 = E. ( 38
)
A simplification of the above iterative scheme is possible by accounting that for any curl free field e we have:

∀n = 0 G 0 n .(k 0 e i n ) = e i n (39)
and

G 0 n .(k 0 e i n ) = 0 (40)
if n = 0. Using this simplification, we obtain for the iterative scheme:

e i+1 n = e i n -G 0 n .( k n * e i n -ω n ( v n * e i n )) (41)
and which starts from the first term:

e 1 n = 0 for n = 0, e 1 0 = E. ( 42
)
At convergence, one can compute the temperature by:

T n = ω n . e n (43)
that leadings to g and then the effective diffusivity (see Eq. ( 10) and Eq. ( 11)). In Eq. ( 41), the convolution product between k n and e i n , also v n and e i n are effected by means of the FFT algorithm. This means that one inverse FFT and two FFT are needed at each iteration. One inverse FFT is needed to compute e i n in real space. Therefore, the product of e n with k n and v n is made in real space, in which the product is local. Thereafter two FFT are applied to compute ke and v.e in Fourier space and to apply the product with the Green tensor.

In [START_REF] Moyne | Two-equation model for a diffusive process in porous media using the volume averaging method with an unsteady-state closure[END_REF] the velocity field v is computed by solving the Stokes problem for which the stress based iterative scheme ( 22) is used together with [START_REF] Hsu | Thermal dispersion in a porous medium[END_REF] and [START_REF] Idris | Microstructural effects on the flow law of power law fluids through fibrous media[END_REF]. The velocity field being known in Fourier space from [START_REF] Idris | Microstructural effects on the flow law of power law fluids through fibrous media[END_REF], the inverse FFT is taken to obtain v in real space and is thereafter stored to make the product with e at each iteration in [START_REF] Moyne | Two-equation model for a diffusive process in porous media using the volume averaging method with an unsteady-state closure[END_REF]. The details about the numerical integration of the algorithm are given below:

Initialization: e 1 n = 0 for n = 0, e 1 0 = E Iteration:

e i n is known e i n = FFT -1 ( e i n ) j n = k n e i n q n = -v n .e i n j n = FFT(j n ) q n = FFT(q n ) e i+1 n = e i n -G 0 n .( j n + ω n q n ) (44)
The iterative scheme is stopped when the convergence criterion is satisfied. This criterion is based on the equilibrium of the flux:

n P n . j n + ω n q n 2 1/2 < ǫ ( 45
)
where ǫ is the precision chosen equal to 10 -3 in our calculations. This criterion represents a measure of the residue of the linear equation which have to be solved.

As shown in [START_REF] Michel | A computational scheme for linear and non-linear 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] for the case of elastic composite or [START_REF] Monchiet | A polarization-based fast numerical method for computing the effective conductivity of composites. accepted for publication in[END_REF] for the case of conductivity, basic schemes based on primal or dual variables encounter convergence rate problems at high contrast. Furthermore, with the presence of the convection term in the integral equations, the basic scheme may also diverge by increasing the value of pressure gradient. Indeed, the convergence is ensured if all the eigenvalues of the linear operator giving e i+1 n as function of e i n in Eq. ( 41) are comprised between -1 and 1. The recurrence relation (41) can be put into the form:

e i+1 n = m δ mn I -G 0 n .Z nm . e i m ( 46
)
where

Z nm = k(ξ n -ξ m )I -ω n ⊗ v(ξ n -ξ m ) (47)
and δ mn is the Kronecker symbol. Since e i n is compatible, it is invariant by the transformation e i n → P n . e i n for any n = 0. Moreover, the Green operator G 0 n is proportional to the projector P n with the constant 1/k 0 . It follows that Eq. ( 46) can be rewritten into the equivalent form:

       e i+1 n = P n . m δ nm I - 1 k 0 Z nm . e i m ∀n = 0 e i+1 0 = e i 0 (48)
A sufficient condition for having the convergence of the iterative scheme is that the operator

δ nm I - 1 k 0 Z nm (49)
has a norm less than 1. In the case of conductive phase (the convection term being zero), the convergence has been already examined in [START_REF] Milton | The Theory of Composites[END_REF]. It has been proved that the above condition is ensured for:

k 0 > 1 2 k max ( 50
)
where k max is the higher value taken by the diffusivity (here the higher value between the diffusivity of the fluid and of the solid). The following value for k 0 can be chosen in order to obtain the better rate of convergence:

k 0 = 1 2 (k f + k s ) (51)
where k f and k s denote the diffusivity of the fluid and of the solid respectively. The rate of convergence decreases by increasing the contrast (defined as the ratio between the diffusivity of the phases) and for infinity values of the contrast the iterative scheme diverges. For the dispersive heat transfer problem, the derivation of convergence conditions is more difficult due to the presence of the velocity which depends on the microstructure geometry. Obviously, the eigenvalues of the operator Z nm could be evaluated numerically. However the computation of Z nm is very long and memory costly (particularly for high dimension problems such as 3D microstructures). In [START_REF] Nguyen | A fourier based numerical method for computing the dynamic permeability of periodic porous media[END_REF], this operator is not computed since the product between Z nm and e i m is effected in the real space in which this product becomes local. Despite this incapacity to obtain simple convergence conditions, the optimal reference modulus k 0 can be determined numerically on a lower grid. Note however that the velocity field, which enters in the recurrence relation, linearly depends on the pressure gradient. This suggests that the case of high values of the gradient of pressure is comparable to that of high values of the contrast for which the iterative scheme leads to large number of iterations. To improve the convergence issue, we propose in the next section a more sophisticated iterative scheme based on the polarization.

A polarization based iterative scheme

A polarization based iterative scheme has been provided in [START_REF] Monchiet | A polarization-based fft iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF][START_REF] Monchiet | A polarization-based fast numerical method for computing the effective conductivity of composites. accepted for publication in[END_REF] in the case of elastic composites. It improves the convergence of simple iterative schemes based on the strain or stress in the domain of high contrasts. In this section we propose to adapt this approach for the heat transfer problem. To this end, let us first introduce the polarization:

τ = (k -k 0 )e. ( 52
)
It follows that the gradient of temperature and the flux are given by:

e = (k -k 0 ) -1 τ , j = τ + k 0 e (53)
The equilibrium of the flux with the convection term can be read:

G 0 n . j n + ω n q n = 0 (54)
while the compatibility for e can be also read:

D 0 n . e n = 0 (55)
where D 0 is the dual Green operator whose coefficients are given, for n = 0, by:

D 0 n = k 0 Q n ( 56
)
where Q n is defined in [START_REF] Kuwahara | A numerical study of thermal dispersion in porous media[END_REF] and D 0 n = 0 when n = 0. Additionally, the average of gradient of pressure must be equal to E, which reads in Fourier space:

e 0 = E (57)
To summarize, we search the value of the polarization τ which comply to relations (53) with the conditions (54), ( 55) and (57). The polarization must then be computed with the following iterative process:

Initialization: τ 1 n = 0 Iteration: τ i n is known τ i n = FFT -1 ( τ i n ) e i n = (k -k 0 ) -1 τ i n j i n = τ i n + k 0 e i n q i n = -v n .e i n e i n = FFT(e i n ) j i n = FFT(j i n ) q i n = FFT(q i n ) τ i+1 n =      τ i n -αK 0 . G 0 n .( j i n + ω i n q i n ) +α D 0 n . e i n
, for n = 0

τ i+1 0 = τ i 0 + αK 0 .( e i 0 -E), for n = 0 (58)
When the convergence is reached, that corresponding to τ i+1 n = τ i n , the conditions (54), (55) and (57) are recovered for the flux and gradient of temperature. This new iterative scheme uses one inverse FFT and three FFT at each iteration that is one FFT more than the simple one based on the primal variable. However, it is expected to reduce significantly the number of iterations at convergence with this new algorithm. The capacity of these methods of resolutions is evaluated in the next section. Let us now discuss about the choice of k 0 and coefficients α. When the convection term is null (that corresponding to J = 0), the convergence of the algorithm is guaranteed if (see [START_REF] Monchiet | A polarization-based fast numerical method for computing the effective conductivity of composites. accepted for publication in[END_REF]):

α ∈]0, 2[, k 0 ∈] -∞, 0[ (59)
Considering now the convection term, the analysis of convergence encounters the same difficulties that for the simple iterative scheme and it is some not guaranteed that the conditions (59) are still valid when the convection term is added. Consequently, the optimal choice of k 0 must then be evaluated numerically.

For the accelerated scheme, three conditions are tested to stop the iterations: (i) a first for the equilibrium of the flux, (ii) a second for the compatibility of the gradient of temperature and (iii) a condition is written for the equality between the applied macroscopic gradient of temperature and the average of the local gradient of temperature.

(i)

n P n . j n + ω n q n 2 1/2 < ǫ (ii) n Q n . e n 2 1/2 < ǫ (iii) e 0 -E < ǫ (60)
Note that the conditions (ii) and (iii) are verified at any step of the primal iterative scheme. For the polarization iterative scheme the iteration is stopped when the three conditions (i), (ii) and (iii) are satisfied.

Applications

The flow between two parallel planes

We consider a benchmark problem where the porous media is constituted of parallel pores of height h regularly distributed along the direction Ox 2 . The volume fraction of the pores is denoted f . The fluid flow is generated by applying the pressure gradient J 1 along the direction Ox 1 . An analytic solution is given in appendix A and is used to evaluate the accuracy of the FFT numerical solutions. For the numerical implementation with FFT, a squared unit cell of dimension 2H = 1 along each space direction is considered and a rectangular pore is centered at the origin (see Fig. 1). The volume fraction is of the pores is f = h/H. For the applications, we put f = 0.5 and the diffusivity contrast is k f /k s = 10. First of all, we analyse the convergence of the two FFT based iterative schemes: the basic scheme formulated with the gradient of temperature (cf. ( 44)) and the accelerated scheme based on the polarization (cf. ( 58)). As already mentioned in the last section, the convergence of these iterative schemes depends on the contrast between the phase diffusivity but also depends on the value of the applied pressure gradient. Unfortunately, it has not been possible to derive the conditions giving the convergence for the reference diffusivity k 0 , the latter must be determined numerically. Since the optimal value of k 0 is nearly independent on the number of wave vectors, it can be determined on a lower size grid. In the present example, we use k 0 = 12 for the basic scheme and k 0 = -3 for the accelerated scheme. Also, the value α = 3/2 considered in [START_REF] Monchiet | A polarization-based fast numerical method for computing the effective conductivity of composites. accepted for publication in[END_REF] has been taken. With these choices, it possible to solve the thermal diffusion problem on a large range of pressure gradients. In this section and in the next parts of the paper, we work with the dimensionless pressure gradient J * defined by:

J * = (2H) 3 η f k s J (61)
On Fig. 2, we represent the number of iterations as function of the applied gradient of pressure J * 1 for the basic and the accelerated scheme. Although for low values of J * 1 the number of iterations at convergence is lower for the basic scheme, it can be observed that the accelerated scheme has a better rate of convergence for higher values of J * 1 . Also, it must be noted that the number of iterations with the basic scheme linearly increases with the applied gradient of pressure J * 1 (in the log frame). However, with the accelerated scheme, the number of iterations is quasi independent of the gradient of pressure. Since it has a better rate of convergence, only the accelerated scheme is used for the applications provided in the next sections of the paper. On Fig. ?? we represent the effective diffusivity K 11 (re-scaled by k s ) as function of the gradient of pressure J * 1 obtained with the basic and accelerated scheme. The results are compared with the exact solution. Note that for the numerical implementation, N = 64 wave vectors are considered along each space direction. Note also that the effective diffusivity K 11 displays two regimes :

at low values of the pressure gradient, the effective diffusivity is not influenced by the pressure gradient and is equal to that computed at J * 1 = 0. -in the range of higher values of pressure gradient (J * 1 > 10 3 ), the convective term becomes predominant in the heat transfer problem. The dependence of the effective diffusivity with the applied gradient of pressure J * 1 is linear in the log-log frame. This suggests that only two points must be computed to determine the asymptotic behavior in the domain of high pressure gradient. Since the three curves are very closed on Fig. (3), a zoom is provided to display their differences. A good agreement between the exact and the numerical solution is observed. The differences observed with the exact solution can be reduced by increasing the number of wave vectors. For completeness, we then provide, on tables 1 and 2, the values of effective diffusivity K 11 /k s computed with the applied gradient of pressure J * 1 = 1000 and J * 1 = 10000 and for various values of the number of wave vectors. Table 1 displays the results obtained with the basic scheme while table 2 provides the same results computed with the accelerated scheme. In each table, we also provide (in the last line) the exact value of the effective diffusivity. The errors in percents, are given between the parenthesis. It can be noted that the numerical solutions converge uniformly to the exact solution by increasing the number of wave vectors. It is observed that for the lower values of the applied pressure gradient, the primal iterative scheme gives a better estimation of the effective permeability (see table 1 and 2 for J * 1 = 10 3 ) whereas at higher values of the applied pressure gradient (see 1 and 2 for J * 1 = 10 4 ) the accelerated scheme provides a more accurate solutions. In the case of elasticity, a hierarchization of the FFT solutions has been proved in [?]. It has been show that the iterative schemes based on the primal or dual variables always provide a more accurate solution than that obtained with an accelerated scheme. This is also true for the heat transfer problem without the convection term. When this term is accounted in the computation of the macroscopic diffusivity, the hierarchization of the FFT solutions is uncertain and depends on the velue of the applied pressure gradient. Table 2 Effective diffusivity K 11 /ks computed with the accelerated scheme as function of the gradient of pressure J * 1 and the number of wave vectors. Solutions obtained with the accelerated FFT scheme.

J * 1 = 1000 J * 1 =

Flow past an array of cylinders

We now consider the flow through a periodic network of cylinders with circular cross section aligned in the direction Ox 3 . The radius of the cylinders is denoted R. For the numerical implementation of the FFT schemes, a squared unit cell of dimension 2H = 1 and containing a circle of radius R located at the origin has been considered (see figure 4). The problem is two dimensional and N = 256 wave vectors are considered in each space directions. Also, only the accelerated scheme has been considered for this application.

2R 2H 2H

Fig. 4 Unit cell of the periodic porous solid with aligned cylinders of circular cross sections.

On Fig. 5 we represent the variations of the effective diffusivity K 11 as function of the ratio R/H. In this illustration, the flow is generated by imposing the gradient of pressure J * 1 . The contrast in the local diffusivity of the solid cylinder and the fluid is k f /k s = 10 where k s and k f referred to the solid and the fluid respectively. The computation are performed by considering various values of the gradient of pressure J * 1 . It can be observed that, for any value of the applied pressure gradient, the effective diffusivity decreases by increasing the radius of the cylinder, since its diffusivity is lower that of the fluid. Also, the gradient of pressure increases the effective diffusivity of the porous medium. Now, we investigate the dependence of the effective diffusivity K 11 with the direction of flow. Indeed, in the last example, the diffusivity K 11 is computed by imposing the macroscopic gradient of temperature E 1 and a gradient of pressure in the same direction. Now, the macroscopic gradient of temperature E 1 is still applied to compute K 11 but the flow is generated along different directions in the plane Ox 1 x 2 . To this end, we put for the component of the gradient of pressure:

J * 1 = J * cos(θ), J * 2 = J * sin(θ) (62)
where J * is the norm of the applied pressure gradient J * and θ lies in the range [0, 2π]. From Fig. 6 to Fig. 8 we provide the variations of the effective diffusivity with the direction of flow in a polar frame. On each figures various values of the norm of the pressure gradient J * have been considered. It must be observed that the diffusivity is affected by the flow direction, particularly when the flow is parallel to the applied gradient of temperature. When the gradient of pressure is applied in the direction transverse to the temperature gradient (θ = 90 • ) the effective diffusivity is not affected by J * . It must be also noted that the influence of the flow direction is more important when the radius of the cylinder is small. For larger values of the cylinder (for instance, R/H = 0.8 on Fig. 8) these effects are greatly reduced.

Porous cell with randomly distributed and oriented fibers

In this section, we apply the method to 3D microstructures made up of interconnected random fibers. One example of such microstructure is provided on Fig. 9. To investigate the dispersion effect, the effective diffusivity is computed for various microstructure realization. The total number of fibers is N = 10 with the radius R/H = 0.05 (where 2H is the dimension of the cubic unit cell).

In order to obtain an interconnected porous medium, each cylinder intersects another by their respective axes but the orientation and position of the cylinders are chosen arbitrarily. The total volume of the solid is then uncontrolled and this also implies some dispersion in the porosity value. The calculations are performed by taking 256 wave vectors along each space directions. For the Stokes flow problem, a stress based iterative scheme is used, this involving 6 × 256 × 256 × 256 ∼ 100.10 6 DOFs. For the heat transfer problem, the method of resolution uses the polarization that is a vector, the total number of DOFs is 3 × 256 × 256 × 256 ∼ 50.10 6 . A contrast k f /k s = 10 is taken as in previous examples. A gradient of pressure J * 3 is prescribed in order to generate the flow in direction Ox 3 . The effective diffusivity K 33 is computed by imposing the gradient of temperature E 3 . The results are presented on table 3 and on Figs. 10 and11. The values of effective diffusivity K 33 (normalized by k s ) are provided on table 3 for three values of the pressure gradient. In the last column, are provided the values of the porosity. The values computed for the effective diffusivity are also given on Figs. 10 and 11 on which is also provided the cumulated average value of K 33 /k s . It is observed that the dispersion in the results are greatly increased with the value of the gradient of pressure. This effect will be attributed to both the low value of porosity and the fibers distribution influence. Indeed at higher values of gradient of pressure it has been already found in previous examples that the convection has a prominent influence on effective diffusivity. The dispersions in the results are then greatly increased between the different occurrences. A finely analysis of these uncertainty and on the choice of the representative volume element would be performed following [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF][START_REF] Koutsourelakis | Stochastic upscaling in solid mechanics: an exercice in machine learning[END_REF]. This is not the main goal of the present application.

Conclusion

This paper provides a numerical method to compute the effective diffusivity of porous material saturated by a viscous fluid. The diffusivity is evaluated in the framework of periodic homogenization and consists in solving two unit cell problems successively: the Stokes flow problem under an applied pressure gradient, the heat transfer problem due to both conduction and convection. The resulting effective diffusivity then depends on the applied gradient of pressure. Two methods of resolution based on FFT and Green operators are provided to solve the heat transfer problem. A first simple one, is based on the gradient of temperature and uses the associated primal Green tensor. An accelerated scheme based on the polarisation and using both the primal and dual Green tensors is also proposed to increase the rate of convergence. The accuracy and capacities of these two iterative schemes are evaluated through various 2D and 3D examples. First, the case of flows in parallel pores is considered as benchmark problem to compare the basic and the accelerated scheme and to check their accuracy with the exact solution. As expected, the accelerated scheme has been found to converge more rapidly and is only considered for the other applications. As a second application, the effective diffusivity is computed in the case of a regular array of cylinders. Particular attention is made on the role of the flow direction on the homogenized diffusivity. In a last example, we solve the problem for a complex 3D microstructure made up of interconnected cylinders of random orientations and which can mimic more realistic some fibers made porous materials.

A Poiseuille solution

In the Poiseuille problem where the pressure gradient acting along the flow direction J = Je 1 , the local velocity field of the fluid admits following parabolic form

v 1 (x 2 ) = J 2η f x 2 2 -h 2 |x 2 | ≤ h, v 1 (x 2 ) = 0 h ≤ |x 2 | ≤ H (63)
In equation ( 63), h is the pore height and H the height of the unit cell (see figure 1). The macroscopic velocity field is:

V 1 = - Jh 3 3Hη f ( 64
)
For the heat transfer problem, a macroscopic temperature gradient is prescribed in the same direction of the flow, i.e. E = Ee 1 . Adopting the following general form of temperature field in the fluid and solid phase

T = Ex 1 + ψ(x 2 ), (65)
the differential equations for ψ(x 2 ) in both the solid and the fluid phase are:

Fluid phase:

k f ψ ′′ (x 2 ) = v 1 (x 2 ) + ks J * h 3 3η f H E = J * ksE H 3 x 2 2 2 - h 2 2 + h 3 3H
Solid phase: ksψ ′′ (x 2 ) = J * ksEh 3 3H 4

(66)

in which we have used the dimensionless pressure gradient:

J * = JH 3 ksη f (67)
After integration, we obtain (accounting for the symmetry of the problem):

Fluid phase: ψ(x 2 ) = J * ksE H 3 k f x 2 2 h 2 4 - x 4 2 24 - h 3 6H x 2 2 + C 1 Solid phase: ψ(x 2 ) = - J * ksE H 3 ks h 3 6H x 2 2 + C 2 + C 3 |x 2 | (68)
The integration introduces three constants C 1 , C 2 and C 3 which have to be determined. The continuity of the temperature and the flux at the interface between the solid and the fluid (at x 2 ± h) leads to:

C 3 = - J * h 3 ksE 3H 3 ks , C 1 = C 2 - J * h 4 ksE 24H 4 4h -5H k f + 4(2H -h) ks (69)
The condition < T > V = 0 yields to:

C 2 = J * h 3 ksE H 5 (5h -6H)h 2 45k f + 2H 3 + 3Hh 2 -2h 3 18ks (70)
Note however that C 2 does not enter into the calculation of the homogenized permeability since it appears in (68) as a constant temperature. The calculation of the permeability from relation [START_REF] Barrere | Détermination numérique de la perméabilité en milieu poreux périodique tridimensionnel[END_REF] leads to:

K 11 = (1 -f )ks + f k f + (J * ) 2 k 2 s f 6 f (51 -84f + 35f 2 ) 945k f + (1 -f ) 3 27ks (71)
Here only the component of K 11 is provided since a macroscopic gradient of pressure along the direction Ox 1 has been prescribed to the unit cell. The component K 22 is not affected by the fluid flow and remains unchanged from the classic solution:

K 22 = 1 -f ks + f k f -1 (72)

2 Fig. 1

 21 Fig.1Unit cell of the periodic porous solid with parallel pores.

Fig. 2

 2 Fig. 2 Number of iterations at convergence of the iterative schemes as function of the applied (dimensionless) pressure gradient J * 1 . The volume fraction is f = 0.5 and the contrast is k f /ks = 10.

Fig. 3

 3 Fig.3Effective diffusivity of the porous medium with parallel pores. The volume fraction is f = 0.5 and the contrast is k f /ks = 10. Comparison of the FFT solutions with the exact one.

Fig. 5

 5 Fig. 5 Effective permeability (re-scaled by ks) as function of the radius of the cylinder for a contrast k f /ks = 10 and for various values of the dimensionless pressure gradient J * 1 .

Fig. 6 Fig. 7 Fig. 8

 678 Fig. 6 Variations of the effective diffusivity with flow direction for R/H = 0.2

Fig. 9

 9 Fig. 9 Exemple of a unit cell of the porous medium with randomly distributed and oriented fibers

Fig. 10 Fig. 11

 1011 Fig. 10 Effective diffusivity for random fibers and for different occurrences (line with circles) and cumulated value of effective diffusivity (line with diamonds) for J * 3 = 10. The number of fibers is N = 10 and the contrast is k f /ks = 10

Table 1

 1 Effective diffusivity K 11 /ks computed with the basic scheme as function of the gradient of pressure J * 1 and the number of wave vectors. Solutions obtained with the simple FFT scheme.

	10000

Table 3

 3 Effective diffusivity for random fibers and for different occurrences. The number of fibers is N = 10 and the contrast is k f /ks = 10

	N 0			K 33 /ks			Volume fraction
			J * 3 = 10	J * 3 = 1000	J * 3 = 2000
	1		1.1257	4.2988		12.5952		0.0311
	2		1.1087	4.1956		12.4108		0.0308
	3		1.1291	5.3233		17.2091		0.0314
	4		1.1305	2.6669		6.7024		0.0321
	5		1.1200	3.5004		9.6789		0.0316
	6		1.1114	4.4589		13.3691		0.0306
	7		1.1297	5.9579		19.3942		0.0314
	8		1.1138	3.0092		8.1635		0.0293
	9		1.1294	2.9117		8.0133		0.0319
	10		1.1442	2.4569		5.8313		0.0323
	K 33 /k s
	10
	9
	8					Value at each test Cumulated average value
	7
	6
	5
	4
	3
	2
	1
	0
	1	2	3	4	5	6	7	8	9	10
				Number of occurences


Hai-Bang Ly
email: ly@univ-mlv.fr

Vincent Monchiet
email: vincent.monchiet@univ-paris-est.fr

Daniel Grande
email: grande@icmpe.cnrs.fr

Computation of permeability with Fast Fourier Transform from 3-D digital images of porous microstructures

Keywords: Porous Media, Permeability, Homogenization, Fast Fourier Transform, Viscous Fluid

come

Computation of permeability with Fast Fourier Transform from 3-D digital images of porous microstructures

Hai Bang Ly, Vincent Monchiet, Daniel

Introduction

The determination of permeability of porous media in connection with microstructure morphology has been addressed with homogenization techniques based on asymptotic series expansion methods [START_REF] Auriault | Study of macroscopic behavior of a deformable porous medium[END_REF][START_REF] Sanchez-Palencia | Nonhomogeneous Media and Vibration Theory[END_REF][START_REF] Levy | Fluid flow through an array of fixed particles[END_REF] or on energy principle and volume averaging [START_REF] Whitaker | Flow in porous media I: A theoretical derivation of Darcy's law[END_REF][START_REF] Allaire | Homogenization of the Stokes flow in a connected porous medium[END_REF]. Analytic methods have been first developed in the literature for solving the associated unit cell problem. For instance, some works use expansions along eigenfunctions [START_REF] Sangani | Slow flow past periodic arrays of cylinders with application to heat transfer[END_REF][START_REF] Sangani | Slow flow through a periodic array of spheres[END_REF][START_REF] Wang | Stokes flow through an array of rectangular fibers[END_REF][START_REF] Wang | Stokes flow through a rectangular array of circular cylinders[END_REF], however such approaches are limited to some simple geometries corresponding to the flow through regular arrays of cylinders or spheres. The use of high-performance computers open number of possibilities for considering more realistic geometries, based on reconstructed porous cell or directly based on digital images obtained from X-ray microtomography. This also encourages the development of robust and efficient numerical tools for computing the effective permeability. For instance, standard numerical methods based on Finite Elements (FEM) or Finite Volumes have been often considered for computing the permeability of porous media [START_REF] Barrere | Détermination numérique de la perméabilité en milieu poreux périodique tridimensionnel[END_REF][START_REF] Cioranescu | Homogenization of the Stokes problem with non-homogeneous slip boundary conditions[END_REF][START_REF] Alcocer | Permeability of periodic porous media[END_REF][START_REF] Alcocer | Permeability of periodic arrays of cylinders for viscoelastic flows[END_REF][START_REF] Sawicki | Numerical investigation of the fluid flow through rotating porous media at both the microscopic scale and macroscopic scale[END_REF][START_REF] Bernard | Permeability assessment by 3D interdendritic flow simulations on microtomography mappings of Al ŰCu alloys[END_REF][START_REF] Stylianopoulos | Permeability calculation in three dimensional isotropic and oriented fiber network[END_REF][START_REF] Malinouskaya | Wave propagation through saturated porous media[END_REF][START_REF] Tahir | Influence of fiber orientation on the transverse permeability of fibrous media[END_REF]. The computation of permeability by means of the Boundary Element Method (BEM) has has been also proposed [START_REF] Borne | Harmonic Stokes flow through periodic porous media: a 3D Boundary Element Method[END_REF] or using the Lattice Boltzmann method [START_REF] Koivu | Transport properties of heterogeneous materials. Combining computerised Xray micro-tomography and direct numerical simulations[END_REF][START_REF] Spaid | Lattice Boltzmann methods for modeling microscale flow in fibrous porous media[END_REF][START_REF] Koponen | Permeability of Three-Dimensional Random Fiber Webs[END_REF]. When the microstructure is defined by voxels, the FEM can be used by considering a regular mesh with cubic elements. However the use of FFT type algorithms is more adapted for reducing the memory saving. Indeed, the method do not need the storage of stiffness or interaction matrices but only the storage of the nodal variables and the Fourier coefficients of the Green tensors. Moreover, the case of incompressible constituents (the fluid) does not requires specific treatments while, the FEM uses more sophisticated interpolation (for instance the "MINI" elements [START_REF] Fortin | Old and new finite element for icompressible flows[END_REF][START_REF] Arnold | A stable finite element for the Stokes equations[END_REF][START_REF] Arbogast | A computatonal method for approximating a Darcy-Stokes system governing a vuggy porous medium[END_REF]) which generally increases the size of the system. In this paper, we apply the method based on Fast Fourier Transform for computing the permeability of complex 3-D microstructures defined by digital images [START_REF] Monchiet | A FFT-based method to compute the permeability induced by a Stokes slip flow through a porous medium[END_REF]. Various problems and porous structures with increasing complexity are considered. First, the case of a flow around a cylinder and a sphere is considered as benchmark problems for which comparison with FEM or existing literature data are provided. Then, we construct more realistic porous media defined as assemblies of overlapping spherical voids to evaluate the capacity of this method.

The FFT method in short

We consider a unit cell containing an interconnected rigid solid Ω s and an interconnected porous open set Ω f saturated by a fluid. The flow is described by the Stokes equations and is generated by an applied pressure gradient denoted G. The unit cell problem reads:

µ∆v -∇p -G = 0 ∈ Ω f div v = 0 ∈ Ω f v = 0 ∈ S (1)
with the periodicity conditions for the velocity and the pressure. Let us denote by σ, the stress defined by:

σ = µ(∇v + ∇ T v) -pI (2)
where I is the two order identity tensor.

The cornerstone of the FFT numerical method is the reformulation of the system of partial differential equation 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 by making a continuation by continuity of the local fields within the porous medium. Particularly, we consider a null velocity field within the solid phase. Alternatively, a null velocity field can be also recovered by taking the Stokes equation within Ω s and a dynamic viscosity µ s that is large enough behind that of the fluid. This viscosity can be interpreted as a penalty coefficient in order to recover the condition v = 0 within Ω s . By doing so, the Stokes problem can be recast in an integral form [START_REF] Monchiet | A FFT-based method to compute the permeability induced by a Stokes slip flow through a porous medium[END_REF]:

σ = Ω. f -∆ 0 : d(σ) (3)
where the notation with the "hat" indicates that we use the Fourier transform of the corresponding quantity. In (3), d(σ) represents the Fourier transform of the strain rate tensor computed with the stress σ. The strain rate being defined by:

d(σ) = A(x) * [σ -pI] (4)
with :

A(x) = 1 2µ I f (x) + 1 2µ s I s (x) (5)
where I f (x) and I s (x) are the characteristic functions of the fluid and the solid phase. I f (x) is equal to 1 if the point x belongs to the pore space and is equal to 1 otherwise while I s (x) = 1 -I f (x).

In equation ( 3), ∆ 0 is the complementary Green operator associated with the reference medium of dynamic viscosity µ 0 . Both ∆ 0 and Ω are explicitly known in the Fourier space. The Green tensor ∆ 0 is defined by:

∆ 0 = 2µ 0 k ⊥ ⊗ k ⊥ + k ⊥ ⊗k ⊥ (6)
for n = 0 and ∆ 0 n = 0 for n = 0 and where k and k ⊥ are given by:

k = 1 |ξ| 2 ξ ⊗ ξ, k ⊥ = I -k (7)
The components of the third order tensor Ω are:

Ω ijk (ξ) = i |ξ| [2ξ i ξ j ξ k -δ ij ξ k -δ ik ξ j -δ jk ξ i ], ξ = ξ/ |ξ| (8)
In equation ( 3), f is given by:

f =        G in Ω f - c f cs G in Ω s (9)
where c f and c s denote the volume fraction of the fluid and the solid phase respectively. f can be interpreted as a source term in the Stokes equation which is equal to the the applied pressure gradient in the fluid phase. A constant value for f is considered in the solid phase and can be interpreted as the drag force due to the fluid flow around the solid Ω s . This drag force is computed in order to obtain the equilibrium of the force applied to the unit cell (see [START_REF] Monchiet | A FFT-based method to compute the permeability induced by a Stokes slip flow through a porous medium[END_REF] for more details).

The solution of the integral equation ( 3) is solved by means of the following iterative process:

σ i+1 = σ i -∆ 0 : d(σ i ) (10)
which can be interpreted as a fixed point iterative method associated with the linear system (3). The iterative scheme (10) is initialized with :

σ i=1 = Ω. f (11)
Once the convergence is achieved, one can compute the velocity field from the strain rate tensor d:

v = d : Ω (12)
for all ξ = 0. The velocity field is defined by its Fourier coefficients for all values of ξ except for ξ = 0. It means that the velocity field is defined up to an added constant that represents its mean value over the volume of the cell. This constant is identified by considering the condition v = 0 in the solid phase. Let us denote by v * the velocity field computed from equation [START_REF] Alcocer | Permeability of periodic porous media[END_REF] for any ξ = 0 and v * = 0 for ξ = 0. This velocity field has a null volume average over Ω and it is constant within Ω s . The total velocity field v is then the sum of v * and a constant velocity V which represents the average of v over the volume Ω:

v = V + v * ( 13
)
Since v is null within Ω s , the macroscopic velocity field is given by V = -v * in Ω s . From a practical point of view, the condition v(x) = 0 is not rigorously verified everywhere in Ω s . The macroscopic velocity is then computed by considering one point taken within Ω s or by taking the average of v * over the volume Ω s . When one point in Ω s is considered, this point must not be taken so closed to the surface of the solid. Indeed, at the vicinity of the surface, the occurrence of the Gibbs phenomenon, that is inherent with the use of Fourier series, could affect the accuracy of the solution.

For the implementation of the iterative scheme we put 1/µ s = 0. Obviously a very high value is considered for the penalty coefficient, however the numerical application shows that with the value 1/µ s = 0 a good rate of convergence is observed. Moreover, we choose µ 0 = 2µ that is optimal for the convergence of the algorithm.

Discretization of the problem

Following the approach developed for composite materials [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], the problem is discretized along a regular grid consisting of N × N × N voxels. The nodal variables are the components of the stress that depends on 6 independent components. The constituents (the solid and the fluid) distribution within the unit cell is accounted in the iterative process by the characteristic function I f which appears in relation [START_REF] Whitaker | Flow in porous media I: A theoretical derivation of Darcy's law[END_REF]. At each node of the regular grid, the value 0 or 1 is assigned depending on either this point is in the solid or it is in the fluid phase. After discretization, the convolution product between I f and the stress is made by means of the FFT algorithm. In equation ( 10), all the components of the Green tensor are represented by their Fourier coefficients in the real space in terms of the discrete wave vectors ξ n given by

ξ n = 2πnζ, n = -N..N -1, ζ i = 1 h i (14)
and h 1 , h 2 , h 3 are the dimensions of the cell along each directions of the cartesian frame. Obviously, the problem is discretized along each space direction, this would involve the use of two indices n 1 and n 2 for 2-D problems and three for 3-D ones. However, only one indice n is used for simplicity. The product between the components of the strain rate tensor and those of the Green operator just consists in the Hadamard product of tables of dimension N × N × N.

The discretized algorithm is:

1 : compute I f n = I f (x n ) 2 :
compute Ω n = Ω(ξ n ) following Eq. ( 8)

3 : compute ∆ 0 n = ∆ 0 (ξ n ) following Eq. (6) 4 : compute f n = -1 cs GI f n 5 : f n = FFT(f n ) 6 : initialize σ 1 n = Ω n . f n 7 :
while test < ǫ do 8 :

σ i n = FFT -1 ( σ i n ) 9 : p i n = tr(σ i n )/3 10 : d n = 1 µ I f n (σ i n -p i n I) 11 : d i n = FFT(d i n ) 12 : v i n = d i n : Ω n 13 : test = n v i n -v i-1 n 2 14 : σ i+1 n = σ i n -∆ 0 n : d n 15 : end while (15)
For simplicity, we have used the following notations for any variable F (x) and its Fourier transform F (ξ):

F n = F (x n ) and F n = F (ξ n ).
At step 7, we use a convergence test based on the norm of the difference between the local velocity computed at the iteration i and i -1 and we choose ǫ = 10 -4 . It is noteworthy that another test has been considered in previous works to check the convergence of the FFT algorithms. For instance, a convergence test based on the compatibility for the strain has been used [START_REF] Monchiet | A polarization based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF]. However, since comparisons are provided with the solution obtained with the FEM, we have chosen a convergence test applicable to both methods. From a practical point of view, the permeability is computed by taking µ = 1 and we apply the pressure gradient G = -e i for i = 1, 2, 3 where (e 1 , e 2 , e 3 ) is the unit vector of the cartesian basis. With this choice, the average of the velocity field, V i =< v i > Ω is exactly V i = K ij when G = -e j is prescribed to the system. In any case considered in the next of the paper, the total number of nodes is increased and the permeability is computed for different values of N in order to evaluate the accuracy of the solution. Considering the maximum possible memory occupancy of our computer cluster, the resolution cannot exceed 1024 × 1024 × 1024 voxels when using the FFT method. Calculations were performed using MATLAB 2011 on an AMD Opteron processor, 2.5GHz, 1024GB RAM and 64 cores.

Benchmark problems

Stokes flow through arrays of cylinders and spheres has been often considered in many studies in the literature. For instance, the permeability of porous media with cylinders has been investigated by Sparrow and Loeffler [START_REF] Sparrow | Longitudinal laminar flow between cylinders arranged in a regular array[END_REF], Banerjee and Hadaller [START_REF] Banerjee | Longitudinal laminar flow between cylinders arranged in a triangular array by a variational technique[END_REF], Sangani and Acrivos [START_REF] Sangani | Slow flow past periodic arrays of cylinders with application to heat transfer[END_REF], Drummond and Tahir [START_REF] Drummond | Laminar viscous flow through arrays of parallel solid cylinders[END_REF], Larson and Higdon [START_REF] Larson | Microscopic flow near the surface of twodimensional porous media, Part I: axial flow[END_REF][START_REF] Larson | Microscopic flow near the surface of twodimensional porous media, Part II: transverse flow[END_REF], Wang [START_REF] Wang | Stokes flow through a rectangular array of circular cylinders[END_REF] and Idris [START_REF] Idris | Microstructural effects on the flow law of power law fluids through fibrous media[END_REF]. The numerical implementation of the Stokes problem through an array of rigid spheres has been provided by Sangani and Acrivos [START_REF] Sangani | Slow flow through a periodic array of spheres[END_REF], Barrere et al. [START_REF] Barrere | Détermination numérique de la perméabilité en milieu poreux périodique tridimensionnel[END_REF], Chapman and Higdon [START_REF] Chapman | Oscillatory Stokes flow in periodic porous media[END_REF], Kadaksham et al. [START_REF] Kadaksham | Permeability of periodic arrays of spheres[END_REF]. Note that some works deal with the extension of the self-consistent in the field of periodic homogenization method to derive analytic expressions for the permeability [START_REF] Berdichevsky | Perform permeability prediction by self consistent method and finite element simulations[END_REF][START_REF] Boutin | Study of permeability by periodic and self consistent homogenization[END_REF] but the formula are restricted to some particular microstructures corresponding to porous media with rigid cylinders or spheres. Obviously, a unit cell containing a rigid cylinder or sphere cannot mimic realistically a porous solid and these two examples are only considered here as benchmark problems. The solutions obtained with the FFT method are compared with the results obtained by Sparrow and Loeffler [START_REF] Sparrow | Longitudinal laminar flow between cylinders arranged in a regular array[END_REF], Wang [START_REF] Wang | Stokes flow through a rectangular array of circular cylinders[END_REF] for the cylinder and those obtained by Sangani and Acrivos [START_REF] Sangani | Slow flow through a periodic array of spheres[END_REF] for the sphere. The computations are performed on a dimensionless unit cell (having the dimension 1 along each direction), leading then to the computation of a dimensionless permeability. For the FEM, we use cubic elements and the interpolation is made with the MINI elements following a method previously proposed [START_REF] Arbogast | A computatonal method for approximating a Darcy-Stokes system governing a vuggy porous medium[END_REF]. In the cubic elements, the solid or the fluid phase is assigned. All the solid voxels are eliminated in the resulting linear system which advantageously reduce the size of the problem. The adherence condition on the surface of the solid is then accounted with penalty coefficients. The total number of voxels considered with the FEM is 192 × 192 × 192 voxels, which corresponds to the maximum possible memory occupancy. A comparison with available data is provided in the case of the flow through an array of rigid cylinders or spheres, the corresponding unit is represented in Figure 1. The transversal and longitudinal components of the permeability tensor are successively com-puted by applying a pressure gradient parallel or orthogonal to the cylinder axis. The results for the transversal and longitudinal permeability are shown in figures 2 and 3, respectively. The FFT solutions are compared with those of Wang [START_REF] Wang | Stokes flow through a rectangular array of circular cylinders[END_REF], Sparrow and Loeffler [START_REF] Sparrow | Longitudinal laminar flow between cylinders arranged in a regular array[END_REF]. A good agreement is observed between our numerical values and those reported in the literature. Comparison of the FFT solution, the FEM solution and with the results obtained by Wang [START_REF] Wang | Stokes flow through a rectangular array of circular cylinders[END_REF]. R is the radius of the cylinder or the sphere, and L is the half-length of the unit cell.

Radius FFT (512 3 ) FFT (128 3 ) FEM (128 [START_REF] Sparrow | Longitudinal laminar flow between cylinders arranged in a regular array[END_REF], and Wang [START_REF] Wang | Stokes flow through a rectangular array of circular cylinders[END_REF].

The results obtained with the FFT method and the FEM are now qualitatively compared on figures 1 for the transverse flow through the array of cylinders and for various values of the Radius. The FFT-based data are provided for the resolution 128 × 128 × 128 and for 512 × 512 × 512 while, the results obtained with FEM are given for the resolution 128 × 128 × 128. We now consider numerical applications in the case of a regular array of rigid spheres. The convergence of the FFT solution with the grid refinement is provided in figure 4 for the radius R = 0.15, R = 0.25 and R = 0.35, and it is also compared with FEM solutions. A good agreement between the FFT and the FEM methods is observed for the largest value of the radius (R = 0.35).

There is a slight difference for the radius R = 0.25, that is more significant for R = 0.15. These differences are essentially attributed to the discretization error. Indeed, the sphere being approximated by means of cubic elements, the error is higher for the smaller radius. This is a well known limitation of the methods which use a regular grid. The FFT solution (with 512 × 512 × 512 voxels) and the FEM one (with 192 × 192 × 192 voxels) are compared with the data provided by Sangani and Acrivos [START_REF] Sangani | Slow flow through a periodic array of spheres[END_REF] in Table 1. Various values of the radius of the sphere have been considered in the range [0.1, 1]. It is observed that the FFT solutions are in a better agreement with the data of Sangani and Acrivos [START_REF] Sangani | Slow flow through a periodic array of spheres[END_REF] than the FEM solution. The error between the FFT method and the results of Sangani and Acrivos [START_REF] Sangani | Slow flow through a periodic array of spheres[END_REF] vary from 0.6% (in the case R = 0.1) to 0.13% (in the case R = 1). This was expected, since the FFT solution was computed with a more refined grid than the FEM solution. Obviously, the results would be compared with the same discretization, however here we compare the solutions obtained by these two methods and the maximum possible memory occupancy of the computer. This clearly shows that the FFT is computationally more attractive for reducing the memory and hence its application to the problems of high dimensions. Such comparison has been already provided in the literature but in the context of composite materials [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF]. The calculation times needed for solving the FFT method are provided on table 3. The calculation time depends on the number of iterations needed to obtain the convergence. Particularly, it must be observed that for the flow around the rigid sphere, the number of iteration increases with the radius of the sphere. 0.9 2920 (∼ 49 min) Table 3 Computation time (in seconds) with the FFT-based algorithm in the case of the flow around the rigid sphere of R.

Advanced numerical results

In this section, we consider more realistic porous cells which can mimic for instance metal foams or porous polymer solids [START_REF] Droumaguet | Engineering functional doubly porous PHEMA-based materials[END_REF]. Various configurations are obtained by an assembly of overlapping spherical pores. For most studies on composite materials, the hypothesis of non-overlapping inclusions is consid-ered, while in the present case, it is necessary to consider overlapping spheres in order to obtained a network of interconnected pores. In this section we are only interested by the evaluation of the capacity of the FFT method to handle the Stokes flow problem through complex porous microstructures and to provide some comparisons with the FEM method. The next step will concern the application to idealized reconstructed unit cells or to digital images obtained from microtomography. As for example, polymeric PHEMA-based porous materials are elaborated with spherical particulates and the porosity is generated with porogenic solvents (see for instance Le [START_REF] Droumaguet | Engineering functional doubly porous PHEMA-based materials[END_REF]. These porous materials have demonstrated some potential or promising results in the development of devices for drug delivery or scaffolds for tissue engineering applications (see for example Liu et al. 2004, Lee et al. 2007). The use of idealized tridimensional porous cell made up of overlapping spherical pores is adapted for this kind of materials and the mass transfer properties must be studied on reconstructed porous cell knowing some the characteristic size of the particulates and the porosity. As another example, the acoustic properties of metal foams are studied can be studied by performing calculations on idealized unit cells (see for instance Hoang and Perrot [START_REF] Hoang | Solid films and transports in cellular foams[END_REF][START_REF] Hoang | Identifying local characteristic lengths governing sound wave properties in solid foams[END_REF]). Metal foams are often used to increase sound absorption in aeronautics and building construction. The determination of the dynamic permeability is the cornerstone of the acoustic approaches of porous media and can be determined using the FFT approach developed in Nguyen et al. [START_REF] Nguyen | A Fourier based numerical method for computing the dynamic permeability of porous media[END_REF]. As a first example, we represent, in figures 5 and 6, a regular array of spherical overlapping voids arranged along a cubic lattice. Figure 5 represents, at the left, a unit cell having the dimension 2 and containing a body-centered cubic (BCC) array of spherical cavities with the radius R = 0.5. Each spheres intersect another one along the diagonals of the cube. The resulting irreducible unit cell of the porous medium is represented in figure 5, at the right. The porosity for the BCC array is φ = 0.939. Another example is given in figure 6. At the left, we represent a face-centered cubic (FCC) structure of spherical voids. At the right, the irreducible unit cell of the porous medium is represented for the FCC structure. The radius of the spherical voids is R = 0.375. The porosity for the FCC array is φ = 0.858. The results for the permeability are shown in figure 7 and8. The permeability, computed with the FFT method and the FEM, is plotted as function of the grid refinement. We now consider a more complex porous medium (Fig. 9) defined by spherical voids with different sizes : a larger void located at the center of the cubic cell having the radius R 1 = 0.52 and 376 smaller voids having the radius R 2 = 0.06.

Each smaller void overlaps with its neighboring void and some small voids also overlaps the larger one. By doing so, the fluid domain is then ensured to be interconnected. The porosity for this microstructure is φ = 0.915. In figure 10, we consider a second example corresponding to a random distribution of the smaller voids. The unit cell contains 536 voids with a radius comprised between 0.06 and 0.08 randomly distributed, and a large void still located at the center of the cell and having the radius R = 0.52. The porosity for this last example is φ = 0.871. The process used to construct the unit cell is the following: We build one eighth of the unit cell that we duplicate by symmetry with respect to three orthogonal planes (Ox1x2, Ox1x3 and Ox2x3 in the cartesian frame). By doing so, the resulting microstructure is then periodic. The construction of the eighth of the unit cell is made in two steps. In the first one, we put the larger sphere of radius R = 0.52 at the center of the unit cell. In the second step we put all the smaller voids and we use numerical tests. The radius of each small void is randomly considered in the range [0.06, 0.08] and its position is randomly chosen in the range [0, 0.5] × [0, 0.5] × [0, 0.5]. Numerical tests are used to ensure that for each occurrence, the smaller void cross at least one other smaller void or the larger void. This verification is needed in order to obtain an interconnected microstructure.

The dimensionless permeability for these two microstructures is represented in figures 11 and 12. In the case of the regular arrangement of smaller voids (and also for the BCC and the FCC arrays), the permeability is isotropic and is defined by only one coefficient. However, in the particular case corresponding to random distribution of smaller voids, the permeability tensor is anisotropic. In fact, in figure 12, we only represent the component K 11 by applying a pressure gradient along the direction Ox 1 . The microstructure is however quasi-isotropic, the off-diagonal coefficients of the permeability tensor are close to 0 and the terms on the diagonal take very close values. The lower convergence observed on figure 12 will probably be attributed to the tortuosity of the porous medium with random pores compared to the case of regularly aligned pores. This tortuosity would introduce some complexity in the local field that an accurate description would require higher order Fourier series expansions (that is obtained by increasing of the resolution of the image). The application of the FEM to the case of random distributions of overlapping sphere fails due to the voxelization of the image. Indeed, for such microstructure the smaller voids must contain a sufficient number of voxels such that the problem would be numerically solved. If a least the number of vowels is not sufficient, the fluid flow circulation is blocked in the smaller voids. When the FFT method is used, the results are only provided for a resolution equal or higher than 128 3 . When the FEM is used the resolution of the linear system is made with an iterative method whose convergence greatly depends on the geometry of the porous medium. Particularly, for the porous cell with a random distribution of spherical voids, the time computation has found to be prohibitive. For completeness, in Tables 2345, we also provide the values of the permeability for the flow through different microstructures depicted in figures 5, 6, 9 and 10 and for the values of resolution 128 3 , 256 3 , 384 3 and 512 3 . Just as indication, we also provide the relative error along formula given in Table 2, to evaluate the accuracy of the solution and the rate of convergence with respect to the resolution. It is observed that the relative error decreases when the resolution increases, that is expected. However, at a given spatial resolution, the error significantly differs from one microstructure to another. In particular, the the relative error is the lowest for the BCC array of pores (see Fig. 5) that is close to 0.1 percent for the spatial resolution of 512 3 . The relative error is the highest for the last example corresponding to the random distribution of small 0.004078 0.01115 Table 6 Values of permeability and relative error as function of the resolution for the problem with a large pore and smaller aligned pores voids depicted on Fig. 10 that is close to 1.2 percents for the spatial resolution of 512 3 . It is also interesting to note that that the relative error given on table 4 and 5 are quite similar, even though the results corresponds to a regular or a random distribution of spherical voids. It seems that the tortuosity, which is more important for the random distribution, do not really affect the rate of convergence with the spatial distribution. The calculation times needed are provided on table 8 for the resolution 512 × 512 × 512.

Conclusion

In the present paper, we have employed a numerical FFT-based method for computing the permeability of periodic porous media from digital images which come from modern imagery devices. Different configurations have been 8 Computation time (in seconds) with the FFT-based algorithm and for the different problems and with the resolution 512 × 512 × 512 voxels tested to show the capacity of the method and the accuracy of its solution. First, benchmark problems, consisting of regular arrays of rigid cylinders or spheres have been considered. The results have been compared to those coming from the literature and show a good accuracy of the FFT-based solution. Advanced numerical results have been then considered. Some microstructures are constructed with an assemblies of overlapping spherical voids that can mimic more realistically the real microstructure of some human-made porous materials. First, porous media are obtained by considering regularly distributed spheres along a body-centered cubic (BCC) or face-centered cubic (FCC) structures. More complex microstructure are then studied by considering porous media having different pore size which are regularly or randomly distributed. In each case, the results show a good convergence with the grid refinement used by the FFT algorithm.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Unit cell containing a rigid cylinder or sphere

Fig. 3 .

 3 Fig.3. Dimensionless longitudinal permeability for the flow through a regular array of cylinders. Comparison of the FFT solution with the results obtained by Sparrow and Loeffler[START_REF] Sparrow | Longitudinal laminar flow between cylinders arranged in a regular array[END_REF], and Wang[START_REF] Wang | Stokes flow through a rectangular array of circular cylinders[END_REF].

Fig. 4 .

 4 Fig. 4. Dimensionless permeability for the flow through a regular array of sphere of radius R = 0.15, R = 0.25 and R = 0.35. Comparison between the FFT and the FEM solution.

Fig. 5 .

 5 Fig. 5. Body-centered array of spherical voids. At the left: the voids distribution, at the right: the rigid porous solid.

Fig. 6 .

 6 Fig. 6. Face-centered array of spherical voids. At the left: the voids distribution, at the right: the rigid porous solid.

Fig. 7 .

 7 Fig. 7. Dimensionless permeability for the flow through a body-centered array of spherical voids.

Fig. 8 .

 8 Fig. 8. Dimensionless permeability for the flow through a face-centered array of spherical voids.

Fig. 9 .

 9 Fig. 9. Porous material with 376 regularly aligned small voids and one larger void located at the center of the unit cell. At the left: the pore positions, at the right: the rigid porous solid.

Fig. 10 .

 10 Fig. 10. Porous material with 536 randomly distributed small voids and one larger void located at the center of the unit cell. At the left: the pore positions, at the right: the rigid porous solid.

Fig. 11 .

 11 Fig.11. Dimensionless permeability for the flow through the microstructure given in Fig.9

Fig. 12 .

 12 Fig.[START_REF] Alcocer | Permeability of periodic porous media[END_REF]. Dimensionless permeability for the flow through the microstructure given in Fig.10

  Comparison of the intrinsic permeability of transverse flow through array of cylinders computed with the FFT method and the FEM.

				3 )
	0.1	0.134310	0.133383	0.131609
	0.2	0.081148	0.080984	0.079120
	0.3	0.051926	0.051386	0.051078
	0.4	0.032864	0.032476	0.032255
	0.5	0.019855	0.019600	0.019366
	0.6	0.010921	0.010692	0.010769
	0.7	0.005143	0.005049	0.004966
	0.8	0.001812	0.001699	0.001782
	0.9	0.000308	0.000279	0.000302
	Table 1

  Dimensionless permeability for the flow through a regular array of sphere with radius lying in the range [0.[START_REF] Auriault | Study of macroscopic behavior of a deformable porous medium[END_REF][START_REF] Auriault | Study of macroscopic behavior of a deformable porous medium[END_REF]. Comparison between the FFT solution, the FEM solution and the results of Sangani and Acrivos[START_REF] Sangani | Slow flow through a periodic array of spheres[END_REF].

	R	FFT	FEM	Sangani and Acrivos [7]
	0.1 0.905109	-	0.911070
	0.2 0.380631	-	0.382190
	0.3 0.207387 0.196734	0.208050
	0.4 0.122899	-	0.123270
	0.5 0.074448 0.070996	0.074668
	0.6 0.044340	-	0.044501
	0.7 0.025136 0.024461	0.025246
	0.8 0.013121	-	0.013197
	0.85 0.009098	-	0.009151
	0.9 0.006140	-	0.006153
	0.95 0.003980	-	0.004003
	1	0.002487	-	0.002520
	Table 2
		Radius Time(s)
		0.1	869 (∼ 14 min )
		0.2	1001 (∼ 17 min)
		0.3	1126 (∼ 19 min)
		0.4	1139 (∼ 19 min)
		0.5	1190 (∼ 20 min)
		0.6	1404 (∼ 23 min)
		0.7	1691 (∼ 28 min)
		0.8	2002 (∼ 33 min)

  Values of permeability and relative error as function of the resolution for the problem with a large pore and smaller randomly distributed pores

	Resolution permeability relative error
	128 3	0.009423	-
	256 3	0.007296	0.12722
	384 3	0.007032	0.01842
	512 3	0.006856	0.01267
	Table 7
	Configurations		Time(s)
	BCC		144149 (∼ 1.7 day)
	FFC		81317 (∼ 0.94 day)
	376 regularly aligned	163826 (∼ 1.9 day)
	536 randomly distributed 643283 (∼ 7.44 days)
	Table