Vincent Monchiet
email: vincent.monchiet@univ-mlv.fr

Guy Bonnet
email: guy.bonnet@univ-mlv.fr

On the inversion of non symmetric sixth-order isotropic tensors and conditions of positiveness of third-order tensor valued quadratic functions

Keywords: Sixth-Order Tensors, Third-Order Tensors, Isotropic Tensors, Microstructured Media

In the present paper we propose new results concerning linear tensorial algebra for third-order and non symmetric isotropic sixth-order tensors in the most general case (i.e. having not the major and minor symmetries). Such tensors are used, for instance, in the theory of microstructured elastic media. A formalism based on an irreducible basis for isotropic sixth-order tensors is introduced, which is useful for performing the classical tensorial operations. Specially, a condensed expression for the product between two isotropic sixth-order tensors is provided, which allows the obtaining of a condition on these tensors for being invertible and a closed form expression of the inverse of such a tensor. Finally, the condition of positiveness of third-order tensor-valued quadratic functions is derived. For instance, such conditions are required for computing the elastic energy of microstructured media.

Introduction

Generalized continuum theories use tensors of-order 2, 3... as variables as shown for instance by Toupin [START_REF] Toupin | Elastic materials with couple-stresses[END_REF] or Mindlin [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF], Mindlin and Eshel [START_REF] Mindlin | On first gradients theory in linear elasticity[END_REF]. In these theories, the constitutive equations for elasticity take the form of a linear mapping between tensors of degree 2, 3... which involves the introduction of elastic tensors whose degrees are 4, 5, 6.... The use of such high-order tensors raises a number of fundamental questions, specially concerning their inversion or the condition of positiveness when computing the elastic energy from these tensors. Here we restrict our analysis to the theories which use third-order tensors. For instance, the third gradient theory uses the gradient of strain, κ ijk = ∇ k ε ij (see [START_REF] Mindlin | On first gradients theory in linear elasticity[END_REF]) whereas the theory of microstructured media uses the gradient of microstrain χ ijk = ∇ k ψ ij (see [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF]). When an hypothesis of isotropy and centrosymmetry is used, the elastic constitutive equations for the both theories contain a linear relation between two third-order tensors, denoted a and b, which can be put into the form:

b = A ⊙ 3 a ( 1
)
Where A is a constant sixth-order tensor, namely independent of a. In (1), the symbol ⊙ 3 represents the triple contraction between two tensors x and y: [x ⊙ 3 y] i..j = x i..pqr y pqr..j . For instance, in [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF], the components of b are b ijk = A ijkpqr a pqr . Two particular cases of (1) can be considered: It corresponds to the case of microstructured elastic media, a ≡ ∇ψ and b is called double stress. It follows that A ijkpqr has no more the minor symmetries. In a recent paper [START_REF] Monchiet | Inversion of higher order isotropic tensors with minor symmetries and solution of higher order heterogeneity problems[END_REF], an irreducible basis for isotropic 2n th -order tensors (with n = 3, 4, ...) having minor symmetries has been provided. It has been shown that the formalism introduced in this work is useful for doing the classical tensorial operations and specially the inversion of symmetric isotropic tensors of degree up to 8. However, the isotropic tensors considered in this study possessed the minor symmetries which restrict the domain of application of this formalism to the case of gradient theories. As it will be shown thereafter, a key to study the algebra is the decomposition between spherical and deviatoric parts. In this context, it is worthwhile to notice some studies which deal with the decomposition of third-order tensors and third-order tensor valued functions (see for instance Pennisi & Trovato [4], Pennisi [START_REF] Pennisi | On third order tensor-valued isotropic functions[END_REF], Smyshlyaev & Fleck [START_REF] Smyshlyaev | The role of strain gradients in the grain size effect for polycrystals[END_REF], Smith [START_REF] Smith | Isotropic tensor-valued polynomial function of second and third-order tensors[END_REF]). However, all those studies dealt only with symmetric third-order tensors. In the present paper, we aim at generalizing the results provided in [START_REF] Monchiet | Inversion of higher order isotropic tensors with minor symmetries and solution of higher order heterogeneity problems[END_REF] to the case of a sixth-order tensor which does not posses minor symmetries. A representation of such tensors in a basis constituted of 15 non symmetric isotropic tensors is provided in section 2. This basis is shown to be useful for deriving a closed form expression of the inverse of a non symmetric isotropic sixth-order tensor (namely having not the minor and the major symmetries). In section 3, we derive the condition of positiveness of A in the case of tensors having the major symmetry (symmetry between the first three indices and the last three indices), which can be written:

-(i)
∀a ̸ = 0 : a ⊙ 3 A ⊙ 3 a > 0 (2)
Here, we deal with tensors having the major symmetries (but not the minor symmetries). Condition (2) is needed, for instance, for obtaining the condition of positiveness of the elastic energy of microstructured media.

An irreducible basis for isotropic sixth-order tensors

In [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF], the following decomposition for an isotropic sixth-order tensor A is used:

A = n=15 ∑ n=1 a n T n (3)
where the components of tensors T n for n = 1..15 are given by:

[T 1 ] ijkpqr = δ ip δ jq δ kr , [T 2 ] ijkpqr = δ ir δ jp δ kq , [T 3 ] ijkpqr = δ iq δ jr δ kp [T 4 ] ijkpqr = δ ir δ jq δ kp , [T 5 ] ijkpqr = δ ip δ jr δ kq , [T 6 ] ijkpqr = δ iq δ jp δ kr [T 7 ] ijkpqr = δ ij δ kr δ pq , [T 8 ] ijkpqr = δ ij δ kp δ qr , [T 9 ] ijkpqr = δ ij δ kq δ pr [T 10 ] ijkpqr = δ ik δ jr δ pq , [T 11 ] ijkpqr = δ ik δ jp δ qr , [T 12 ] ijkpqr = δ ik δ jq δ pr [T 13 ] ijkpqr = δ jk δ ir δ pq , [T 14 ] ijkpqr = δ jk δ ip δ qr , [T 15 ] ijkpqr = δ jk δ iq δ pr (4)
When a tensor A possesses the major symmetry, A ijkpqr = A pqrijk , the following conditions for the a n are required: a 2 = a 3 , a 8 = a 13 , a 9 = a 10 , a 11 = a 15 . In addition, the conditions for A to possess the minor symmetries are detailed in [START_REF] Monchiet | Inversion of higher order isotropic tensors with minor symmetries and solution of higher order heterogeneity problems[END_REF].

The 15 tensors T n for n = 1..15 constitute an irreducible basis for any isotropic sixthorder tensors. However, this basis appears to be not useful for obtaining the inverse of A since T n ⊙ 3 T m ̸ = 0 whatever the values of n and m. The use of basis T n for the inversion of a sixth-order tensor leads then to a full linear system of dimension 15. For this reason, we propose a new basis denoted (K n , J m ) which is more convenient for inverting A. This basis is constituted of sixth-order tensors denoted by K n for n = 1..6, which are:

K 1 = 1 6 [ T 1 + T 2 + T 3 -T 4 -T 5 -T 6 ] K 2 = 1 6 [ T 1 + T 2 + T 3 + T 4 + T 5 + T 6 ] - 1 15
[ T 7 + T 8 + T 9 + T 10 + T 11 + T 12 + T 13 + T 14 + T 15 ] K 3 = 1 6 [ T 10 + T 11 -T 13 -T 14 ] + 1 3 [ T 1 + T 4 + T 15 -T 3 -T 6 -T 12 ] K 4 = 1 6 [ T 11 + T 12 -T 14 -T 15 ] + 1 3 [ T 3 + T 5 + T 13 -T 2 -T 4 -T 10 ] K 5 = 1 6 [ T 7 + T 8 -T 13 -T 14 ] + 1 3 [ T 2 + T 5 + T 15 -T 3 -T 6 -T 9 ] K 6 = 1 6 [ T 8 + T 9 -T 14 -T 15 ] + 1 3 [ T 1 + T 6 + T 13 -T 2 -T 4 -T 7 ] (5)
and 9 tensors denoted J m for m = 1..9 which are defined by:

J 1 = 1 10 [ 4T 14 -T 13 -T 15 ] , J 2 = 1 10 [ 4T 15 -T 14 -T 13 ] J 3 = 1 10 [ 4T 13 -T 15 -T 14 ] , J 4 = 1 10 [ 4T 11 -T 10 -T 12 ] J 5 = 1 10 [ 4T 12 -T 10 -T 11 ] , J 6 = 1 10 [ 4T 10 -T 11 -T 12 ] J 7 = 1 10 [ 4T 8 -T 7 -T 9 ] , J 8 = 1 10 [ 4T 9 -T 7 -T 8 ] , J 9 = 1 10 [ 4T 7 -T 8 -T 9 ] (6)
It can be shown that the J n for n = 1..9 and the K n for n = 1..6 have the property

K n ⊙ 3 J m = J m ⊙ 3 K n = 0
whatever the values of n and m. The triple contraction between two elements taken from (J 1 , ..., J 9 ) are given in table 1. The triple contraction between two elements taken from (K 1 , ..., K 6 ) are given in table 2.

⊙ 3 J 1 J 2 J 3 J 4 J 5 J 6 J 7 J 8 J 9 J 1 J 1 J 2 J 3 0 0 0 0 0 0 J 2 0 0 0 J 1 J 2 J 3 0 0 0 J 3 0 0 0 0 0 0 J 1 J 2 J 3 J 4 J 4 J 5 J 6 0 0 0 0 0 0 J 5 0 0 0 J 4 J 5 J 6 0 0 0 J 6 0 0 0 0 0 0 J 4 J 5 J 6 J 7 J 7 J 8 J 9 0 0 0 0 0 0 J 8 0 0 0 J 7 J 8 J 9 0 0 0 J 9 0 0 0 0 0 0 J 7 J 8 J 9
Table 1: The triple contraction between the J n for n = 1..9

⊙ 3 K 1 K 2 K 3 K 4 K 5 K 6 K 1 K 1 0 0 0 0 0 K 2 0 K 2 0 0 0 0 K 3 0 0 K 3 K 4 0 0 K 4 0 0 0 0 K 3 K 4 K 5 0 0 K 5 K 6 0 0 K 6 0 0 0 0 K 5 K 6
Table 2: The triple contraction between the K n for n = 1..6

Denoting by E 6 the space of isotropic sixth-order tensors, the basis (K n , J m ) constitutes an irreducible basis for (E 6 , ⊙ 3 , I) where I = T 1 is the identity for the composition ⊙ 3 . Moreover, it can be observed that (E 6 , ⊙ 3 , I) define a monoid, namely an algebraic structure with a single associative binary operation and an identity element. Every tensor A ∈ E 6 is then decomposed as follows:

A = α 1 K 1 + α 2 K 2 + α 3 K 3 + α 4 K 4 + α 5 K 5 + α 6 K 6 + α 7 J 1 + α 8 J 2 + α 9 J 3 + α 10 J 4 + α 11 J 5 + α 12 J 6 + α 13 J 7 + α 14 J 8 + α 15 J 9 (7)
α 1 , ..., α 15 are the components of A in the (K n , J m ) basis. The base change relations, giving coefficients α n as function of the a n are provided in appendix A.

Denoting by J the sub-space of isotropic sixth-order tensors given by A = α 7 J 1 + ... + α 15 J 9 (or α 1 = α 2 = ... = α 6 = 0 in ( 7)), it can be observed that (J , ⊙ 3 , J) define a submonoid, here J = J 1 + J 5 + J 9 is the unit element of J for the composition ⊙ 3 . Similarly, we denote by K the sub-space of isotropic sixth-order tensors given by A = α 1 K 1 + ... + α 6 K 6 . It can be observed that (K, ⊙ 3 , K) define a submonoid, K =

K 1 + K 2 + K 3 + K 6
being the unit element of K for the composition ⊙ 3 . It follows that any isotropic sixth-order tensors can be decomposed into A = A J + A K where A J = J ⊙ 3 A ⊙ 3 J and A K = K ⊙ 3 A ⊙ 3 K and that the identity for the triple contraction is given by I = K + J. This decomposition will be useful for providing various results and particularly for the product between two sixth-order tensors or the inversion of sixthorder tensors. Note that the decomposition (K n , J m ) gives rise to a generalization of the definition of the spherical and deviatoric parts of a third-order tensors, commonly used for two-order tensors, as it will be discussed in section 3. The basis (K n , J m ) contains two tensors having the major symmetry, which are K 1 and K 2 . Other tensors do not possess the major symmetry. The condition for tensor A having the major symmetry is:

{ α 3 + 2α 5 = α 6 +
Let us consider now a sixth-order tensor A ∈ E 6 and let us denote by α 1 , ..., α 15 its components in the basis (K n , J m ), as in [START_REF] Smyshlyaev | The role of strain gradients in the grain size effect for polycrystals[END_REF]. Let us do the same with a second isotropic sixth-order tensor B and let us denote by β 1 , ..., β 15 its components in (K n , J n ). We aim at computing the components of tensor C defined by C = B ⊙ 3 A. As explained above, any sixth-order tensor C can be decomposed along the two orthogonal subspaces J and K as follows:

C = C J + C K . It follows that C K = B K ⊙ 3 A K and C J = B J ⊙ 3 A J .
The computation of the components of C K and C J can be performed by using tables 1 and 2. We denote by γ n the components of C in the basis (K n , J n ). These components are given by:

                       γ 1 = β 1 α 1 , γ 2 = β 2 α 2 , γ 3 = β 3 α 3 + β 4 α
We now look at the inverse of tensor A, namely we search B which verify the equation

B ⊙ 3 A = A ⊙ 3 B = I
where it is recalled that I = K + J is the identity for the triple contraction. Let us recall that the decomposition of J and K in the (K n , J m ) basis, reads

K = K 1 + K 2 + K 3 + K 6 , J = J 1 + J 5 + J 9 .
Consequently, the linear system giving the components of B can be obtained from (9) in which we put γ 1 = γ 2 = γ 3 = γ 6 = γ 7 = γ 11 = γ 15 = 1 and all other γ n = 0. The computation of the β n for n = 1..15 can be performed by a few simple algebraic equations. It leads to: where Ω and ∆ are given by:

                       β 1 = 1 α 1 , β 2 = 1 α 2 , β 3 = α 6 Ω , β 4 = - α 4 Ω , β 5 = - α 5 Ω , β 6 = α 3 Ω β 7 =
Ω = α 3 α 6 -α 4 α 5 ∆ =
With the following condition for A having an inverse: α 1 α 2 Ω∆ ̸ = 0.

Condition of positiveness of the elastic energy

As it is well known from thermodynamic considerations, it is crucial to ensure that the elastic energy is positive definite1 . Such a condition can be expressed easily for classically used isotropic second or fourth-order tensors, but it is not obvious to derive similar expressions for the case of higher degree elasticity. We now propose to determine the conditions of positiveness of the quadratic form W (a) = a ⊙ 3 A ⊙ 3 a = b ⊙ 3 a where b = a⊙ 3 A. Tensor A is now assumed to possess the major symmetry (A ijkpqr = A pqrijk ). As a consequence, its components in the basis (K n , J m ) comply to relations [START_REF] Toupin | Elastic materials with couple-stresses[END_REF]. Let us first introduce the spherical part of a third-order a, denoted by S(a), and the deviatoric part of a, denoted by D(a):

D(a) = K ⊙ 3 a, [D(a)] ijk = a ijk -[θ 1 ] i δ jk -[θ 2 ] j δ ik -[θ 3 ] k δ ij S(a) = J ⊙ 3 a, [S(a)] ijk = [θ 1 ] i δ jk + [θ 2 ] j δ ik + [θ 3 ] k δ ij (12)
where the components of θ n for n = 1, 2, 3 are given by:

[θ 1 ] i = 1 10 [ 4a ipp -a pip -a ppi ] , [θ 2 ] i = 1 10 [ 4a pip -a ipp -a ppi ] [θ 3 ] i = 1 10 [ 4a ppi -a ipp -a pip ] (13)
and it is recalled that J = J 1 + J 5 + J 9 and K = K 1 + K 2 + K 3 + K 6 . Effecting the same operations with the third-order tensor b, it appears that W (a) can be decomposed into W (a) = W D (a) + W S (a) where:

W D (a) = D(b) ⊙ 3 D(a), W S (a) = S(b) ⊙ 3 S(a) ( 14
)
where b is given by [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF]. It follows that the positiveness of W (a) is ensured if and only if D(b) ⊙ 3 D(a) ≥ 0 and S(b) ⊙ 3 S(a) ≥ 0. We now introduce the six deviatoric sub-parts of a third-order tensor a, denoted by D n (a) = K n ⊙ 3 a for n = 1..6, whose components are given by:

[D 1 (a)] ijk = ε ijk χ, [D 2 (a)] ijk = κ ijk , [D 3 (a)] ijk = ε ijp [χ 1 ] pk [D 4 (a)] ijk = ε ijp [χ 2 ] pk , [D 5 (a)] ijk = ε ikp [χ 1 ] pj , [D 6 (a)] ijk = ε ikp [χ 2 ] pj (15)
where ε ijk is the permutation symbol. The scalar χ and tensors χ 1 , χ 2 ,κ are given by:

κ ijk = 1 6 { [D(a)] ijk + [D(a)] kij + [D(a)] jki + [D(a)] kji + [D(a)] ikj + [D(a)] jik } χ = 1 6 ε ijk [D(a)] ijk [χ 1 ] ij = 1 6 ε ipq ( [D(a)] pqj + [D(a)] jqp ) + 1 6 ε jpq ( [D(a)] pqi + [D(a)] iqp ) [χ 2 ] ij = 1 6 ε ipq ( [D(a)] jpq + [D(a)] pjq ) + 1 6 ε jpq ( [D(a)] ipq + [D(a)] piq ) (16)
χ 1 and χ 2 are two symmetric, traceless, second-order tensors whereas κ is a symmetric traceless third-order tensor (namely κ ijk is invariant by any permutation of its indices i, j, k and κ ijj = κ iji = κ jji = 0). Tensors χ 1 , χ 2 , κ are respectively defined by 5, 5 and 7 independent coefficients. Similarly, the triple contraction between the J n and a produces the 9 spherical sub-parts of a, denoted S n (a) = J n ⊙ 3 a. Their components are given by:

[S 1 (a)] ijk = [θ 1 ] i δ jk , [S 2 (a)] ijk = [θ 2 ] i δ jk , [S 3 (a)] ijk = [θ 3 ] i δ jk [S 4 (a)] ijk = [θ 1 ] j δ ik , [S 5 (a)] ijk = [θ 2 ] j δ ik , [S 6 (a)] ijk = [θ 3 ] j δ ik [S 7 (a)] ijk = [θ 1 ] k δ ij , [S 8 (a)] ijk = [θ 2 ] k δ ij , [S 9 (a)] ijk = [θ 3 ] k δ ij (17)
Where θ 1 , θ 2 , θ 3 are given by (13).

As quoted in the previous section, I is the identity for the composition ⊙ 3 , namely: a = I ⊙ 3 a. Since the decomposition of I in the basis (K n , J m ) is

I = K 1 + K 2 + K 3 + K 6 + J 1 + J 5 + J 9 , it
                                   D 1 (b) = α 1 D 1 (a) D 2 (b) = α 2 D 1 (a)
Again, table 1 and 2 must be used to derive expressions above. Taking into account (18), (19) with definitions (15) and ( 17), it follows that:

W D (a) = 6α 1 χ 2 + α 2 κ ⊙ 3 κ + φ 1 χ 1 ⊙ 2 χ 1 + φ 2 χ 2 ⊙ 2 χ 2 + 2φ 3 χ 1 ⊙ 2 χ 2 W S (a) = ψ 1 θ 1 ⊙ 1 θ 1 + ψ 2 θ 2 ⊙ 1 θ 2 + ψ 3 θ 3 ⊙ 1 θ 3 + 2ψ 4 θ 2 ⊙ 1 θ 3 +2ψ 5 θ 1 ⊙ 1 θ 3 + 2ψ 6 θ 1 ⊙ 1 θ 2 (20)
Where coefficients φ n and ψ n are given by:

φ 1 = 2α
The condition of positiveness of the elastic potential is then given by: α 1 ≥ 0, α 2 ≥ 0, φ 1 ≥ 0, φ 2 ≥ 0, Ψ 1 ≥ 0, Ψ 2 ≥ 0, Ψ 3 ≥ 0 (22

)
where Ψ 1 , Ψ 2 , Ψ 3 are the eigenvalues of the 3 × 3 matrix:

Ψ =      ψ 1 ψ 6 ψ 5 ψ 6 ψ 2 ψ 4 ψ 5 ψ 4 ψ 3      (23)

Conclusion

In the present paper, we provide an irreducible basis for sixth-order isotropic tensors. This basis is constituted of 15 elements denoted K n for n = 1, 6 and J n for n = 1, 9 and it appears to be useful for doing the classical tensorial operations. For instance, a closed form expression of the inverse of an isotropic sixth-order tensor has been derived in section 3. Moreover, the projection of a third-order tensor along the tensors of the new basis (K n , J m ), gives rise to a generalization of a spherical and deviatoric parts commonly used for two-order tensors to the case of third-order tensors. This decomposition has been subsequently used for deriving the condition of positiveness of the quadratic form a ⊙ 3 A ⊙ 3 a. For instance, this last result can be used for giving the condition of positiveness of the elastic energy of microstructured elastic media.

  tensors a and b are symmetric with respect to their first two indices, namely: a ijk = a jik and b ijk = b jik . It corresponds to the case of gradient elastic theory (since the strain is a symmetric tensor), a ≡ ∇ε and b is called hyperstress. It follows that a and b are defined by only 18 independent components. The sixth-order tensor A has then the following minor symmetries: A ijkpqr = A jikpqr and A ijkpqr = A ijkqpr . -(ii) tensors a and b have no symmetry and are defined by 27 independent coefficient.

  follows that any third-order tensor a can be decomposed into:a = D 1 (a) + D 2 (a) + D 3 (a) + D 6 (a) = D(a) + S 1 (a) + S 5 (a) + S 9 (a) = S(a) (18)Note that a is defined by 27 independents components whereas tensors D 1 (a), D 2 (a), D 3 (a), D 6 (a), S 1 (a), S 5 (a), S 9 (a) are respectively defined by 1, 7, 5, 5, 3, 3, 3 independents components. The decomposition (18) can be also applied for b and the linear relation (1) can then be put into the form:

D 3 (

 3 b) = α 3 D 3 (a) + α 4 D 4 (a) D 6 (b) = α 5 D 5 (a) + α 6 D 6 (a) S 1 (b) = α 7 S 1 (a) + α 8 S 2 (a) + α 9 S 3 (a) S 5 (b) = α 10 S 4 (a) + α 11 S 5 (a) + α 12 S 6 (a) S 9 (b) = α 13 S 7 (a) + α 14 S 8 (a) + α 15 S 9 (a)

  2α 4 , α 7 + 3α 10 + α 13 = 3α 8 + α 11 + α 14 α 8 + α 11 + 3α 14 = α 9 + 3α 12 + α 15 , 3α 9 + α 12 + α 15 = α 7 + α 10 + 3α 13

  α 12 α 14 -α 11 α 15 ∆ , β 8 = α 8 α 15 -α 9 α 14 ∆ , β 9 = α 9 α 11 -α 8 α 12 ∆

	β 10 =	α 10 α 15 -α 12 α 13 ∆	, β 11 =	α 9 α 13 -α 7 α 15 ∆	, β 12 =	∆ α 7 α 12 -α 9 α 10	(10)
	β 13 =	α 11 α 13 -α 10 α 14 ∆	, β 14 =	α 7 α 14 -α 8 α 13 ∆	, β 15 =	α 8 α 10 -α 7 α 11 ∆

  α 7 α 12 α 14 + α 8 α 10 α 15 + α 9 α 11 α 13 -α 7 α 11 α 15 -α 8 α 12 α 13 -α 9 α 10 α 14

  3 + α 5 , φ 2 = α 4 + 2α 6 , φ 3 = α 3 + 2α 5 = 2α 4 + α 6 ψ 1 = 3α 7 + α 10 + α 13 , ψ 2 = α 8 + 3α 11 + α 14 , ψ 3 = α 9 + α 12 + 3α 15 ψ 4 = α 8 + α 11 + 3α 14 = α 9 + 3α 12 + α 15ψ 5 = 3α 9 + α 12 + α 15 = α 7 + α 10 + 3α 13 ψ 6 = α 7 + 3α 10 + α 13 = 3α 8 + α 11 + α 14

When the case of isotropy is considered, the elastic energy in[START_REF] Mindlin | Micro-structure in linear elasticity[END_REF] can be decomposed into two independent terms corresponding to a two-order and a third-order tensor valued quadratic functions. Here we only focus our attention on the third-order valued quadratic function

Appendix A. Appendix : Base change relations

The components α i of a 6 th order tensor along the basis (J n , K m ) from its components a j along the bases (T n ) are given by:

a 1 + a 5 + a 13 + 3a 14 + a 15 ; α 8 = a 3 + a 6 + a 13 + a 14 + 3a 15 α 9 = a 2 + a 4 + 3a 13 + a 14 + a 15 ; α 10 = a 2 + a 6 + a 10 + 3a 11 + a 12 α 11 = a 1 + a 4 + a 10 + a 11 + 3a 12 ; α 12 = a 3 + a 5 + 3a 10 + a 11 + a 12 α 13 = a 3 + a 4 + a 7 + 3a 8 + a 9 ; α 14 = a 2 + a 5 + a 7 + a 8 + 3a 9 α 15 = a 1 + a 6 + 3a 7 + a 8 + a 9 (A. [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF]


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

Guy Bonnet
email: guy.bonnet@univ-paris-est.fr

A polarization based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast

Keywords: Elastic composites, Fast Fourier transform, Polarization, Voids, Rigid inclusions

It is recognized that the convergence of FFT based iterative schemes used for computing the effective properties of elastic composite materials drastically depends on the contrast between the phases. Particularly, the rate of convergence of the strain based iterative scheme strongly decreases when the composites contain very stiff inclusions and the method diverges in the case of rigid inclusions. Reversely, the stress based iterative scheme converges rapidly in the case of composites with very stiff or rigid inclusions, but leads to low convergence rates when soft inclusions are considered and to divergence for composites containing voids. It follows that the computation of effective properties is costly when the heterogeneous medium contains simultaneously soft and stiff phases. Particularly the problem of composites containing voids and rigid inclusions cannot be solved by the strain or the stress based approaches. In this paper, we propose a new polarization-based iterative scheme for computing the macroscopic properties of elastic composites with an arbitrary contrast which is nearly as simple as the basic schemes (strain and stress based) but which has the ability to compute the overall properties of multiphase composites with arbitrary elastic moduli, as illustrated through several examples.

Introduction

An alternative approach to FEM and BEM methods [START_REF] Helsing | An integral equation method for elastostatics of periodic composites[END_REF][START_REF] Gusev | Representative volume element size for elastic composites: a numerical study[END_REF][START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF] was proposed in the middle of the nineties [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] for the computation of effective properties of linear elastic periodic composite. The local cell problem is solved by means of an iterative scheme which uses the periodic Green's tensor for the strain and exact images of the microstructures. The main advantages of this approach over other existing ones can be summarized as follows: (i) the method is based on a recurrence relation and does not lead to the resolution of a linear system, (ii) the discrete Fourier transform and its inverse are computed with the Fast Fourier Transform (FFT) which significantly increases the performance of the method, (iii) the memory needed for solving the problem is significantly reduced compared to other methods. The approach has been afterwards extended to non linear composites [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF]. It have been used by the authors themselves but also by other researchers for computing the overall response of visco-elasto-plastic composites (see for instance [START_REF] Bilger | Effect of a nonuniform distribution of voids on the plastic response of voided materials: a computational and statistical analysis[END_REF][START_REF] Bilger | Bounds and estimates for the effective yield surface of porous media with a uniform or a nonuniform distribution of voids[END_REF][START_REF] Lebensohn | Orientation image-based micromechanical modelling of subgrain texture evolution in polycrystalline copper[END_REF][START_REF] Idiart | Infinite-contrast periodic composites with strongly nonlinear behavior: Effective-medium theory versus full-field simulations[END_REF][START_REF] Lee | Modeling the Viscoplastic Micromechanical Response of Two-Phase Materials using Fast Fourier Transforms[END_REF]). 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 elastic solution into Neumann series, along the lines of a method which was first introduced for composite conductors by Brown [START_REF] Brown | Solid mixture permittivities[END_REF] and later by Kroener [START_REF] Kroner | Statistical Continuum Mechanics[END_REF]. The convergence of this series depends on the elastic properties of the phases but also on the choice of a reference medium. The conditions which ensure the convergence of the Neumann series have been derived in [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF] (see also [START_REF] Milton | The Theory of Composites[END_REF] for a study of the convergence in the case of the thermal conduction). It has been found that, from a theoretical point of view, the strain based iterative scheme does not converge in both limited cases of voids and rigid inclusions. New conditions for the convergence of the strain based iterative scheme are derived in the present paper. An interesting modification of the results provided in [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF] is given in the paper: it is shown that the convergence can be ensured for composites with voids, as it is confirmed by comparison with numerical calculations. From another point of view, the application of the strain based method for composites containing very stiff inclusions is computational costly and the limited case of rigid inclusions cannot be solved by the method. In the domain of high stiffness of inclusions, the stress based iterative scheme [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] is shown to be better suited. However, in the region of low stiffness, the calculation times become very large and the case of void cannot be solved by this dual method. It suggests that the basic schemes are not adapted for solving some specific problems for which the heterogeneous medium contains simultaneously both stiff and soft phases (some examples of such specific problems are given in the last part of the paper). To circumvent the incapacity of the original approach [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF], other iterative schemes have been developed during the ten past years. An "accelerated" scheme has been proposed by Eyre and Milton [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF] for conductivity problems and has been generalized to elasticity in [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF] (see also [START_REF] Vinogradov | An accelerated FFT algorithm for thermoelastic and non-linear composites[END_REF] for its application to nonlinear periodic composites). The accelerated scheme improves the rate of convergence of the strain based approach for large values of the contrast. However, similarly to the original scheme [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF], the number of iterations at convergence still linearly depends on the contrast (see [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF]). As a consequence the method still diverges for rigid inclusions. To address the problems of composites with infinite contrast, alternative methods, based on augmented Lagrangians, have been proposed by Michel et al. [START_REF] Michel | A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast[END_REF][START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF]. In [START_REF] Michel | A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast[END_REF], a Lagrangian based on the local strain potential was introduced to treat the case of voids. In [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF], a Lagrangian based on the local stress potential is used to solve the problem of an elastic composite reinforced by rigid inclusions. More recently, Brisard and Dormieux [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF] have proposed a new method based on the Hashin-Shtrikman variational principle [START_REF] Hashin | On some variational principles in anisotropic and nonhomogeneous elasticity[END_REF][START_REF] Willis | Bounds and self-consistent estimates for the overall properties of anisotropic composites[END_REF] and used the conjugate gradient method for solving the problem. We emphasize that the use of conjugate gradient has been also used in the recent paper of Zeman et al. [START_REF] Zeman | Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients[END_REF] for improving the convergence of the strain based iterative scheme.

In the present paper we propose a simple iterative scheme, based on the polarization, which uses the same ingredients as the basic schemes. The paper contains three parts. In section 2, the equations of the local problem are recalled and new compact notations for representing compatibility and equilibrium equations are introduced: they use a representation by means of the Walpole basis [START_REF] Walpole | Elastic behavior of composite materials: theoretical foundations[END_REF]. That representation appears also adequate for putting the Green's tensors into a more compact form. These notations are used to propose new convergence tests for the basic schemes. The second part of the paper is dedicated to recalling the basic schemes. The rate of convergence and the limitations of the methods are discussed and illustrated through a few numerical calculations. In the last part of the paper, we present the new polarizationbased method. Its relevance over the basic schemes is shown through several applications for which the studied heterogeneous medium contains simultaneously very soft and very stiff inclusions and also voids and rigid inclusions.

The cell problem for linear elastic composites with prescribed macroscopic strain or stress

We consider a periodic composite defined by a parallelepipedic unit cell and three (two for plane strain 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 (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 α (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 ). Due to periodicity, the strain and stress fields are expanded into Fourier series:

ε(x) = ∑ n ε(ξ n ) exp(iξ n .x) σ(x) = ∑ n σ(ξ n ) exp(iξ n .x) (3)
where the ξ n are the discrete wave vectors arranged along a discrete network having a period 2π/b i in the direction x i . In (3), ε(ξ) and σ(ξ) denote the Fourier transforms of the strain and stress fields:

ε(ξ) = 1 V ∫ V ε(x) exp(-iξ.x)dx with : ε(0) = E σ(ξ) = 1 V ∫ V σ(x) exp(-iξ.x)dx with : σ(0) = Σ (4)
The difference between quantities in Fourier space and in real space is denoted by adding the variable, for example ε(ξ) for the Fourier transform of ε(x) with the same symbol "ε" in both cases. Now, a new tensorial formalism is introduced, which gives a geometrical interpretation of the compatibility and equilibrium equation in the Fourier space. This formalism uses the Walpole basis and will also have the advantage to give compact forms to Green's operators which are used within FFT iterative solution. In a first step , it is important to notice that all operators such as compatibility operator, Green's tensors,.. are ξ dependent fourth order tensors which are transversely isotropic around wave vector ξ. Therefore, the Walpole basis [START_REF] Walpole | Elastic behavior of composite materials: theoretical foundations[END_REF] can be used to represent all these tensors. Let us therefore introduce the fourth order tensors P(ξ) and Q(ξ) defined for every ξ ̸ = 0 by:

Q(ξ) = E 1 (ξ) + E 3 (ξ), P(ξ) = E 2 (ξ) + E 4 (ξ)
in which the E n (ξ) for n = 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 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 (5)
where k and k ⊥ are given by:

k = 1 |ξ| 2 ξ ⊗ ξ, k ⊥ = i -k ( 6
)
where i is the two order identity tensor. Since the E n (ξ) for n = 1..4 are projectors, P(ξ) and Q(ξ) are also two projectors:

∀ξ ̸ = 0 :              P(ξ) : P(ξ) = P(ξ), Q(ξ) : Q(ξ) = Q(ξ) P(ξ) : Q(ξ) = Q(ξ) : P(ξ) = 0 I = P(ξ) + Q(ξ) (7)
The action of Q(ξ) on a second order tensor a(ξ), having a zero volume average, gives the projection of a(ξ) along the plane of normal ξ denoted Q ξ . The action of P(ξ) on a second order tensor a(ξ) provides the "out of plane" components, namely along P ξ the supplementary of Q ξ . The two orthogonal subspaces P ξ and Q ξ have both dimension 3 for three dimensional problems. In 2D cases, E 3 (ξ) = 0, the subspaces P ξ and Q ξ are of dimension 2 and 1 respectively. The compatibility and equilibrium conditions are differential equations which become linear algebraic equations in Fourier space. They can then be represented by algebraic equations which use tensors P(ξ) and Q(ξ). For instance, the compatibility for every periodic field ε(x) is:

∀ξ ̸ = 0 : P(ξ) : ε(ξ) = ε(ξ) or Q(ξ) : ε(ξ) = 0 (8)
From another hand, the equilibrium condition for the periodic stress field σ(x), given by (3), is:

∀ξ ̸ = 0 : P(ξ) : σ(ξ) = 0 or Q(ξ) : σ(ξ) = σ(ξ) (9)
It is worthwile mentioning that, in this last case, it can be surprising to use a fourth order tensor to express equilibrium equations, while equilibrium equations lead primarily to a third order tensor. However, it is easy to show that a periodic field σ is in equilibrium iff P(ξ) : σ(ξ) = 0. The subspace of uniform fields constitute a particular case in the sense that such fields comply with both compatibility and equilibrium equations. Note also that a similar description, using decomposition into orthogonal subspaces, was also proposed by Milton [START_REF] Milton | The Theory of Composites[END_REF] in the context of conduction.

To summarize, the elastic solution of ( 1) is searched under the form (3), in which ε(ξ n ) and σ(ξ n ) comply with conditions ( 8) and ( 9).

3 The basic iterative schemes

The strain and stress based iterative schemes

An iterative scheme has been proposed by Moulinec et al. [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] for computing the solution of the linear elastic problem (1), on the basis of previously used similar schemes introduced for conduction problems. The approach introduces a reference medium of rigidity C 0 to re-write the system of equations (1) as a Lippmann-Schwinger type equation. The solution of that integral equation is expanded along Neumann series [START_REF] Brown | Solid mixture permittivities[END_REF][START_REF] Kroner | Statistical Continuum Mechanics[END_REF]. Each term of the series is obtained by the following recurrence relations:

                         ε i (x) = F -1 (ε i (ξ)) σ i (x) = C(x) : ε i (x) σ i (ξ) = F(σ i (x)) convergence test ε i+1 (ξ) = ε i (ξ) -Γ 0 (ξ) : σ i (ξ) (10)
which start from the first term given by:

ε 1 (ξ) =        0 ∀ξ ̸ = 0 E for ξ = 0 (11)
In [START_REF] Helsing | An integral equation method for elastostatics of periodic composites[END_REF], Γ 0 (ξ) is the periodic "strain Green's tensor" , associated to the reference medium of rigidity C 0 . When an isotropic reference medium is considered, namely C 0 = 3k 0 J + 2µ 0 K with J = i ⊗ i/3 and K = I -J (i and I being the identity for second and fourth order tensors), the Green's tensor is given, for every ξ ̸ = 0, by:

Γ 0 (ξ) = 1 λ 0 + 2µ 0 E 2 (ξ) + 1 2µ 0 E 4 (ξ) (12)
with Γ 0 (ξ) = 0 when ξ = 0. In [START_REF] Idiart | Infinite-contrast periodic composites with strongly nonlinear behavior: Effective-medium theory versus full-field simulations[END_REF], F and F -1 denote the Fourier Transform and its inverse. In [START_REF] Helsing | An integral equation method for elastostatics of periodic composites[END_REF], the macroscopic strain, E, is prescribed whereas the macroscopic stress, Σ, is obtained at convergence of the iterative scheme [START_REF] Helsing | An integral equation method for elastostatics of periodic composites[END_REF] as follows: Σ = σ(ξ = 0). The strain field ε i (ξ) is compatible at each step, Q(ξ) : ε i (ξ) = 0, while the equilibrium for ε i (ξ) is verified only at convergence of the iterative scheme.

The following convergence test can then be used in [START_REF] Helsing | An integral equation method for elastostatics of periodic composites[END_REF]:

∥P(ξ) : σ i (ξ)∥ ∥σ i (ξ)∥ < ϵ ( 13
)
where ∥ • ∥ denotes the Frobenius norm. The value ϵ = 10 -3 has been used in this paper.

Dually, the solution of the inhomogeneity problem (1), can be computed by using the following iterative scheme:

                         σ i (x) = F -1 (σ i (ξ)) ε i (x) = S(x) : σ i (x) ε i (ξ) = F(ε i (x)) convergence test σ i+1 (ξ) = σ i (ξ) -∆ 0 (ξ) : ε i (ξ) (14)
which starts from the first term given by:

σ 1 (ξ) =        0 ∀ξ ̸ = 0 Σ for ξ = 0 (15)
The macroscopic stress Σ is prescribed over the unit cell. It means that the macroscopic strain E is now considered as an unknown and is obtained at convergence. In ( 14), ∆ 0 (ξ) is the "stress Green's tensor", given, for ξ ̸ = 0, by:

∆ 0 (ξ) = C 0 -C 0 : Γ 0 (ξ) : C 0 = 2µ 0 (3λ 0 + 2µ 0 ) λ 0 + 2µ 0 E 1 (ξ) + 2µ 0 E 3 (ξ) (16)
with ∆ 0 (ξ) = 0 when ξ = 0. In ( 14), the stress field σ i (ξ) in ( 14) is equilibrated at each step while the strain field ε i (ξ) is compatible only at convergence of the iterative scheme. The following convergence test is then used:

∥Q(ξ) : ε(ξ)∥ ∥ε(ξ)∥ < ϵ ( 17
)
In which the value ϵ = 10 -3 will be used in the following.

Estimation of the convergence of the basic schemes

The condition giving the convergence of the strain based iterative shceme has been derived in [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF]. Here we present a modification of that result which is specially interesting for voided composites. Next, this result will be extended to the stress based iterative scheme.

The convergence of the strain based iterative scheme

In [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF], the authors investigate the conditions giving the convergence of the Neumann series. It is well known that the Neumann series converges only if the eigenvalues of the operator are strictly comprised between -1 and 1. From another point of view, the original iterative scheme used by these authors can be written as a Neumann series in the operator

-Γ 0 (x) * (C(x) -C 0 ).
Studying the eigenvalues of this operator, they found that the convergence cannot be ensured when the composite material comprises voids, because this case produces eigenvalues of the operator equal to 1.

In the present paragraph we propose to establish a modification of that result.

For this, we focus our attention on the convergence of the iterative scheme (10) instead of the Neumann series. Let us therefore define by ϕ and by e ϕ respectively an eigenvalue and the associated eigenvector of the linear operator I -Γ 0 (x) * C(x). They comply with:

ϕe ϕ (x) = e ϕ (x) -Γ 0 (x) * C(x) : e ϕ (x) (18)
As shown in appendix A the convergence of such an iterative scheme is ensured if the eigenvalues comply with: -1 < ϕ ≤ 1. In [START_REF] Michel | Nonuniform transformation field analysis[END_REF], it can be observed that e ϕ (x) is compatible and has a null volume average if ϕ ̸ = 1. However, for every compatible field ε(x), such that < ε(x) > V = 0, one has:

ε(x) = Γ 0 (x) * C 0 : ε(x) (19)
Then, equation ( 18) can be also put into the form:

Γ 0 (x) * [ C(x) -(1 -ϕ)C 0 ] : e ϕ (x) = 0 (20)
Fhis relation (which is also obviously true for ϕ = 1) indicates that (C(x)-(1ϕ)C 0 ) : e ϕ (x) is an equilibrated field. As explained in [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF], e ϕ (x) is then the solution of an elasticity problem corresponding to an elastic heterogeneous material with stiffness

C ′ (x) = C(x) -(1 -ϕ)C 0 ,
and subjected to a zero macroscopic strain. Such a problem has the trivial solution e ϕ (x) = 0 if C ′ (x) is positive or negative definite. At this stage, the procedure is identical to that conducted in [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF] with the substantial modification that the value ϕ = 1 is now allowed. It leads to the following result (the reader will refer to the paper of Michel et al. [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF] for more details about the procedure):

k max < 2k 0 , k min ≥ 0 and µ max < 2µ 0 , µ min ≥ 0 (21)
with:

k min = min[k 1 , k 2 , ..., k N ], k max = max[k 1 , k 2 , ..., k N ] and µ min = min[µ 1 , µ 2 , ..., µ N ], µ max = max[µ 1 , µ 2 , ..., µ N ].
It is worthwhile mentioning that the value k min = 0 and µ min = 0 are acceptable contrarily to the result of [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF].

In fact, the very reason of such a difference with the result of [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF] is that two iterative schemes can be used : the one corresponding to [START_REF] Helsing | An integral equation method for elastostatics of periodic composites[END_REF] and an alternative scheme leading to a Neumann series in -Γ 0 (x) * (C(x) -C 0 ). Studying the spectra of the operator used in the alternative scheme seems to lead to no convergence in the case of void-containing materials. However, it has been proved by [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] that both schemes are strictly equivalent. As a consequence, the previous proof that scheme [START_REF] Helsing | An integral equation method for elastostatics of periodic composites[END_REF] converges even for void containing materials proves that the the alternative scheme converges also in this case, as show also numerous numerical tests.

The convergence of the stress based iterative scheme

Our purpose is now to establish a similar result in the case of the stress based iterative scheme [START_REF] Lee | Modeling the Viscoplastic Micromechanical Response of Two-Phase Materials using Fast Fourier Transforms[END_REF]. To this aim, the eigenvalues of the linear operator I -∆ 0 (x) * S(x) are now studied. As effected in the previous section, we introduce an eigenvalue ϕ and the associated eigenvector e ϕ such that:

ϕe ϕ (x) = e ϕ (x) -∆ 0 (x) * (S(x) : e ϕ (x)) (22)
If ϕ ̸ = 1, the eigenvector complies with equilibrium and has a null volume average. It follows that e ϕ complies with:

e ϕ (x) = ∆ 0 (x) * S 0 : e ϕ (x) (23)
Consequently, ( 22) can be written as:

∆ 0 (x) * [ (S(x) -(1 -ϕ)S 0 ) : e ϕ (x) ] = 0 (24)
Similarly to the case of the strain based iterative scheme, one observes that (S(x) -(1 -ϕ)S 0 ) : e ϕ (x) is a compatible field. So, e ϕ is the solution of an elasticity problem corresponding to an elastic heterogeneous material with the compliance S ′ (x) = (S(x) -(1 -ϕ)S 0 ) : e ϕ (x) and subjected to a zero macroscopic stress. There exists a non trivial solution for e ϕ (x) if S ′ (x) is not strictly positive or negative definite. Again, applying the procedure described in [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF], we obtain:

1 k 0 > 1 2k min , 1 2k max ≥ 0 and 1 µ 0 > 1 2µ min , 1 2µ max ≥ 0 (25)
Which shows that the case of rigid inclusion corresponding to the case of k max or µ max being infinite corresponds to admissible bounds for the inverses of these moduli. In the next subsection, we propose to present a few numerical examples to show the features of the convergence, more specially for infinite contrasts, and showing that improvements are necessary if one wants to deal with microstructures containing simultaneously voids and rigid inclusions.

Numerical study of the rate of convergence of the basic schemes

The rate of convergence of the basic schemes is examined through numerical calculations in the case of composites reinforced by long fibers arranged along a hexagonal array (the unit cell being defined in figure 1). The unit cell is a rectangle and contains an inclusion centered at x = 0 and four quarters of inclusions located at the corner of the cell. In our computations, the radii of the fibers are chosen as R = 0.4 (the area of the unit cell being taken as S = 1) and N = 64 wave vectors are taken along each space direction. Both phases, matrix and inclusions, are assumed to be incompressible and we denote the "contrast" c as c = µ I /µ M . Plane strain conditions are assumed and the unit cell is subjected to the macroscopic strain E 11 = -E 22 = 1 when the strain based method is used. For the illustrations of the stress based iterative scheme, the macroscopic stress Σ 11 = -Σ 22 = 1 is applied over the unit cell. The exact Fourier transform is thereafter replaced by the discrete Fourier transform which is computed by means of the FFT algorithm. In [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF], the stress-strain relation is applied in the real space at each point along a regular grid. For many inclusion-matrix problems, such a description is only approximated. As an example, if cylindrical fibers are concerned, the cross section of the fiber must be described within the real space by using a regular mesh, which approximates the form of the inclusion. In [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF], exact expressions for the Fourier transforms of the characteristic functions I α (x) have been used. These Fourier transforms are named "shape coefficients" and are the basis of the approximation for the elastic properties of elastic media obtained by Nemat-Nasser and Hori [START_REF] Nemat-Nasser | Micromechanics: Overall Properties of Heterogeneous Materials[END_REF]. These "shape coefficients" are defined by:

I α (ξ) = ∫ Vα exp(-iξ.x)dx (26)
The use of shape functions improves the convergence of the iterative scheme.

On figure 2, we represent the logarithm of the number of iterations (log(N )), needed for obtaining the convergence of the strain based iterative scheme [START_REF] Helsing | An integral equation method for elastostatics of periodic composites[END_REF], as a function of the logarithm of the "contrast" (log(c)). On this figure, we compare the original approach [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] (full line with points), and the method which uses the shape functions (full line with circles). The results show that the con- 2007) [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] Fig. 2. Number of iterations at convergence as a function of the "contrast". Comparison between the initial strain based iterative scheme [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] and the method using shape functions [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF].

vergence is obtained when c = 0 for both approaches, with or without using the shape coefficients. However, one observes that the use of the shape functions improves the convergence of the basic scheme, and this, whatever the value of the contrast. When the value of c = µ I /µ M is higher than 1, it can be observed that the number of iterations linearly increases with the contrast (in the loglog frame). The strain based iterative scheme diverges therefore for an infinite "contrast". The case of a composite with rigid inclusions cannot be handled by the method. Note that the shear modulus of the reference medium used in [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] and also used here is defined by: µ 0 = (µ M + µ I )/2. With that choice, the convergence of the iterative scheme is found to be obtained with the lowest number of iterations. When the shape coefficients are used, it has been found that the better choice for µ 0 does not coincide exactly with µ 0 = (µ M + µ I )/2. The reference shear modulus is chosen on the form: µ 0 = (1 -ω)µ M + ωµ I where the value ω = 0.5 is taken for c > 1 but the value ω = 0.45 is taken for c < 1. Note that the convergence is still achieved when the value ω = 0.5 is used and for a contrast c < 1, however the choice ω = 0.45 instead of ω = 0.5 leads to a gain of convergence rate. On Figure 3 we represent the logarithm of the number of iterations at convergence as a function of log(c) = log(µ I /µ M )for the strain based scheme (full line with points), and for the stress based scheme (full line with circles). Both iterative schemes are computed with a representation of the elastic tensors using the shape functions. It can be observed that the stress based approach converges rapidly in the domain of high values of the contrast c. However, for soft inclusions, namely in the region defined by c < 1, it can be seen that the number of iterations increases linearly when the logarithm of c decreases. It means that the case of materials containing voids cannot be treated by the stress based method. The optimal reference medium, leading to the lowest number of iterations, has been found to be defined by µ 0 = ((1 -ω)/µ M + ω/µ I ) -1 . The value ω = 0.5 is considered when c < 1 while ω = 0.45 when c > 1. To summarize, the strain based approach is better suited for composites containing soft inclusions, namely for which µ I < µ m . However, when the composite is made up of stiff inclusions (µ I > µ M ), the stress based iterative scheme is better suited. Nevertheless, it suggests that the basic schemes are not adapted for computing the effective properties of composites which contain both stiff and soft inclusions. Moreover the case of composites with rigid inclusions and voids cannot be handled by any of these methods. In the next section, we propose therefore a new iterative scheme which improves the convergence, by comparison with the basic schemes.

A new iterative scheme based on the polarization

Description of the algorithm

Let us introduce in the system of equations (1) the "polarization tensor"τ (x) = (C(x) -C 0 ) : ε(x):

                                             ε(x) = 1 2 (∇u(x) + ∇ t u(x)), ∀x ∈ V div(σ(x)) = 0, ∀x ∈ V ε(x) = (C(x) -C 0 ) -1 : τ (x), ∀x ∈ V σ(x) = τ (x) + C 0 : ε(x), ∀x ∈ V u(x) -E.x periodic σ(x).n antiperiodic < τ (x) > V = T (27)
In this problem, the uniform polarization T is now prescribed over the unit cell. It follows that the macroscopic strain E and stress Σ are now both considered as unknowns. However, due to the linearity of the equations, the relations giving the macroscopic stress and strain fields as functions of the uniform polarization T can be expressed as:

E = ε(ξ = 0) = A : T , Σ = σ(ξ = 0) = B : T ( 28
)
Note that Σ = T + C 0 : E. As a consequence, it follows that the fourth order tensors A and B are linked by B = I + C 0 : A. The effective elastic tensor C hom and the effective compliance S hom ), are then defined by:

C hom = A -1 + C 0 or S hom = S 0 : (I -B -1 ) (29)
The solution of ( 27) can be computed by using the following iterative scheme:

                               τ i (x) = F -1 (τ i (ξ)) ε i (x) = (C(x) -C 0 ) -1 : τ i (x) ε i (ξ) = F(ε i (x)) σ i (ξ) = C 0 : ε i (ξ) + τ i (ξ) convergence test τ i+1 (ξ) = τ i (ξ) -αC 0 : Γ 0 (ξ) : σ i (ξ) -β∆ 0 (ξ) : ε i (ξ) (30)
where α and β are coefficients which will be chosen thereafter in order to obtain the best rate of convergence. In (30), tensors Γ 0 (ξ) and ∆ 0 (ξ) are respectively the periodic Green's tensors for strain and stress defined by equations ( 12) and [START_REF] Michel | A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast[END_REF]. In this iterative scheme, the stress field σ i (ξ) and the strain field ε i (x) are respectively equilibrated and compatible fields only at convergence of the iterative scheme. Indeed, the convergence is achieved when τ i+1 (ξ) = τ i (ξ) in (30) and then, when:

αΓ 0 (ξ) : σ i (ξ) + βS 0 : ∆ 0 (ξ) : ε i (ξ) = 0 (31)
However, the projection of equation (31) along Q(ξ) and P(ξ) : C 0 provides Q(ξ) : ε i (ξ) and P(ξ) : σ i (ξ) respectively. Then, it means that ε i (ξ) and σ i (ξ) are compatible and equilibrated respectively. The iterative scheme is stopped when:

max ( ∥Q(ξ) : ε i (ξ)∥ ∥ε i (ξ)∥ , ∥P(ξ) : σ i (ξ)∥ ∥σ i (ξ)∥ ) < ϵ (32)
and the first term of the series is defined by:

τ 1 (ξ) =        0 ∀ξ ̸ = 0 T for ξ = 0 (33)
It is worthwhile mentioning that this polarization based iterative scheme has some similarities with the accelerated scheme used in [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF]. More precisely, it is shown in appendix A that the accelerated scheme coincides with (30) for α = 2 and T = (< C(x) > V -C 0 ) : E. Note that, with that choice, the prescribed polarization T diverges for rigid inclusions while, in our approach, the uniform polarization T is finite, being prescribed over the unit cell, its value being chosen independently of the elastic properties of the heterogeneous medium. In addition, a specific feature of the polarization based scheme is that it involves always finite linear relations between local stress, strain and polarization fields (second and fourth lines of (30)) even if the elasticity and compliance tensors are not finite. Note also that a dual form, of the iterative scheme (30), can also be used. This iterative scheme is based on the eigenstrain η and reads:

                               η i (x) = F -1 (η i (ξ)) σ i (x) = (S(x) -S 0 ) -1 : η i (x) σ i (ξ) = F(σ i (x)) ε i (ξ) = η i (ξ) + S 0 : σ i (ξ) convergence test η i+1 (ξ) = η i (ξ) -αS 0 : ∆ 0 (ξ) : ε i (ξ) -βΓ 0 (ξ) : σ i (ξ) (34)
In (34), the stress is computed in the real space and uses the compliance tensor S(x) whereas in (30) the rigidity C(x) is used. The convergence test is still given by (31). The first term of the series is associated to a uniform eigenstrain, denoted N , prescribed over the unit cell. Using this "eigen-strainbased" scheme, it has been found that its rate of convergence is comparable to the one of the polarization scheme. For this reason, only results for the polarization scheme are presented in the following.

The convergence of the polarization-based iterative scheme

In the iterative scheme (30), the values of coefficients α, β and of the elastic moduli of the reference medium, λ 0 , µ 0 , will be chosen in order to obtain the convergence. In other words, they must be chosen in order to ensure that eigenvalues of the operator

I -αC 0 : Γ 0 (x) * (C(x) : δC -1 (x)) -β∆ 0 (x) * δC -1 (x) (35) are comprised in the interval ] -1, 1]. In (35), we have introduced δC -1 (x) = (C(x) -C 0 ) -1 .
Let us introduce an eigenvalue ϕ of that operator and its associated eigenvector e ϕ (x).They comply with:

ϕe ϕ (x) = e ϕ (x) -αC 0 : Γ 0 (x) * (C(x) : δC -1 (x) : e ϕ (x)) +α∆ 0 (x) * (δC -1 (x) : e ϕ (x)) (36)
If ϕ ̸ = 0, one has < e ϕ (x) >= 0 and for any such second order tensor:

e ϕ (x) = (∆ 0 (x) : S 0 + C 0 : Γ 0 (x)) * e ϕ (x) (37)
It follows that (36) can be put into the form:

C 0 : Γ 0 (x) * σ ϕ (x) + ∆ 0 (x) * ε ϕ (x) = 0 ( 38
)
with:

σ ϕ (x) = αC(x) : δC -1 (x) : e ϕ (x) -(1 -ϕ)e ϕ (x) ε ϕ (x) = βδC -1 (x) : e ϕ (x) -(1 -ϕ)S 0 : e ϕ (x) (39)
Owing to property (38) and following an argument similar to the one used previously for (31), ε ϕ (x) and σ ϕ (x) are compatible and equilibrated fields. They are solutions of an elasticity problem related to a heterogeneous material having a stiffness C ′ (x) given by:

C ′ (x) = [ αC(x) : δC -1 (x) -(1 -ϕ)I ] : [ βδC -1 (x) -(1 -ϕ)S 0 ] -1 (40)
and subjected to a zero polarization. That problem does not lead to the trivial solution e ϕ (x) = 0, if C ′ (x) is not strictly positive or negative definite. Let us recall that the convergence of the iterative scheme is ensured if

-1 < ϕ ≤ 1.
So, a sufficient condition for having the convergence is that C ′ (x) must be strictly positive or negative definite only within the interval ] -∞, -1] and ]1, +∞[. The elastic properties of each phase and the reference modulus being isotropic, C ′ (x) is also isotropic and can be put into the form

C ′ (x) = 3k ′ (x)J + 2µ ′ (x)K,
where k ′ (x) and µ ′ (x) are given by:

k ′ (x) = k 0 αk(x) -(1 -ϕ)(k(x) -k 0 ) βk 0 -(1 -ϕ)(k(x) -k 0 ) µ ′ (x) = µ 0 αµ(x) -(1 -ϕ)(µ(x) -µ 0 ) βµ 0 -(1 -ϕ)(µ(x) -µ 0 ) (41)
So, a sufficient condition for having the convergence is that k ′ (x) and µ ′ (x) must be strictly positive or negative definite only within the interval ]-∞, -1] and ]1, +∞[. Functions k ′ (x) and µ ′ (x) are zero if:

ϕ = 1 - αk(x) k(x) -k 0 , ϕ = 1 - αµ(x) µ(x) -µ 0 (42)
and are singular if:

ϕ = 1 - βk 0 k(x) -k 0 , ϕ = 1 - βµ 0 µ(x) -µ 0 ( 43
)
It can be easily verified that, if α, β, k 0 , µ 0 comply with:

0 ≤ α < 2, -2 < β ≤ 0, -∞ < µ 0 < 0, -∞ < k 0 < 0 (44)
the values of ϕ given by ( 42) and (43) are comprised in the interval ] -1, 1] and then k ′ (x) and µ ′ (x) are strictly positive or negative definite within the interval ] -∞, -1] and ]1, +∞[. The inequalities in (44) give then a sufficient condition for having the convergence of the polarization based iterative scheme. It is worth noting that the conditions of convergence of the basic schemes (see [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] and ( 25)) depend on the properties of the elastic moduli of the heterogeneous medium while the conditions (44) do not depend anymore on these properties. Thus, it is expected that the convergence of the polarization based iterative scheme will be reached independently of the values of the contrast. This is investigated in the next section through several numerical applications. In addition, it is also interesting to notice that the elastic properties of the reference medium are negative, a feature which is common with the "accelerated" scheme [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF] 4.3 Numerical investigation of the rate of convergence

Application to a two phase composite

We now investigate the rate of convergence of the polarization based iterative scheme through some numerical examples. As a first illustration of the method, we still consider the example of the section 3.3. The composite is made up of long fibers which are periodically arranged along a hexagonal array (see figure 1). The radius of the fibers is R = 0.4 (where it is recalled that the area of the unit cell is S = 1) and N = 64 wave vectors are used in each direction. The matrix and the inclusions are assumed to be incompressible and we denote by c = µ I /µ M the contrast. Plane strain conditions are considered and a prescribed uniform polarization is defined by its components: T 11 = -T 22 = 1, other ones being null. Du to the similarity between the conditions on α and -β for the conditions ensuring convergence of the polarization scheme, it is set β = -α. For an incompressible composite, the conditions in (44) reduce to µ 0 < 0 and 0 ≤ α < 2. We propose to validate that condition with numerical calculations. For this, we consider fixed values of α within the interval ]0, 2[ and we measure the rate of convergence of the iterative scheme as a function of the choice of µ 0 . For instance, on figure 4 we represent the number of iterations at convergence of the iterative scheme as a function of the shear modulus of the reference medium, µ 0 , for a fixed α = 3/2 and for the two limit cases of rigid inclusions and voids. It can be observed that the convergence is attained whatever the value of µ 0 < 0 (excepted for µ 0 = -∞). It is also observed that each curve has a minimum which defines the optimal reference medium. At a distance of this minimum, the logarithm of the number of iterations increases linearly with the logarithm of µ 0 and it suggests that the iterative scheme diverges for µ 0 = 0 and µ 0 = -∞.

The determination of the optimal reference medium is of a great importance but is a very difficult task since it depends on several parameters such as the contrast, the resolution, the geometry and the size of the inclusions. The opti-- mal reference medium is very sensitive to the contrast when the basic schemes (strain based or stress based) are considered. For example, the optimal reference medium defined for the strain based iterative scheme depends linearly on the contrast and tends to infinity in the case of rigid inclusion. In the present work, the polarization scheme uses additional parameters, namely α and β, which have to be optimized in order to obtain the lowest number of iterations at convergence. In the absence of theoretical arguments yielding to the optimal choice of (α, β, µ 0 , k 0 ), they must be determined numerically. For example, in figure 5 are reported the interval for µ 0 /µ M (µ M being the shear modulus of the matrix) corresponding to the optimal value of µ 0 for several values of the contrast varying from 10 -6 to 10 6 and for α = -β = 3/2. It can be observed that, contrarily to the basic approaches, the optimal reference medium of the polarization based method does not diverge in the limited cases of voided and rigid inclusions. When other values of the coefficients α are considered within the interval ]0, 2[, the results are similar to the ones given on figures 4 and 5. However, the choice α = 3/2 has been found to lead to the better rate of convergence of the iterative scheme. So, in the absence of any other arguments, the choice α = 3/2 is kept for the other examples treated in the following. We now compare the new iterative scheme with the basic ones. In Figure 6 we represent the number of iterations at convergence as a function of the contrast (in a log-log frame), the optimal value of µ 0 being determined numerically (see figure 5). On the same figure the results obtained with the strain and the stress based approaches were also reported . The polarization based approach clearly appears to be more efficient since the convergence is reached after a lower number of iterations and this, whatever the value of the contrast. It has been found that the convergence is obtained after 12 iterations when log(µ I /µ M ) = ±6. In figure 7, is plotted the prediction for the effective shear modulus as a function of the number of wave vectors (used for computing the Fourier transform and its inverse). The contrast has been chosen as µ I /µ M = 10. It can seen that, as observed currently by users of basic schemes, the strain based iterative scheme leads to an upper type estimate of the homogenized elastic tensor whereas the stress based approach leads to a lower type estimate of that modulus. By comparison, the prediction obtained with the polarization based method is always comprised between the estimates of the basic methods. As a consequence, the number of wave vectors needed for ensuring the convergence is lower.

Application to three phases composites

To check the ability of the polarization based method to deal with highly contrasted microstructures, we investigate the rate of convergence for three-phases composites which contain both soft and stiff constituents. Two examples are considered and are represented by their unit cell on figure 8. Those microstructures are depicted below: • Microstructure a: the composite is made up of an incompressible elastic matrix whose shear modulus is µ M which contains circular inclusions being arbitrary distributed within the unit cell. Two kinds of inclusions, with the same proportion, are considered: (i) circular inclusions with the shear modulus µ I1 = pµ M , (ii) circular inclusions with modulus µ I2 = µ M /p (the inclusions are also considered incompressible). The unit cell is squared (with the area S = 1) and the inclusions have the same radius R = 0.05. The volume fraction of stiff and soft inclusions is f I1 = f I2 = 0.0785. The contrast between the phases is controlled by the parameter p which varies from 1 to 10 6 in our applications. Note that the particular case p = 1 corresponds to a homogeneous material whereas the limited case p = +∞ corresponds to a composite which contains both cavities and rigid inclusions. • Microstructure b: the composite is made of an incompressible elastic matrix whose shear modulus is µ M and containing circular inclusions having a circular interphase. A squared unit cell is still considered (with the area S = 1), the radius of the inclusions (the cores) is R = 0.05 and the external radius of the interphase (the circular shells) is R = 0.075. The volume fraction of the inclusions is f I1 = 0.118 and the volume fraction of the interphase is f I2 = 0.147. The inclusions and the interphase are both chosen incompressible. The shear modulus of the inclusions is µ I1 = pµ M whereas the shear modulus of the interphase is µ I2 = µ M /p. The contrast is again controlled by the parameter p whose value is taken from 10 -6 to 10 6 .

In our calculations a grid 256 × 256 is considered, the macroscopic loading is the same as in the last example. The number of iterations at convergence is plotted in figure 9 as a function of the parameter p, in the case of the microstructure "a", constituted of soft and stiff inclusions. The logarithm of the number of iterations increases linearly with log(p) for the two basic schemes. As expected, the case of a composite made up of both voids and rigid inclusion (corresponding to p = +∞) cannot be handled by these methods. When the polarization based iterative scheme is used, the number of iteration increases with p but reaches moderate values. Indeed, for higher values of p, the number of iterations converges at N = 19 iterations. In figure 10 is plotted the number of iterations at convergence as a function of p for the second microstructure involving interphases surrounding circular inclusions (microstructure "b" on figure 8). When the basic schemes are used, it is observed again that the computation times become very long when p is large (corresponding to a very soft interphase surrounding a stiff inclusion) and when p is low (the case of a soft inclusion with a very stiff interphase). Again, it is observed that the convergence is attained with the polarization based method and this for any value of the parameter p. It can be also noticed that the number of iterations for having the convergence, for this kind of microstructure, is higher than for the first microstructure (microstructure a on figure 8). When log(p) = -6 the number of iterations at convergence is N = 52 whereas N = 57 iterations are needed when log(p) = 6. This diminution of the rate of convergence can be attributed to the fact that the soft and the stiff phases are connected for the microstructure "b", leading to high local fluctuations of the mechanical fields near the interface between core and shell. Some indications on the numerical parameters which were used in the iterative schemes are of interest. When the strain based iterative scheme is used, the reference modulus given in [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF] has been considered: µ 0 = (µ min + µ max )/2 = µ M (p + 1/p)/2. When the stress based iterative scheme is used, the reference modulus is chosen as 1/µ 0 = (1/µ min + 1/µ max )/2 = (p + 1/p)/(2µ M ). When the polarization based iterative scheme is used, the optimal reference medium has been determined numerically. In the case of the composite made up of stiff and soft inclusions (microstructure "a") it has been found that the optimal shear modulus is µ 0 = -µ m and then is independent of the parameter p.

In the case of the composite with circular interphases (microstructure "b") it has been found that µ 0 = -µ M does not correspond to the optimal reference medium (squares on figure [START_REF] Helsing | An integral equation method for elastostatics of periodic composites[END_REF]). In this case, the use of an optimized reference medium (stars on figure [START_REF] Helsing | An integral equation method for elastostatics of periodic composites[END_REF]) leads to a gain on the rate of convergence. It suggests that the optimal reference not only depends on the elastic properties of constituents but also on their spatial distribution within the unit cell. The question relative to the determination of the optimal reference medium needs to be further investigated in order to be able to predict the relation between the microstructure ant this optimum.

Conclusion

A new iterative scheme has been proposed for the computation of the macroscopic properties of elastic composites. The main features and advantages of the present method can be summarized as follows:

• The method uses the same ingredients as the basic schemes: the periodic Green tensors for the strain and the stress, the FFT algorithm for computing the Fourier Transform and its inverse... No additional treatment or numerical tools (Lagrangian, conjugate gradient,...) are required for its implementation. • The approach uses an iterative scheme for solving an elasticity problem expressed in terms of the polarization. In that problem a uniform polarization is prescribed over the unit cell instead of a uniform macroscopic strain or stress in the classical approaches. • It has been shown that the condition ensuring the convergence of this new iterative scheme is independent of the elastic moduli of the heterogeneous material. This result has been proven theoretically and tested numerically through several examples. More specifically, it has been shown that the problem of composites containing both voids and rigid inclusions can be solved by the method with a very good rate of convergence. The application of the method for determining the overall properties of composites with soft or stiff interphases has been addressed in the next part of the paper. The method can be extended to the context of non linear composites; this point will be the subject of a forthcoming work.

B Connection with the accelerated scheme

The recurrence relation of the accelerated scheme is [START_REF] Michel | A computational scheme for linear and nonlinear composites with arbitrary phase contrast[END_REF]:

ε i+1 (x) = ε i (x)

Fig. 1 .

 1 Fig. 1. Unit cell of the composite reinforced by long fibers arranged along a hexagonal array.

Fig. 3 .

 3 Fig. 3. Number of iterations at convergence as a function of the "contrast": comparison between the strain based and the stress based formulations.

Fig. 4 .

 4 Fig. 4. Number of iterations at convergence as a function of the ratio µ 0 /µ M for α = -β = 3/2.

Fig. 5 .

 5 Fig. 5. Value of the optimal reference medium as a function of the logarithm of the contrast

20

 20

Fig. 6 .

 6 Fig.6. Number of iterations at convergence as a function of the contrast: comparison between the strain, the stress and the polarization based approaches.

Fig. 7 .

 7 Fig. 7. The effective shear modulus for a contrast c = 10. Comparison of the estimations obtained from the strain, the stress and the polarization approaches.

Fig. 8 .

 8 Fig. 8. The unit cell of the three phase composite. Case a: composite with soft and stiff circular inclusions. Case b: composite reinforced with circular inclusions surrounded by circular interphases.

Fig. 9 .

 9 Fig. 9. Number of iterations as a function of the ratio p = µ I1 /µ M = µ M /µ I2 in the case of three phases composite constituted of both stiff and soft circular inclusions.

Fig. 10 .

 10 Fig. 10. Number of iterations as a function of the ratio p = µ I1 /µ M = µ M /µ I2 in the case of coated circular inclusions.

τ

  -2(C(x) -C 0 ) -1 : C 0 : [ Γ 0 (x) * (C 0 : ε i (x)) -ε i (x) + E + Γ 0 * (C(x) : ε i (x)) ] (B.1)which is initialized by ε 1 (x) = E. We multiply the terms at the right and at the left of the equality by (C(x)-C 0 ). It leads to:τ i+1 (x) = τ i (x) -2C 0 : [ Γ 0 (x) * (C 0 : ε i (x)) -ε i (x) + E +Γ 0 (x) * (C(x) : ε i (x)) ] (B.2)Using now the following equality:C 0 : [ Γ 0 (x) * (C 0 : ε i (x)) -ε i (x) + E ] = -∆ 0 (x) * ε i (x) (B.3)and introducing σ i (x) = C(x) : ε i (x), it gives:τ i+1 (x) = τ i (x) -2C 0 : Γ 0 (x) * σ i (x) + 2∆ 0 (x) * ε i (x) (B.4)Which reads, in the Fourier space: i+1 (ξ) = τ i (ξ) -2C 0 : Γ 0 (x) :σ i (ξ) + 2∆ 0 (ξ) : ε i (ξ) τ i+1 (ξ) = T (B.5)In the recurrence relation (B.1), the macroscopic strain is prescribed over the unit cell. It follows that T =< τ (x) > V = (< C(x) > V -C 0 ) : E.

A Condition of convergence of a recurrence relation

Let us consider a recurrence relation on the form:

and starting form a 1 = A. denoting by ϕ k and e ϕ k for k = 1..K the eigenvalues and the associated eigenvectors of the linear operator L. It follows that the iterate i + 1, one has:

where the A k are such that:

Three case are distinguished:

the factor (ϕ k ) i increases at each iteration and the iterative scheme diverges if

and the component of a i along e ϕ k remains unchanged and equal to

the factor (ϕ k ) i decreases at each iteration and tends to zero when i → ∞.

Consequently, the convergence is ensured if -1 < ϕ k ≤ 1. Moreover, assuming that the K ′ first eigenvalues are equal to 1 while -1 < ϕ k < 1 for k = K ′ + 1..K, when the limit i → +∞ is taken, one has:

Which shows that a ∞ = 0 if all the eigenvalues are strictly comprised between 1 and -1. There must exist some eigenvalues equal to 1 if one hopes that the iterative scheme does not converge to a ∞ = 0.


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

Guy Bonnet
email: guy.bonnet@univ-paris-est.fr

Guy Lauriat
email: guy.lauriat@univ-paris-est.fr

A FFT-based method to compute the permeability induced by a Stokes slip flow through a porous medium

Keywords: Porous media, homogenization, effective properties, permeability, viscous fluid Milieux poreux, homogénéisation, propriétés effectives, perméabilité, fluide visqueux

The paper presents a FFT based-method for obtaining the permeability of a periodic micro-porous medium. The periodic medium is constituted of a rigid solid matrix saturated by a viscous fluid. The flow obeys the Stokes equations and a slip condition at the surface of the rigid skeleton is considered. The permeability is obtained from the homogenization of periodic media and an extension of the FFT method used for composite elastic media. The incorporation of the slip condition is made trough the introduction of an interphase between the fluid and the solid. The method of solution uses the classical expansion along Neumann series of the velocity field of the periodic problem and Green's tensor in Fourier space.

To cite this article :

Résumé

Cet article présente une méthode basée sur la FFT pour la détermination de la perméabilité d'un milieu poreux périodique. Le milieu poreux est constitué d'une matrice solide indéformable saturée par un fluide visqueux. L'écoulement obéit aux équations de Stokes et une condition de glissement sur la surface du squelette solide est également considérée. La perméabilité est alors obtenue en résolvant les équations issues de l'homogénéisation de milieux périodiques et par une extension de la méthode FFT utilisée pour les matériaux composites élastiques périodiques. La prise en compte de la condition de glissement est effectuée en introduisant une interphase entre le fluide visqueux et le solide. La méthode repose sur une représentation du champ de vitesses en séries de Fourier et sur l'expression explicite du tenseur de Green dans l'espace de Fourier. Pour citer cet

Introduction

The determination of the effective properties of porous media, and particularly the permeability, is of fundamental importance in several practical problems in geomechanics, biomechanics, among others, and it has been the subject of many studies. Fluid flow through porous media is usually modeled using a volume-averaged approach and Darcy's law is valid for fluid flows with pore Reynolds numbers less than one. For forced flows, Darcy's law reads:

Q i = - k ij µ ∇ j P ( 1
)
where Q i and ∇ i P are respectively the components of the flux and of the pressure gradient. µ is the dynamic viscosity and k ij the components of the permeability tensor of the considered medium. The determination of the permeability of a periodic porous medium knowing the geometry of the microstructure has been studied by many authors, see for instance [START_REF] Sangani | Slow flow past periodic arrays of cylinders with application to heat transfer[END_REF], [START_REF] Barrère | Détermination numérique de la perméabilité en milieu poreux périodique tridimensionnel[END_REF], [START_REF] Wang | Stokes flow through a rectangular array of circular cylinders[END_REF], [START_REF] Kadaksham | Permeability of periodic arrays of spheres[END_REF]. The methods to obtain a solution of the problem at the microscopic scale include expansion into series of eigenfunctions, finite element methods, etc. However, in many instances, as in the case of rarefied gases when the Knudsen number is small, a slip may occur at the surface of the skeleton. A semi-numerical method for the determination of the permeability in the case of a slip of flow the velocity field has been proposed by [START_REF] Wang | Stokes slip flow through square and triangular arrays of circular cylinders[END_REF] but such a method is limited to specific inclusions. The aim of the present paper is to propose a FFT-based method to compute the permeability of a periodic porous medium including the slip effect. The Fourier based-method has been firstly introduced by Moulinec and Suquet [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear properties of composites[END_REF] for the determination of the effective properties of elastic composites on the basis of previous works [START_REF] Brown | Solid mixture permittivities[END_REF][START_REF] Kroner | Statistical Continuum Mechanics[END_REF]. However the case of infinite contrast such as the case of rigid inclusions cannot be handled easily by the first method and so, the authors have proposed an alternate scheme based on augmented Lagrangian method [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF]. The case of rigid inclusions can also be solved by using the dual formulation [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF] which is extended to the case of porous media in the present paper. The consideration of a slip condition introduces a discontinuity in the velocity field at the surface of the rigid inclusions. Following [START_REF] Benveniste | Imperfect soft and stiff interfaces in twodimensional elasticity[END_REF], this discontinuity is introduced by considering an interphase of lower viscosity, between the fluid and the rigid inclusions. In fine, it is shown that the obtained iterative scheme takes a form similar to the one used for linear composites. The differences lie in the consideration of the first term introduced in this iterative scheme and in the recovering of the velocity field in our case, while elasticity problems need only deformation and stress fields. The methodology used to build this first term is depicted in section 4 and consists in introducing body forces in the rigid inclusion. An application to the case of aligned rigid cylinders is proposed in the last part of the paper.

Statement of the problem

Following [START_REF] Ene | Equations and surface phenomena for the flow in a porous medium[END_REF][START_REF] Auriault | Study of macroscopic behavior of a deformable porous medium[END_REF][START_REF] Levy | Fluid flow through an array of fixed particles[END_REF][START_REF] Sanchez-Palencia | Nonhomogeneous Media and Vibration Theory[END_REF] and [START_REF] Allaire | Homogenization of the Navier-Stokes equations with a slip boundary condition[END_REF][START_REF] Cioranescu | Homogenization of the Stokes problem with nonhomogeneous slip boundary condition[END_REF] for slip flow, we consider a periodic medium made up of a rigid skeleton and of an incompressible viscous fluid obeying the Stokes equations with a slip condition:

     µ f ∆v(x) + ∇p(x) = 0, div v * (x) = 0, ∀x ∈ V f v * n (x) = 0, v * t (x) = c∇ n v t (x) ∀x ∈ Γ ( 2
)
where c is the "slip coefficient" function of the nature of the fluid and the solid, v * (x) and p(x) are the local velocity field and pressure, respectively. Γ defines the surface of the solid and n is the normal unit vector taken on the surface Γ. being periodic. At the local scale, the components of the velocity field v * i (x) linearly depends on the components of the macroscopic pressure gradient :

v * i (x) = A ij (x)G j ( 3
)
where A ij (x) depends on the position vector, x. As shown in the framework of the homogenization applied to porous media [START_REF] Auriault | Study of macroscopic behavior of a deformable porous medium[END_REF], the average value over the unit cell of the velocity field reads:

< v * i (x) > V = 1 V ∫ V v * i (x)dV =< A ij (x) > V G j = - k ij µ G j ( 4
)
where k ij are the components of the permeability tensor. So, the computation of the permeability requires to determine the local velocity v * (x). That can be done by using a FFT based-method. However, the slip condition implies to introduce a discontinuity in the velocity field at the surface Γ of the inclusion.

To take into account the slip condition, we introduce an interphase located between the rigid inclusions and the viscous fluid (see figure 1). This interphase is defined by a layer which has a constant thickness denoted e. The dynamic viscosity of this layer is denoted µ e . The local problem ( 2) is replaced by :

µ(x)∆v * (x) + ∇p(x) = 0, div v * (x) = 0, v * (x) = 0 ∀ x ∈ Γ (5)
where the velocity field is null on the the surface Γ, µ(x) = µ f in the fluid and µ(x) = µ e in the interphase. The dynamic viscosity µ e is chosen on the form µ e = (eµ f )/c which allows to introduce a slip on the form given by ( 2) when e/a → 0, a being a representative lengthscale of the inclusions (see [START_REF] Benveniste | Imperfect soft and stiff interfaces in twodimensional elasticity[END_REF] for more details about the method).

The stress formulation of the problem and its solution

In this section, a "stress based" formulation of the problem is depicted. Its solution is given in the Fourier space. We suppose first that the interphase has a finite thickness and the limited case, e → 0 is examined secondly. In the fluid, the local stress can be put into the form σ = Σ + σ * where σ * ij is a periodic stress field and Σ ij is given in terms of the macroscopic pressure gradient, Σ ij = G k x k δ ij . In the rigid solid, the stress is undefined. However, we propose to make a continuation by continuity of σ * and Σ within the solid. We also introduce:

M(x) =              M f in the fluid M s = 0 in the solid M e in the interphase with: M α = 1 2µ α K (α = f, s, e) (6)
where µ s → +∞ within the solid, K = I-J where I is the fourth order identity tensor and J ijkl = δ ij δ kl /3. The local problem (5) can thus be cast under the following form:

div(σ * ) + f (x) = 0, ε * (x) = M(x) : σ * (x), ε * (x) = grad s (v * ) (7)
where ε * are the components of the strain rate and f i = Σ ij,j . Introducing a reference medium defined by the dynamic viscosity µ 0 and the fourth-order tensor M 0 = K/(2µ 0 ) (see [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF] for more details) it is easy to show that the periodic stress field complies with :

σ * (ξ) = -Ω(ξ).f (ξ) -∆ 0 (ξ) : [ ε * (ξ) -M 0 : σ * (ξ) ] (8)
where ∆ 0 (ξ) is the periodic fourth-order Green's operator for the stress associated to the reference medium M 0 . Ω(ξ) is a third order tensor whose components are given by:

Ω ijk (ξ) = i |ξ| 4 [ (δ ij ξ k + δ ik ξ j + δ jk ξ i )|ξ| 2 -2ξ i ξ j ξ k ] (9)
where |ξ| denotes the norm of the wave vector |ξ| = (ξ i ξ i ) 1/2 . The periodic strain field, ε * (ξ), is defined by: [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF] In the last relation I α (ξ) * σ * (ξ) denotes the convolution product between I α (ξ) and σ * (ξ). The slip condition at the surface Γ is accounted for by taking the limit e/a → 0. It follows that:

ε * (ξ) = ∑ α 1 2µ α [ I α (ξ) * σ * (ξ) ] with : I α (ξ) = 1 V ∫ Vα exp(iξ.x)dx
ε * (ξ) = 1 2µ f [ I f (ξ) * σ * (ξ) + γI g (ξ) * σ * (ξ) ] (11)
with:

I g (ξ) = a V ∫ Γ exp(iξ.x)dx (12)
in which we have introduced γ = c/a. Note that, in [START_REF] Kroner | Statistical Continuum Mechanics[END_REF], the main difference with the case of elastic composites lies in the presence of : (i) the body force f i whose main consequence on the iterative scheme is explained in the following section, (ii) the slip condition introduced by I g (ξ).

The FFT based iterative scheme

Relation ( 8) is the Fourier transform of an integral equation of "Lippman-Schwinger-Dyson" type which can be obtained by the use of the equivalent explicit recurrence process:

σ * i+1 (ξ) = -Ω(ξ).f (ξ) -∆ 0 (ξ) : [ ε * i (ξ) -M 0 : σ * i (ξ) ] ( 13
)
starting from a first term σ * 1 which will be defined later. Using the following equality:

∆ 0 (ξ) : M 0 : σ * i (ξ) = σ * i (ξ) + Ω(ξ).f (ξ) ( 14
)
the recurrence relation ( 13) can also be written in the form:

σ * i+1 (ξ) = σ * i (ξ) -∆ 0 (ξ) : ε * i (ξ) ( 15
)
in which the "body forces", f i , do not appear. The iterative scheme [START_REF] Sanchez-Palencia | Nonhomogeneous Media and Vibration Theory[END_REF] takes the same form as in the case of an elastic composite. However, the main difference lies in the choice of σ * i at the first step in [START_REF] Sanchez-Palencia | Nonhomogeneous Media and Vibration Theory[END_REF]. This stress tensor must be periodic and must equilibrate the macroscopic gradient of pressure, of components G i , within the fluid. It is not therefore divergence free as for elasticity. However, the associated strain rate field may not comply with the compatibility equation. Such a field is chosen, in the Fourier space, as σ * 1 (ξ) = -Ω(ξ).f (ξ) where f i = G i within the fluid but f i remains undefined within the inclusion. If we consider a uniform "body force" f i within the inclusion and considering the equilibrium of the unit cell, < f i > V = 0, it is easy to check that f i is given by:

f i (x) =        G i for x ∈ V f - 1 -c s c s G i for x ∈ V s (16)
or:

f i (ξ) =        - G i c s I s (ξ) ξ i ̸ = 0 0 ξ i = 0 (17)
where c s defines the volume fraction of inclusions. Note that the "body force" within the inclusion is necessary in order to fix the inclusion. Note also that there are other possibilities to fix the inclusion, by considering a line force for example. The local strain is known at convergence of the iterative process and then, one can compute the components of the local velocity field as follows:

v * i (ξ) = 2i |ξ| 2 ε * ij (ξ)ξ j (18)
At this stage, the velocity field is defined by its Fourier coefficients for all values of ξ i except for ξ i = 0. It means that the velocity field is defined up to an added constant. The reason for such a result is that the condition that v * i is null at the surface of the solid has never been taken into account. The velocity field has to be computed from the one given by (18) but one has to specify that v * i = 0 at one point taken within the rigid inclusion.

Application

As an application, we consider a periodic porous medium made up of rigid cylinders (of radius a) aligned in the plane of their cross sections and arranged along a squared lattice of area 4b 2 . In the case of a rigid cylinder located at x i = 0, the shape coefficients I s (ξ) and I g (ξ) are given by:

I s (ξ) = πa 2b 2 |ξ| J 1 (a|ξ|), I g (ξ) = πa 2 2b 2 J 0 (a|ξ|) (19)
where J 0 and J 1 are the Bessel functions of the first kind, and zero and first order, respectively. The flow of the incompressible fluid is due to a prescribed macroscopic pressure G i = G 1 δ i1 . Due to the squared symmetry of the unit cell, the second order permeability tensor is isotropic and can be depicted by a scalar k. The normalized permeability k/(4b 2 ) is plotted in figure 3 as a function of the ratio a/b for three values of γ (where it is recalled that γ = c/a, c being the slip coefficient). As a validation, the FFT based solution (full line) is compared in figure 3, to a Finite Element solution (circles) for γ = 0. The FEM solution has been obtained by using the COMSOL software, the computation being performed with 3584 triangular elements. The porosity is given by ϕ = 1 -πa 2 /(4b 2 ). The maximum theoretically achievable volume fraction of the cylinders is π/4 ≃ 0.785, when a = b, the porosity is then given comprised between ϕ = 0.215 and ϕ = 1. It can be observed that the permeability takes the same value for a = 0 and a = b, whatever the value of

Conclusion

A FFT based-method for computing the permeability of a periodic microporous media saturated by a viscous fluid has been proposed in the present paper. The components of the permeability are obtained by determining the velocity field over the unit cell due to a prescribed macroscopic gradient pressure. The slip effect on the permeability, such as in the case of rarefied gases or micropores, is also studied in this paper. The methodology used to introduce this slip condition is to consider an interphase between the fluid and the solid whose thickness tends to zero. As an illustration purpose, the method is applied for the computation of the permeability of a porous medium made up of aligned rigid cylinders. The results show that the permeability increases with the slip coefficient. It has been also shown that the effect of the slip condition increases with the cylinders radius and is maximal when the inclusions are closely spaced. The method for taking into account interface effects described in this paper can be obviously translated to the case of elasticity for studying nano-composites.

I8n ( 2 )Figure 1 .

 21 Figure 1. Periodic unit cell of the porous medium.

Figure 2 .

 2 Figure 2. The unit cell for periodic array of aligned rigid cylinders

γ.Figure 3 .

 3 Figure 3. Normalized permeability k/(4b 2 ) as function of the ratio a/b.

Figure 4 .

 4 Figure 4. Ratio k(γ = 1)/k(γ = 0) as function of the ratio a/b.