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

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

Interfacial models in viscoplastic composites materials

Keywords: Imperfect Interface, Interphase, Viscoplasticity, Plasticity, Composites

The aim of the present work is to extend the concept of interphase and equivalent imperfect interface in the context of viscoplasticity. The interphase takes the form of a thin curved layer of constant thickness, made up of a rigid viscoplastic material located between two other surrounding materials. We aim at representing this interphase by an interface with appropriately devised interface conditions. To reach this objective, a Taylor expansion of the relevant physical fields in the thin region is used. It is shown that, depending of the degree of stiffness of the layer with respect to the neighboring media, this interphase can be replaced by an idealized imperfect interface involving the jump of the velocity field or the traction vector. The first kind of interface model, applicable to soft interphases, is the "spring-type" interface across which the traction are continuous but the velocity field exhibits a discontinuity which is given in term of the traction by a power-law type relation. Moreover, it is shown that the constant of the model can be expressed in terms of the material parameters of the interphase. When the interphase is stiffer than the two surrounding media, one obtain a "stress-type" interface across which the velocity is continuous and a jump condition for the traction is given by a generalization of the so-called Young Laplace model to viscoplastic solids.

Introduction

The properties of a solid can be significantly affected at or near an interface, as it has been shown from atomistic simulations [START_REF] Miller | Size-dependent elastic properties of nanosized structural elements[END_REF], [START_REF] Diao | Atomistic simulation of the structure and elastic properties of gold nanowires[END_REF][START_REF] Diao | Atomistic simulations of the yielding of gold nanowires[END_REF]. In the last decades, many works deal about the incorporation of these interfacial effects into the overall properties of elastic composites (see for instance Hashin [START_REF] Hashin | Extremum principles for elastic heterogeneous media with imperfect interfaces and their application to bounding of effective moduli[END_REF], Chen et al. [START_REF] Chen | Fibrous nanocomposites with interface stress: Hill's and Levin's connections for effective moduli[END_REF], Duan et al. [START_REF] Duan | Size-dependent effective elastic constants of solids containing nano-inhomogeneities with interface stress[END_REF]). Two kinds of imperfect interface models are commonly used: -the "spring interface" model for which the displacement field is discontinuous and given in terms of the traction vector.

-the "stress interface" model [START_REF] Gurtin | A continuum theory of elastic material surfaces[END_REF][START_REF] Povstenko | Theoretical investigation of phenomena caused by heterogeneous surface tension in solids[END_REF], in which the jump of the traction vector is given by a generalized Young Laplace equation. A rigorous methodology of construction of imperfect interfaces models has been firstly proposed by Sanchez-Palencia [START_REF] Sanchez-Palencia | Comportement limite d'un problème de transmission à travers une plaque faiblement conductrice[END_REF], Pham Huy and Sanchez-Palencia [START_REF] Huy | Phénomènes de transmission à travers des couches minces de conductivité élevée[END_REF] for thermal conduction problem. It consists in replacing a thin interphase by an idealized imperfect interface. This approach has been later generalized in the case of linear elasticity by Karlbring and Mochvan [START_REF] Klarbring | Derivation of a model of adhesively bonded joints by the asymptotic expansion method[END_REF][START_REF] Klarbring | Asymptotic modelling of adhesive joints[END_REF], Bövik [START_REF] Bövik | On the modelling of thin interface layers in elastic and acoustic scattering problems[END_REF], Hashin [START_REF] Hashin | Thin interphase/imperfect interface in elasticity with application to coated fiber composites[END_REF] and Benveniste et al. [START_REF] Benveniste | Imperfect soft and stiff interphases in twodimensional elasticity[END_REF][START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF] (see [START_REF] Benveniste | Imperfect soft and stiff interphases in twodimensional elasticity[END_REF][START_REF] Rubina | A Cosserat shell model for interphases in elastic media[END_REF] for a comprehensive list of publications on the subject). In particular, it has been found that the "spring" or "stress" type interfaces appears as the two limit cases of stiff or soft interphases. Recently, some works dealt about the consideration of imperfect interface in the context of non linear composites materials [START_REF] Zhang | Effect of surface energy on the yield strength of nanoporous materials[END_REF][START_REF] Zhang | Effect of surface stress on the asymmetric yield strength of nanowires[END_REF][START_REF] Zhang | Effect of surface/interface stress on the plastic deformation of nanoporous materials and nanocomposites[END_REF], [START_REF] Dormieux | An extension of Gurson model incorporating stresses effects[END_REF]. For instance, Dormieux et al. [START_REF] Dormieux | An extension of Gurson model incorporating stresses effects[END_REF] postulate a plastic interface model to investigate the effect of interfacial stress on the yield strength of nanoporous materials. In the present study, we derive interface models for viscoplastic materials. To reach this objective, we propose to generalize the Taylor expansion approach initiated by Bövik [START_REF] Bövik | On the modelling of thin interface layers in elastic and acoustic scattering problems[END_REF] and formalized in Benveniste [START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF]. A thin viscoplastic three dimensional curved layer between two surrounding media is considered. This interphase is replaced by an idealized interface involving the jump of both velocity field and traction vector. The two particular cases of a stiff or soft interphase are investigated and allow us to derive a generalization of the "spring" and "stress" interface models in the non linear context. Finally, it is shown that the interface model postulated in [START_REF] Dormieux | An extension of Gurson model incorporating stresses effects[END_REF] appears as a special case of our approach once proper restrictions are applied on the magnitude of the material properties of the interphase.

Statement of the problem

Consider a three dimensional medium constituted of two media separated by a thin layer, called interphase or coating, whose thickness is denoted by h. The thin interphase is delimited by two parallel perfectly bonded interfaces, denoted by S 1 and S 2 , the "middle" surface, at half distance between S 1 and S 2 , being denoted by S 0 (see figure 1). All points within the interphase are represented in the curvilinear coordinates system (x 1 , x 2 , x 3 ). The iso-x 3 surfaces define surfaces which are parallel to S 0 . The surface x 3 = h/2 defines the interface between coating and medium "2", namely S 2 whereas the surface x 3 = -h/2 defines the interface S 1 . The interphase is assumed to be made up of a homogeneous rigid-viscoplastic material whose properties differ from the ones of the two surrounding media. The equations of the model are summarized below:

div(σ) = 0 d = 1 2 [ grad(v) + grad T (v) ] d = ∂ψ ∂σ , ψ(σ) = σ 0 ε0 1 + n ( σ eq σ 0 ) 1+n , tr(d) = 0 (1)
or, dually:

σ = pi + 2σ 0 3ε 0 ( d eq ε0 ) m-1 d (2)
In the above expressions, v is the velocity field, d and σ are the local strain rate and stress tensors. By p, we denote the local pressure, p = tr(σ)/3, i represent the second order identity tensor. σ eq and d eq are the equivalent "von Mises" stress and strain rate and σ 0 and ε0 are two material parameters. The two exponents n and m are classically related by n = 1/m and are defined such that: 0 ≤ m ≤ 1 and 1 ≤ n ≤ +∞. The particular case of an incompressible linear viscous material corresponds to n = m = 1 whereas the limited case n = +∞ (and then m = 0) corresponds to a rigid ideally-plastic material.

We aim now at replacing this interphase by an idealized surface and at deriving the constitutive equations giving the jump of the velocity and of the traction across S 0 . As previously mentioned, the approach uses the Taylor expansion and is detailed in section 4.

Preliminaries

Let us introduce the unit normal vector n taken on each iso-x 3 surfaces and the associated two orthogonal projectors π and π ⊥ defined by:

π = i -n ⊗ n, π ⊥ = n ⊗ n ( 3
)
where i is the second order identity tensor. The strain rate is decomposed into:

d = d s + g s ⊗ n + n ⊗ g s + g n π ⊥ (4)
where:

d s = π.d.π, g s = π.d.n, g n = d : π ⊥ Due to the incompressibility, one has g n = -tr(d s )
, consequently the strain rate can also be written:

d = d s -tr(d s )π ⊥ + g s ⊗ n + n ⊗ g s (5)
The equivalent strain rate can be put into the form:

d eq = (D 2 s + G 2 s ) 1/2 (6)
with:

D s = [ 2 3 ( d s : d s + tr(d s ) 2 ) ] 1/2 , G s = [ 4 3 g s : g s ] 1/2 (7)
Similarly, the stress field σ is decomposed into:

σ = σ s + t s ⊗ n + n ⊗ t s + t n π ⊥ (8) with: σ s = π.σ.π, t s = π.σ.n, t n = σ : π ⊥
There are two kinds of quantities which enter into the expression of the strain rate (5) and of the stress (8):

-quantities which are continuous across a perfectly bonded interface: d s , t s and t n , -quantities which are discontinuous across a perfect interface: σ s and g s . We aim at expressing σ s and g s as function of d s , t s and t n . To do that, we decompose the stress-strain relation as follows:

                           σ s = 2σ 0 3ε 0 ( d eq ε0 ) m-1 d s + pπ t s = 2σ 0 3ε 0 ( d eq ε0 ) m-1 g s t n = - 2σ 0 3ε 0 ( d eq ε0 ) m-1 tr(d s ) + p (9)
From the second relation in ( 9), one has:

t s .t s = σ 2 0 3ε 2 0 ( d eq ε0 ) 2m-2 G 2 s ( 10
)
Introducing into the above equation

T s = √ 3t s .t s ( 11
)
it can be observed that G s is solution of the following non linear equation in z:

(z 2 + D 2 s ) m-1 - ( εm 0 σ 0 ) 2 T 2 s z 2 = 0 (12)
A closed form solution of the above equation can be found in the special case of a linear viscous material (m = 1) and of a perfectly plastic material (m = 0) (see sections 5 and 7) or in the case of specific rational values of m =1 2 and m = 1 3 . In any case 0 < m < 1, we formally denote the solution 1 by z = Z where Z depends on D s , T s and on the material parameters ε0 , σ 0 . The equivalent strain reads:

d eq = (D 2 s + Z 2 ) 1/2
, and from the second relation in [START_REF] Dormieux | An extension of Gurson model incorporating stresses effects[END_REF], one has:

g s = 3ε 0 2σ 0 ( (D 2 s + Z 2 ) 1/2 ε0 ) 1-m t s ( 13
)
The pressure p, can be eliminated from the first and the last relation in [START_REF] Dormieux | An extension of Gurson model incorporating stresses effects[END_REF]. It leads to:

σ s = 2σ 0 3ε 0 ( (D 2 s + Z 2 ) 1/2 ε0 ) m-1 (d s + tr(d s )π) + t n π (14)
In ( 13) and ( 14), the quantities g s and σ s are only expressed in terms of quantities d s , t s and t n .

The Taylor expansion approach

As explained in the introduction of the paper, an imperfect interface can be seen as the limit case of a thin stiff or soft interface. The connection between the properties of the interphase and of the equivalent interface model has been amply studied by Benveniste [START_REF] Benveniste | Imperfect soft and stiff interphases in twodimensional elasticity[END_REF] in the case of a two dimensional elastic curved surface and has been later generalized in the three dimensional context [START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF].

In order to establish the connection between the properties of the thin interphase and of the imperfect interface, two methods have been applied in past studies. The first one uses the asymptotic expansion. This method has been initiated by Sanchez-Palencia [START_REF] Sanchez-Palencia | Comportement limite d'un problème de transmission à travers une plaque faiblement conductrice[END_REF] and Pham Huy and Sanchez-Palencia [START_REF] Huy | Phénomènes de transmission à travers des couches minces de conductivité élevée[END_REF] for conductivity, and applied by Klarbring et al. [START_REF] Klarbring | Derivation of a model of adhesively bonded joints by the asymptotic expansion method[END_REF][START_REF] Klarbring | Asymptotic modelling of adhesive joints[END_REF], Avila-Pozos et al. [START_REF] Avila-Pozos | Asymptotic model of orthotropic highly inhomogeneous layered structure[END_REF] and Benveniste [START_REF] Benveniste | Imperfect soft and stiff interphases in twodimensional elasticity[END_REF] in the context of linear elasticity. The second method has been initiated by Bövik [START_REF] Bövik | On the modelling of thin interface layers in elastic and acoustic scattering problems[END_REF] and formalized by Hashin [START_REF] Hashin | Thin interphase/imperfect interface in elasticity with application to coated fiber composites[END_REF] and Benveniste [START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF]. It uses the Taylor expansion and adequate surface differential operators which make possible to derive the general expressions of the jumps of the physical fields. In the special cases of a plane layer [START_REF] Klarbring | Asymptotic modelling of adhesive joints[END_REF][START_REF] Avila-Pozos | Asymptotic model of orthotropic highly inhomogeneous layered structure[END_REF] and of a two dimensional curved layer [START_REF] Benveniste | Imperfect soft and stiff interphases in twodimensional elasticity[END_REF], the results obtained by the asymptotic expansion method are consistent with the one derived form the Taylor expansion approach [START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF].

Although it has not yet been proven, it is likely that this equivalence between the two methods remains true in the general case of an anisotropic three dimensional curved layer.

Of course, the use of the Taylor expansion based method suppose that the physical fields which are expanded are continuously differentiable, defined and regular everywhere in the thin layer. Such conditions are not needed when the asymptotic expansion method is use, it constitutes the main difference between the two methods. However, the advantage of the Taylor expansion approach over the asymptotic one is that it leads to a compact representation in which the particular cases of "spring interface" and of "stress interface" appear as limit cases.

In the prensent paper we extend the Taylor expansion method to the context of viscoplasticity. The first step consist in expanding the velocity field, v, and the traction vector t at the mid surface, S 0 , about the lower interface S 1 or upper interface S 2 :

v| S 0 = v| S 1 + h 2 grad n (v)| S 1 + o(h 2 ) v| S 0 = v| S 2 - h 2 grad n (v)| S 2 + o(h 2 ) ( 15
)
where grad n (v) = grad(v).n are the derivatives of v along the x 3 direction. In [START_REF] Klarbring | Asymptotic modelling of adhesive joints[END_REF] the series is truncated at the first order, higher terms are neglected in the present work. It supposes that the thickness of the layer is small compared to the characteristic length of the microstructure, as for example the size of the inclusions for a composite material. A combination of the two equations in [START_REF] Klarbring | Asymptotic modelling of adhesive joints[END_REF] gives:

v| S 2 -v| S 1 = h 2 [ grad n (v)| S 1 + grad n (v)| S 2 ] + o(h 2 ) (16)
Similarly, an equivalent equation can be written for the traction vector t:

t| S 2 -t| S 1 = h 2 [ grad n (t)| S 1 + grad n (t)| S 2 ] + o(h 2 ) (17)
In ( 15) and ( 17), we are dealing with continuously differentiable fields, defined and regular everywhere in the thin layer. Following the lines of Benveniste [START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF] we aim at expressing the differential operators, grad n (v) and grad n (t) in ( 16) and ( 17) as functions of the surface derivatives of the velocity, grad s (v) = grad(v).π and of the traction, t. Once those operators are determined, the derivation of the interface equations will be easier. Firstly, the normal and tangential components of grad n (v) are given by:

grad n (v).n = g n = -tr(d s ) grad n (v).π = 2g s -n. grad s (v) (18)
where g s is given by [START_REF] Hashin | Thin interphase/imperfect interface in elasticity with application to coated fiber composites[END_REF]. It follows that:

grad n (v) = 3ε m 0 σ 0 (D 2 s + Z 2 ) (1-m)/2 t s -n. grad s (v) -tr(d s )n ( 19
)
Using now the equilibrium equation, one has:

grad n (t) = grad(σ) : π ⊥ = -grad(σ) : π (20)
In the above expression, the decomposition ( 8) is used, in which σ s is replaced by [START_REF] Klarbring | Derivation of a model of adhesively bonded joints by the asymptotic expansion method[END_REF]. It follows that:

grad n (t) = -grad { 2σ 0 3ε m 0 ( D 2 s + Z 2 ) (m-1)/2 (d s + tr(d s )π) +t s ⊗ n + n ⊗ t s + t n i } : π (21)
In ( 19) and ( 21) the "normal" gradient of v and t are given in term of v and t and their surface derivatives which are equal to their corresponding quantities in the surrounding media "1" and "2".

In the configuration of figure 1, the interphase has been eliminated and replaced by an idealized interface S 0 . Equations ( 16), ( 17) with ( 19) and ( 21) provide the general expressions for the jump of the velocity field and traction across the equivalent imperfect interface S 0 . We aim now at examining the special cases of a soft or stiff interphase for which explicit expressions of the interface model can be derived.

The linear case: m = 1

Consider the particular case of a linear viscous material, n = m = 1. Expressions ( 19) and ( 21) become:

grad n (v) = 1 µ t s -n. grad s (v) -tr(d s )n grad(t).n = -grad {2µ(d s + tr(d s )π) + t s ⊗ n + n ⊗ t s + t n i} : π (22)
in which we have introduced µ = σ 0 /(3ε 0 ). The traction vector and the gradient of the velocity are assumed to takes finite values at the interfaces S 1 and S 2 . Consequently, if µ remains finite when the limit h → 0 is taken, the terms at the right of the equality in ( 16) and ( 17) vanishes when h → 0. The term containing h/2 in both expressions ( 16) and ( 17) does not vanish when h → 0 int the following special cases:

-If µ is assumed to be of the same magnitude of the width h, the interphase is soft and can be replaced by a "spring interface" model which involves a jump of the velocity:

[v] S 0 = 1 μ t s , [t] S 0 = 0 (23)
where μ = hµ. It can be observed that only the tangential components of the velocity field are discontinuous across the interface whereas the normal components remain continuous, this is a consequence of the material incompressibility. The above results are consistent with the one provided in [START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF] when a linear incompressible interphase is considered.

-When µ is of the order of magnitude of 1/h, the interphase is stiff and can be replaced by "stress interface" model which involves a jump of the traction:

[v] S 0 = 0, [t] S 0 = -grad(τ ) : π, τ = 2µ(d s + tr(d s )π) (24)
with μ = µ/h and τ is called interfacial stress. The jump of the traction takes the form of a generalized Young Laplace Law. Again, The above results are consistent with the one provided in [START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF] for an elastic incompressible interphase.

6 The viscoplastic case: 0 < m < 1

The value of the exponent m is now arbitrary taken in the interval ]0, 1[. We first consider the case for which ε0 is proportional to 1/h, and we put ε0 = ε0 /h. σ 0 and ε0 are assumed to remains finite when h → 0. The non linear equation ( 12) become:

(z 2 + D 2 s ) m-1 - 1 h 2m-2 ( εm 0 σ 0 ) 2 T 2 s z 2 = 0 (25)
When low values of h are considered in equation ( 12), the solution for z is obtained by using an asymptotic expansion. The solution z is taken as: z = z -1 /h + z 0 + z 1 h + .... The computation of z -1 , z 0 ... gives:

z -1 = ε0 ( T s σ 0 ) n , z 0 = 0, z 1 = 1 2 (1 -n) D 2 s ε0 ( T s σ 0 ) -n , etc. ( 26
)
Only the higher order term, z -1 is kept in the expansion of z, as a consequence:

Z = ε0 h ( T s σ 0 ) n + o(h) ( 27
)
Replacing this expression in [START_REF] Rubina | A Cosserat shell model for interphases in elastic media[END_REF] and [START_REF] Zhang | Effect of surface energy on the yield strength of nanoporous materials[END_REF], and taking in ( 16) and ( 17) the limit h → 0 leads to:

[v] S 0 = 3ε 0 σ 0 ( T s σ 0 ) n-1 t s , [t] S 0 = 0 (28)
Dually, we suppose that the reference stress σ 0 of the coating is proportional to 1/h and we put σ 0 = σ0 /h, (σ 0 and ε0 are assumed to remains finite when h → 0). When low values of h are considered in equation ( 12), the solution z can be taken into the form: z = z 1 h + z 2 h 2 + .... Only the term z 1 is kept in the expansion of z, and solution for Z can be read:

Z = hε m 0 D 1-m s T s σ 0 + o(1) (29)
Taking in ( 16) and ( 17) the limit h → 0 leads to:

[v] S 0 = 0, [t] S 0 = -grad(τ ) : π with: τ = 2σ 0 3ε 0 ( D s ε0 ) m-1 (d s + tr(d s )π) (30)
As in the linear case, the jump of the traction takes the form of a Young Laplace law. However, the relation giving the interfacial stress, τ as a function of the surface strain rate is given by a surface viscoplastic law under plane stress conditions.

7 The plastic case: m = 0

We now look for the particular case of a rigid perfectly-plastic case corresponding to n = +∞ or m = 0. Equation ( 12) can be explicitly solved and gives:

Z = T s D s √ σ 2 0 -T 2 s (31)
Replacing (31) in ( 19) and ( 21) leads to:

grad n (v) = 3D s √ σ 2 0 -T 2 s t s -n. grad s (v) -tr(d s )n grad(t).n = -grad    2 3 √ σ 2 0 -T 2 s D s (d s + tr(d s )π) +t s ⊗ n + n ⊗ t s + t n i} : π (32)
Since the material parameter ε0 vanishes in the case of a plastic material, the conditions [v] S 0 ̸ = 0 and [t] S 0 = 0 cannot be obtained as the limit case of a thin plastic layer. However, when the reference stress is chosen as σ0 = σ 0 /h, the above equations gives:

grad n (v) = o(1) grad(t).n = - 1 h grad { 2σ 0 3D s (d s + tr(d s )π) } : π + o(1) (33)
The limit h → 0 is taken in ( 16) and ( 17), giving: The obtained model coincides with the one postulated in [START_REF] Dormieux | An extension of Gurson model incorporating stresses effects[END_REF].

[v] S 0 = 0, [t] S 0 = -div s (τ ), τ = 2σ 0 3D s (d s + tr(d s )π) (34

Conclusion

Interface models, commonly used in elasticity, have been generalized, in this paper, to the context of viscoplasticity. The methodology, applied to derive the constitutive equations of the interface, uses the concept of equivalent interphase model and Taylor expansions. It has been shown that, when the thickness of the coating, h, is very small compared with the characteristic length of the microstructure, the coating can be replaced by an idealized interface involving the jump of the velocity field and of the traction vector. Two particular cases has been investigated: -When the interphase is soft compared with the surrounding media, the obtained interface model involves a jump of the velocity field whereas the traction vector remains continuous. The jump of the velocity field across the interface is given in terms of the surface component of the traction vector and takes the form of a power-law type relation.

-When a stiff interphase is considered, the equivalent interface model involves a jump of the traction vector whereas the velocity field remains continuous. The relation giving the jump across the interface takes the form of the Young-Laplace equation. However, the interfacial stress-strain relation is given by a viscoplastic law under plane stress conditions. These interface models can be applied for taking into account the surface effects in viscoplastic composites.

) 2 (

 2 The dissipation reads τ : d s = σ 0 D s . The surface strain rate, d s , can be eliminated in the relation giving τ as a function of d s : 3

Equation (12) gives a positive and a negative real solution but also complex solutions, here the solution Z makes reference to the positives real solution.


Thu-Huong Tran
email: thu-huong.tran@univ-paris-est.fr

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

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

A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media

Keywords: Gradient elasticity, Homogenization, Periodic composites, Series expansion method, Generalized Hill-Mandel lemma

A micromechanics-based approach for the derivation of the effective properties of periodic linear elastic composites which exhibit strain gradient effects at the macroscopic level is presented. At the local scale, all phases of the composite obey the classic equations of three-dimensional elasticity, but, since the assumption of strict separation of scales is not verified, the macroscopic behavior is described by the equations of strain gradient elasticity. The methodology uses the series expansions at the local scale, for which higher-order terms, (which are generally neglected in standard homogenization framework) are kept, in order to take into account the microstructural effects. An energy based micro-macro transition is then proposed for upscaling and constitutes, in fact, a generalization of the Hill-Mandel lemma to the case of higher-order homogenization problems. The constitutive relations and the definitions for higher-order elasticity tensors are retrieved by means of the "state law" associated to the derived macroscopic potential. As an illustration purpose, we derive the closed-form expressions for the components of the gradient elasticity tensors in the particular case of a stratified periodic composite. For handling the problems with an arbitrary microstructure, a FFT-based computational iterative scheme is proposed in the last part of the paper. Its efficiency is shown in the particular case of composites reinforced by long fibers.

Introduction

Conventional homogenization processes applied to periodic composite materials generally use the assumption of scale separation, i.e. the size of microstructural elements remains negligible compared with the size of the studied mechanical components or with the characteristic length of the applied load. Under this hypothesis, the macroscopic quantities are nearly constant at local scale (the scale of the heterogeneities) and the effective properties, obtained from the homogenization procedure, do not introduce the size of the microstructural elements, typically, the size of inclusions or the size of the periodic unit cell. This framework is valid in many situations, specially when the heterogeneities are micrometer sized and when the gradient of the macroscopic strain remains small enough. When such assumptions are not satisfied, the size of the microstructural elements affects the macroscopic behavior of composites. Moreover, these effects can become predominant for some special practical problems in mechanics and of high importance for structural design such as in the case of fracture (see [START_REF] Exadaktylos | Gradient elasticity with surface energy: Mode-I crack problem[END_REF], Altan and Aifantis 1992, [START_REF] Georgiadis | Energy theorems and the J-integral in dipolar gradient elasticity[END_REF] or the localization phenomenon (see [START_REF] Triantafyllidis | A gradient approach to localization of deformation[END_REF][START_REF] Chambon | Plastic continuum with microstructure, local second gradient theories for geomaterials: localization studies[END_REF][START_REF] Engelen | An evaluation of higher-order plasticity theories for predicting size effects and localization[END_REF]) for which very large values of strain gradient may occur. Generalized continuum-mechanics theories introduce a characteristic length of the microstructure, [START_REF] Toupin | Elastic materials with couple-stresses[END_REF][START_REF] Toupin | Theories of elasticity with couple-stress[END_REF][START_REF] Mindlin | Micro-structure in linear elasticity[END_REF], [START_REF] Mindlin | On first strain gradient theory in linear elasticity[END_REF], (see Fleck andHutchinson 1993, Fleck et al. 1994, Mühlhaus and Aifantis 1991 for a plastic version of the strain gradient models). These approaches extend the usual linear elasticity by including the gradient of strain and higher-order derivatives of the strain in the expression of the elastic energy density. These models are often considered as describing more adequately the effect of the microstructure at the macroscale. However they are phenomenological in the sense that the constitutive elastic coefficients are not derived from a micromechanics-based framework. Moreover, such models generally introduce a large number of coefficients whose identification constitutes an important limitation for their use in the modeling of structures. The derivation of the constitutive equations of non local models in a homogenization framework has been the subject of extensive research during the last three decades. Among the first, [START_REF] Diener | Bounds for the non-local effective properties of random media[END_REF][START_REF] Diener | Bounds for the non-local effective properties of random media II[END_REF][START_REF] Diener | Bounds on the non-local effective elastic properties of composites[END_REF] and later Drugan and Willis 1996 derived a non local constitutive model from the variational Hashin-Shtrikman variational principle. Other works, [START_REF] Bouyge | A micromechanically based couple-stress model of an elastic two-phase composite[END_REF], 2002[START_REF] Kouznetsova | Multi-scale constitutive modeling of heterogeneous materials with gradient enhanced computational homogenization scheme[END_REF][START_REF] Kouznetsova | Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy[END_REF][START_REF] Yuan | A micromechanical approach of nonlocal modeling for media with periodic microstructures[END_REF], have used higher-order boundary type conditions for the unit cell in order to derive the constitutive equations of generalized continuum models. These approaches use a generalization of the Hill-Mandel´s macrohomogeneity condition for deriving the higherorder constitutive relations and can deal with non-linear materials while the present paper offers only developments in the linear elastic case. However, as mentioned by [START_REF] Yuan | A micromechanical approach of nonlocal modeling for media with periodic microstructures[END_REF], these definitions produce overevaluated values of the macroscopic internal energy of gradient elastic media and give rise to some unphysical results. Indeed, the macroscopic behavior obtained from this procedure still leads to a macroscopic gradient elastic model even if the unit cell is defined by a homogeneous material. Still in the context of periodic microstructures, Gambin and Kroener 1989, Boutin 1996, Triantafyllidis and Bardenhagen 1996 and later Smyshlyaev and Cherednichenko 2000 have studied the influence of the higher-order terms of the series expansion on the macroscopic behavior of linear elastic composites. This approach of series expansion, which has been initiated by Benssoussan et al. 1978 and[START_REF] Sanchez-Palencia | Nonhomogeneous media and vibration theory[END_REF], has shown to be a rigorous and efficient method for introducing the effect of the macroscopic gradient of strain on the local response of linear composites. The asymptotic expansion method introduces the scale factor ϵ defined as the ratio between the characteristic length of the microstructure and the one of the applied macroscopic loading. When this scale factor ϵ is very small compared with 1, there is a strict separation of the micro and macroscopic scales. In practice, the limit ϵ → 0 is taken in the expression of strains and stresses, and the standard homogenization framework can be used. When the parameter ϵ is close to 1, no homogenization is valid anymore. When ϵ is lower than 1 but not negligible before 1, the solution can be approximated by keeping higher-order terms in the series expansion. All these terms are then obtained by solving a hierarchy of higher-order elasticity problems with prescribed body forces and eigenstrains whose expressions depend on the solution at the lower-order. In Boutin 1996, the author does not interpret the macroscopic response of the composite as of strain gradient type but only considers the higher-order terms of the series in order to evaluate the deviation of the macroscopic stress-strain relation from the classic one, due to the microstructural effects. Later Smyshlyaev and Cherednichenko 2000 also studied the higher-order terms of the asymptotic expansion in the framework of variational principles and showed that the resulting "homogenized model", with the appropriate definitions for the macroscopic quantities, leads to constitutive equations which are in agreement with the gradient elasticity theories. In the present paper we propose a simple method for computing the effective properties of elastic composites which exhibit a gradient effect at the macroscopic scale which depends on a characteristic length of the microstructure. The method follows the "asymptotic based" approaches of Boutin 1996, Smyshlyaev and[START_REF] Smyshlyaev | On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media[END_REF] The transition between the micro and macro scales is effected by defining the macroscopic elastic energy density as the average of the local energy density: this can be interpreted as a generalization of the Hill-Mandel Lemma for higher order homogenization problems, similarly to the works of [START_REF] Bouyge | A micromechanically based couple-stress model of an elastic two-phase composite[END_REF][START_REF] Yuan | A micromechanical approach of nonlocal modeling for media with periodic microstructures[END_REF]. The macroscopic law of the composite is then obtained by means of the "state law" associated to the macroscopic potential. As an illustration, we derive the closed-form solution for the higher-order elastic coefficients in the case of a periodic stratified composite. For handling the problem of an arbitrary geometry of the unit cell, we extend the FFT-based method initially introduced for elasticity by [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF]. The ability of this numerical approach is illustrated in the case of a composite made up of periodically aligned fibers along a regular network.

Higher-order homogenization problems and their solutions

In this section, we recall shortly without new developments the asymptotic expansion method applied to periodic composites and its elastic solution as derived in the work of [START_REF] Boutin | Microstructural effects in elastic composites[END_REF], because this work constitutes the basis of the developments proposed in the present paper. Let us consider a periodic microstructure defined by a parallelepipedic unit cell Y and three vectors denoted by d 1 , d 2 , d 3 which characterize the translational invariance. Each constituent obeys the equations of three-dimensional elasticity:

             ∇ z .σ + f (z) = 0 σ = c(z) : ε ε = 1 2 (∇ z ⊗u + u ⊗ ∇ z ) (1)
where z = (z 1 , z 2 , z 3 ) denotes the vector position. By c(z), we denote the fourth-order elasticity tensor of the periodic composite which is piecewise constant and invariant by any translation according to the vectors d 1 , d 2 , d 3 . The inclusions and the matrix are assumed to be perfectly bonded through their interfaces. Then, the displacement field, u, and the traction vector, t = σ.n, remain continuous across the interface between the matrix and the heterogeneities. The problem is defined by its two characteristic length scales, h and L. The first, h, is characteristic of the microstructural elements (the size of the unit cell, of the inclusions, the distance between two neighboring inclusions...). L denotes a characteristic length of the studied macrostructure or the characteristic size related to the applied load. The ratio between these two lengths defines the scale factor, ϵ = h/L which is assumed to be smaller than 1 but not negligible. The existence of a small parameter in (1) allows to search the solution along asymptotic series. The technique basically consists in four steps:

• The presence of two length scales suggests to introduce the non dimensional space variables x = z/L and y = z/h, where x is the "slow" oscillating variable while y is the "fast" oscillating variable. • The displacement field is expanded as a power series in ϵ:

u(x, y) = L +∞ ∑ i=0 ϵ i u i (x, y) (2)
where the displacements fields u i (x, y) are Y -periodic and non dimensional due to the presence of the lengthscale L before the sum in the previous equation. • The gradient operator is decomposed into two parts :

∇ z = 1 L (∇ x + 1 ϵ ∇ y )
, where the indices x and y indicate that the derivatives are taken respectively with respect to x and y.

• By equating terms with the same power of ϵ, we get a hierarchy of equations for the quantities u 0 (x, y), u 1 (x, y), u 2 (x, y), etc.

The total displacement field takes the form:

u(x, y) = U (x) + ϵX 1 (y) : E(x) + ϵ 2 X 2 (y) :• G(x) + ϵ 3 X 3 (y) :: D(x)... ( 3
)
In which quantities U (x), E(x), G(x) and D(x) are functions of the variable x while tensors X 1 (y), X 2 (y), ... are only functions of the variable y. Tensors X 1 (y), X 2 (y), ... are obtained by solving a hierarchy of linear elasticity problems which are given in appendix A. The physical meaning of quantities U (x), E(x), G(x) and D(x) will be specified in the next section.

The associated strain field is:

ε(x, y) = a 0 (y) : E(x) + ϵa 1 (y) :• G(x) + ϵ 2 a 2 (y) :: D(x)... ( 4
)
where the components of tensors a 0 (y), a 1 (y), etc., are related to the ones of tensors X 1 (y), X 2 (y), etc., by:

a 0 ijpq (y) = 1 2 (δ ip δ jq + δ iq δ jp ) + 1 2 { ∂X 1 ipq ∂y j (y) + ∂X 1 jpq ∂y i (y) } a 1 ijpqr (y) = 1 2 { X 1 ipq (y)δ jr + X 1 jpq (y)δ ir } + 1 2 { ∂X 2 ipqr ∂y j (y) + ∂X 2 jpqr ∂y i (y) } etc. (5)
The local stress field reads:

σ(x, y) = c 0 (y) : E(x) + ϵc 1 (y) :• G(x) + ϵ 2 c 2 (y) :• D(x)... ( 6
)
with c α (y) = c(y) : a α (y). The local fields then require the identification of the components of tensors X α (y) for α ≥ 1.

The micro-macro transition

In this section, we propose a homogenization process based on the elastic energy density which proves to be adapted for taking properly into account higher order terms of the asymptotic series and the microstructural effects at the macroscopic level. As in the standard thermodynamic framework, the macroscopic law is obtained by means of the state law associated to the derived expression of the homogenized elastic potential. It is shown that the resulting macroscopic model complies with the first gradient theory of [START_REF] Toupin | Elastic materials with couple-stresses[END_REF][START_REF] Toupin | Theories of elasticity with couple-stress[END_REF] and [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF][START_REF] Mindlin | On first strain gradient theory in linear elasticity[END_REF], when the series is truncated at the second order, and allows to recover the more general gradient elastic theory established by Green and Rivlin 1964, when all the terms of the series are kept.

Definitions for the stress and hyperstress

Since all tensors X α (y) in ( 3) have a null volume average over a period, it follows that:

< u(x, y) > V = U (x) (7)
U (x) represents the macroscopic displacement field which can be interpreted as the translation of the geometric center of the unit cell. Taking into account the following equalities:

< a 0 (y) > V = I; < a α (y) > V = 0 α = 1, 2, ... ( 8
)
where I denotes the Fourth order identity tensor, it follows that the average of the local strain is given by:

< ε(x, y) > V = E(x) = 1 2 (∇ x ⊗ U (x) + U (x) ⊗ ∇ x ) (9)
E(x) is then interpreted as the macroscopic strain and G(x) and D(x) are respectively the non-dimensional gradient and double gradient of macroscopic strain, the derivative being taken according to the non dimensional macroscopic vector position, x. It follows that, in all the expressions, ϵG(x) and ϵ 2 D(x) must be replaced by: h∇E(x) and h 2 ∇ 2 E(x). Taking into account the expression of the total strain and stress fields (see equations ( 4), ( 6)), the average of the local density energy over the volume of the unit cell can be written in the form: 2

< σ(x, y) : ε(x, y) > V = 1 2 { E(x) : C 0,0 : E(x) + 2hE(x) : C 0,1 : ∇E(x) + h 2 ∇E(x) : C 1,1 : ∇E(x) + 2h 2 E(x) : C 0,2 : ∇ 2 E(x) + 2h 3 ∇E(x) : C 1,2 : ∇ 2 E(x) + 2h 3 ∇E(x) : C 0,3 : ∇ 3 E(x) + . . . } (10)
with:

C α,β i..jp..q =< c α kli..j (y)a β klp..q (y) > V =< c klmn (y)a α mni..j (y)a β klp..q (y) > V ( 11
)
where the C α,β are tensors of order α + β + 4. Note that tensor C β,α can be seen as the "transpose" of C α,β obtained by permuting the first α + 2 indices with the last β +2 indices. The quantity W (x) is called thereafter macroscopic elastic energy. It can be observed that W (x) is a quadratic function according to the macroscopic quantities E(x), ∇E(x), ∇ 2 E(x), etc. As usually in the thermodynamic framework, it is possible to define the dual variables associated to E(x), ∇E(x), ∇ 2 E(x), by means of the state law:

Σ(x) = ∂W ∂E =< σ(x, y) : a 0 (y) > V = C 0,0 : E(x) + hC 0,1 :• ∇E(x) + h 2 C 0,2 :: ∇ 2 E(x) + ... T (x) = ∂W ∂∇E = h < σ(x, y) : a 1 (y) > V = hC 1,0 : E(x) + h 2 C 1,1 :• ∇E(x) + h 3 C 1,2 :: ∇ 2 E(x) + ... M (x) = ∂W ∂∇ 2 E = h 2 < σ(x, y) : a 2 (y) > V = h 2 C 2,0 : E(x) + h 3 C 2,1 :• ∇E(x) + h 4 C 2,2 :: ∇ 2 E(x) + ... etc. (12)
Σ(x) is the macroscopic stress and T (x), M (x), etc., are called hyperstresses. Tensors C α,β are the elastic tensors of the generalized continuum whose determination requires the computation of functions X α (y) for α = 1, 2, 3.... Let us recall that these tensors are determined by solving successively a hierarchy of linear elastic problems also called "higher order homogenization problem" in the literature and which can be found in appendix A. The problems can be solved by standard numerical tools as the Finite Element Method (FEM) or, as proposed in section 5, by a Fast Fourier Transform Based numerical approach, adapted to the context of the computational homogenization of periodic structures.

It must be emphasized that the macroscopic stress is not defined by the classic volume average, Σ(x) ̸ =< σ(x, y) > V . This constitutes an important difference with Boutin 1996, together with the fact that the hyperstresses have not been introduced in Boutin's work. With the above definitions, it follows that:

< σ(x, y) : ε(x, y) > V = Σ(x) : E(x) + T (x) :• ∇E(x) + M (x) :: ∇ 2 E(x) + ... ( 13
)
which constitutes a generalization of the Hill-Mandel lemma to higher-order homogenization framework. For completeness, we have to demonstrate that the stress and hyperstresses [START_REF] Exadaktylos | Gradient elasticity with surface energy: Mode-I crack problem[END_REF] comply with the balance equation for such models of generalized continua. This is the subject of the following sub-section. Before proceeding to this, let us provide the following remarks: (a) When the limit h → 0 of the macroscopic potential ( 10) is taken, all terms associated to the gradient of strain and to higher order derivatives vanish and the energy density reduces to W (x) = 1 2 E(x) : C 0,0 : E(x). The macroscopic elastic law reduces to Σ(x) = C 0,0 : E(x). Obviously, the resulting macroscopic model is of Cauchy type and the standard equations for homogenization are recovered. (b) When all the series are truncated at the order 1 in the expression of the microscopic strain field (4), the macroscopic elastic energy density W (x) computed from [START_REF] Drugan | A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites[END_REF] takes the form:

W (x) = 1 2 E(x) : C 0,0 : E(x) + hE(x) : C 0,1 : ∇E(x) + h 2 2 ∇E(x) : C 1,1 : ∇E(x) (14)
which can be interpreted as the energy density introduced in the first strain gradient theories of [START_REF] Toupin | Elastic materials with couple-stresses[END_REF][START_REF] Toupin | Theories of elasticity with couple-stress[END_REF] and [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF][START_REF] Mindlin | On first strain gradient theory in linear elasticity[END_REF]). Note also that a truncature at the order 2 in the expression of the microscopic strain field (4), leads to a macroscopic energy density which is quadratic with respect to the macroscopic strain, gradient of strain and double gradient of macroscopic strain. Such a potential can be interpreted as the energy density of the second strain gradient theory of [START_REF] Mindlin | Second gradient of strain and surface-tension in linear elasticity[END_REF].

(c) When all the terms of the series are kept in the expression of the microscopic strain field (4), the macroscopic energy density is the one introduced in a more general theory established by [START_REF] Green | Multipolar continuum mechanics[END_REF] which includes the derivatives of the strain at any order.

The macroscopic balance equation

In this section we prove that the definitions of the stress and of the hyperstresses ( 12) are compatible with the balance equation of elastic strain gradient continua. To this end, let us start from the local balance equation with the decomposition

∇ z = 1 L (∇ x + 1 ϵ ∇ y ): ∇ x .σ + 1 ϵ ∇ y .σ + f (x, y) = 0 (15)
The traction σ.n at the boundary of the unit cell is Y -antiperiodic. The volume average of the quantity ∇ y .σ is then equal to zero and:

∇ x . < σ > V +F (x) = 0 ( 16
)
where we have denoted F (x) =< f (x, y) > V and < σ > V is only a function of the macroscopic variable x. In the classic homogenization framework, the series in ( 6) is truncated at the zero order and < σ > V represents the macroscopic stress tensor, Σ(x). In that case, equation ( 16) reduces to ∇ x .Σ(x)+F (x) = 0 which means that the macroscopic law of the composite is of Cauchy type. When higher order terms are kept in the series (6), the quantity < σ > V cannot be interpreted as the macroscopic stress since the definitions ( 12) are now used. However, a representation of the average < σ > V in terms of the macroscopic quantities Σ(x), T (x), M (x) ... is possible. Starting from the definitions of the macroscopic stress [START_REF] Exadaktylos | Gradient elasticity with surface energy: Mode-I crack problem[END_REF] and of the localization tensor a 0 (y) given by ( 5), one has, for α ≥ 0:

1 V ∫ V σ pq a α pqi..jk (y)dV = 1 V ∫ V σ pk X α pi..j (y)dV + 1 V ∫ V σ pq ∂X α+1 pi..jk ∂y q (y)dV (17)
Where X 0 ij is conventionally taken as X 0 ij (y) = δ ij . Using the divergence theorem and accounting for ∂/∂y i = -ϵ∂/∂x i , one has :

1 V ∫ V σ pq ∂X α+1 pi..jk ∂y q (y)dV = 1 V ∫ ∂V σ pq n q X α+1 pi..jk (y)dV +ϵ ∂ ∂x q { 1 V ∫ V σ pq X α+1 pi..jk (y)dV } (18)
The first term at the right of the above equality is null, due to the periodicity of X α+1 i..j (y) and the antiperiodicity of the traction vector σ ij n j . It follows that:

1 V ∫ V σ pq a α pqi..jk (y)dV = 1 V ∫ V σ pk X α pi..j (y)dV +ϵ ∂ ∂x q { 1 V ∫ V σ pq X α+1 pi..jk (y)dV } (19)
Writing now this equation for α = 0, 1, 2... we obtain:

Σ ij (x) = 1 V ∫ V σ ij dV + ∂ ∂x q { ϵ V ∫ V σ pq X 1 pij (y)dV } T ijk (x) = h V ∫ V σ pk X 1 pij (y)dV + ∂ ∂x q { hϵ V ∫ V σ pq X 2 pijk (y)dV } M ijkl (x) = h 2 V ∫ V σ pl X 2 pijk (y)dV + ∂ ∂x q { h 2 ϵ V ∫ V σ pq X 3 pijkl (y)dV } ( 20
)
By combination of the above relations, we obtain:

< σ ij > V = Σ ij (x) - ∂T ijk ∂x k (x) + ∂ 2 M ijkl ∂x k ∂x l (x) -... (21)
The macroscopic balance equation reads then:

∂Σ ij ∂x j (x) - ∂ 2 T ijk ∂x j ∂x k (x) + ∂ 3 M ijkl ∂x j ∂x k ∂x l (x) -... + F i (x) = 0 ( 22
)
Which is exactly the balance equation of the general theory of gradient elastic media developed by [START_REF] Green | Multipolar continuum mechanics[END_REF]. When the series for the strain field is truncated at the first order, the macroscopic potential is given by ( 14) and the effective law only displays a relation between the macroscopic stress Σ(x) and the first hyperstress T (x) with the macroscopic strain, E(x) and the first gradient of strain, ∇E(x). Higher order hyperstresses vanish and, in that case, the macroscopic balance equation reduces to:

∂Σ ij ∂x j (x) - ∂ 2 T ijk ∂x j ∂x k (x) + F i (x) = 0 ( 23
)
Which is exactly the balance equation of first gradient elastic media.

It must be noted that Smyshlaev & Cherednichenko (2000) developed a variational approach for deriving new higher-order effective relations for linear elastic composites which include the microstructural effects and the gradients of the macroscopic strain. By taking the variation of their elastic functional, Smyshlaev & Cherednichenko (2000) observe that balance equation does not reduce to the standard, Cauchy-type, one but introduces new additional terms. As mentioned by the authors themselves, the resulting balance equation may be recognized as the one introduced by the Toupin and Mindlin theories. However these additional terms introduce higher order homogenized stress tensors (which are only second-order tensors) but are not related to hyperstresses which are tensors of order 3, 4, etc. as in the present work.

Closed-form solutions for a stratified composite

As a first illustration we propose to derive the closed-form expressions for the components of gradient elasticity tensors in the case of a stratified composite. The material, depicted on Figure 1, is constituted of two isotropic layers a and b, having the thickness (1 -τ )h and τ h with 0 ≤ τ ≤ 1. The volume fractions of phases a and b are respectively τ and 1 -τ . The elastic moduli of the phase a are λ 1 , µ 1 while those of the layer b are denoted λ 2 and µ 2 . The layers a and b are periodically distributed along direction Oy 1 . The material remains unchanged by any translation along directions Oy 2 and Oy 3 which implies that the local displacements at any given order are only functions of the variable y 1 (which is thereafter denoted y for simplicity). Still for simplicity reasons, the "slow" variable x is omitted in the expression of the macroscopic quantities. Exact expressions of functions X 1 (y) and X 2 (y), which are needed for the computation of effective elastic properties, have been already derived in Boutin 1996. However, due to the presence of misprints in the results, the corrected expressions are provided in appendix B. We now provide the expressions for the higher order elasticity tensors. As for the case of fourth order tensors, a matricial representation is useful for representing tensorial operations on these tensors. The macroscopic strain, E, can be represented by a vector of dimension 6, by using the modified Voigt notations:

E = (E 1 , ..., E 6 ) = (E 11 , E 22 , E 33 , √ 2E 23 , √ 2E 13 , √ 2E 12 ) (24)
Since the gradient of strain is symmetric according to its first two indices

(E ij,k = E ji,k ), it is possible to use a similar representation for E ij,1 , E ij,2 and E ij,3 : (E 1,k , ..., E 6,k ) = (E 11,k , E 22,k , E 33,k , √ 2E 23,k , √ 2E 13,k , √ 2E 12,k ) (25)
with k = 1, 2, 3. Thus, the gradient of strain is represented by a vector of dimension 18:

∇E(x) = (E 1,1 , ..., E 6,1 , E 1,2 , ..., E 6,2 , E 1,3 , ..., E 6,3 ) (26)
The elasticity tensors C 0,0 , C 0,1 , C 1,1 , are represented by matrices of dimensions 6 × 6, 6 × 18 and 18 × 18:

C 0,0 =        C 11 . . . C 16 . . . . . . Sym C 66        6×6 ( 27
)
C 0,1 =        C 111 . . . C
C 1,1 =                               C 1111 . . . C (28)
C 6363                               18×18 (29)
Here we propose to restrict our analysis to these tensors; the identification of higher-order elasticity tensors requires the computation of functions X i (y) for i = 3, 4, ... In the above matricial representation, coefficients C ij , C ijk and C ijkl are computed from components C 0,0 ijpq , C 0,1 ijpqr and C 1,1 ijkpqr as follows:

C IP = αC 0,0 ijpq , C IP r = αC 0,1 ijpqr , C IkP r = αC 1,1 ijkpqr ( 30
)
where indices I, J and coefficient α are defined by:

I =      i if i = j 9 -i -j if i ̸ = j P =      p if p = q 9 -p -q if p ̸ = q α =              1 if i = j and p = q √ 2 if i = j or p = q 2 if i ̸ = j and p ̸ = q (31)
The coefficient α takes the values α = 1, α = √ 2 or α = 2 depending on the values of the set {I, P }. The components of tensor C 0,1 are null for symmetry reasons. The components of tensor C 0,0 remain unchanged from the standard homogenization and are given in Boutin 1996. Components of tensor C 1,1 are computed from [START_REF] Engelen | An evaluation of higher-order plasticity theories for predicting size effects and localization[END_REF] together with equation [START_REF] Bouyge | A micromechanically based couple-stress model of an elastic orthotropic twophase composite[END_REF]. The non null components of C 1,1 are:

C 2222 = C 3333 = (1 -τ ) 2 τ 2 Q 2 4 Q 2 3Q 2 1 C 1212 = C 1313 = (1 -τ ) 2 τ 2 (λ 2 µ 1 -λ 1 µ 2 ) 2 Q 2 3(λ 2 + 2µ 2 ) 2 (λ 1 + 2µ 1 ) 2 Q 2 1 C 2323 = C 3232 = (1 -τ ) 2 τ 2 (λ 2 -λ 1 ) 2 µ 2 1 µ 2 2 Q 3 2 3(λ 2 + 2µ 2 ) 2 (λ 1 + 2µ 1 ) 2 Q 2 1 C 5252 = C 6363 = (1 -τ ) 2 τ 2 (µ 1 -µ 2 ) 2 Q 3 6µ 2 1 µ 2 2 Q 2 2 C 6262 = C 5353 = (1 -τ ) 2 τ 2 2(µ 1 -µ 2 ) 2 Q 5 3µ 2 1 µ 2 2 Q 2 2 C 4242 = C 4343 = (1 -τ ) 2 τ 2 (µ 1 -µ 2 ) 2 Q 2 6 C 1222 = C 1333 = (1 -τ ) 2 τ 2 (λ 2 µ 1 -λ 1 µ 2 )Q 4 Q 2 3(λ 1 + 2µ 1 )(λ 2 + 2µ 2 )Q 2 1 C 1232 = C 1323 = (1 -τ ) 2 τ 2 µ 1 µ 2 (λ 2 -λ 1 )(λ 2 µ 1 -λ 1 µ 2 )Q 2 2 3(λ 1 + 2µ 1 ) 2 (λ 2 + 2µ 2 ) 2 Q 2 1 C 2232 = C 3323 = (1 -τ ) 2 τ 2 µ 1 µ 2 (λ 2 -λ 1 )Q 4 Q 2 2 3(λ 1 + 2µ 1 )(λ 2 + 2µ 2 )Q 2 1 C 2243 = C 3342 = (1 -τ ) 2 τ 2 √ 2(µ 2 -µ 1 )Q 4 Q 2 6Q 1 C 1243 = C 1342 = (1 -τ ) 2 τ 2 √ 2(λ 2 µ 1 -λ 1 µ 2 )(µ 2 -µ 1 )Q 2 6(λ 1 + 2µ 1 )(λ 2 + 2µ 2 )Q 1 C 3243 = C 2342 = (1 -τ ) 2 τ 2 √ 2µ 1 µ 2 (µ 1 -µ 2 )(λ 1 -λ 2 )Q 2 2 6(λ 1 + 2µ 1 )(λ 2 + 2µ 2 )Q 1 C 5263 = C 6352 = (1 -τ ) 2 τ 2 (µ 1 -µ 2 ) 2 Q 3 6µ 2 1 µ 2 2 Q 2 2 C 5362 = C 6253 = (1 -τ ) 2 τ 2 (µ 1 -µ 2 ) 2 Q 6 3µ 2 1 µ 2 2 Q 2 2 (32)
with:

Q 1 = 1 -τ λ 1 + 2µ 1 + τ λ 2 + 2µ 2 Q 2 = 1 -τ µ 1 + τ µ 2 Q 3 = (1 -τ )µ 1 + τ µ 2 Q 4 = (µ 2 -µ 1 )Q 1 + (λ 2 -λ 1 ) µ 1 µ 2 Q 2 (λ 1 + 2µ 1 )(λ 2 + 2µ 2 ) Q 5 = µ 2 (µ 2 + λ 2 ) λ 2 + 2µ 2 τ + µ 1 (λ 1 + µ 1 ) λ 1 + 2µ 1 (1 -τ ) Q 6 = µ 2 λ 2 (λ 2 + 2µ 2 ) τ + µ 1 λ 1 λ 1 + 2µ 1 (1 -τ ) (33)
Figure 2 shows the variations of the components of the gradient elasticity tensor, given by ( 32), as functions of the volume fraction τ for materials whose elastic properties are given below. For brevity, only the "in plane" components have been plotted. It can be observed that all these higher order elastic coefficients become null when τ = 0 and τ = 1. In these two particular cases, the unit cell is constituted of a homogeneous medium, having the elastic properties λ 1 , µ 1 or λ 2 , µ 2 . At the macroscopic scale, the strain gradient effects vanish and the constitutive relations are of Cauchy type with the elastic coefficients λ 1 , µ 1 or λ 2 , µ 2 . For the applications proposed in figures 2 the following Poisson ratios have been used: ν 1 = 1/3 and ν 2 = 1/4. The Young's moduli are given by E 1 = 10GP a and E 2 = 20GP a. The gradient elastic coefficients represented on figures 2 have been normalized by the Young's modulus E 1 . Note that in [START_REF] Bouyge | A micromechanically based couple-stress model of an elastic two-phase composite[END_REF][START_REF] Kouznetsova | Multi-scale constitutive modeling of heterogeneous materials with gradient enhanced computational homogenization scheme[END_REF][START_REF] Yuan | A micromechanical approach of nonlocal modeling for media with periodic microstructures[END_REF], where higher order boundary conditions for the unit cell are used, the macroscopic behavior of the composite remains described by a gradient elasticity model even if at the local scale the material remains homogeneous. It is noteworthy that, at the contrary, the asymptotic expansion method combined with the energy based micro-macro transition described in the previous section leads to results for the gradient elastic properties which appear to be more physically realistic.

A FFT-based computational algorithm for computing the components of the six th -order elasticity tensor

We aim now at deriving a FFT based numerical method for computing the effective elasticity tensors which enter the macroscopic law [START_REF] Exadaktylos | Gradient elasticity with surface energy: Mode-I crack problem[END_REF]. As explained in section 3.1, one has to compute functions X 1 (y), X 2 (y), X 3 (y), ... by solving the hierarchy of elasticity problems provided in section A. To this end, and along the lines of the paper of Moulinec and Suquet 1994, we first propose to derive the Lippmann-Schwinger integral equation associated to the higherorder homogenization problems, which is the basis of the FFT-based iterative scheme. Indeed, the integral equation described in Moulinec and Suquet 1994 is only applicable for the 0-order homogenization problem (corresponding to g 0 (y) = 0 and p 0 (y) = c(y) : E) and must be generalized for handling problems dealing with arbitrary expressions for the body force g(y) and the polarization p(y) (here and in the next of the paper the variable x is omitted for simplicity).

The Lippmann-Schwinger equation associated to higher-order homogenization problems

Let us introduce in (A.1) a homogeneous reference medium with the elasticity tensor c 0 . The stress-strain relation can be rewritten as:

σ(y) = c 0 : ε(y) + τ (y) (34)
Where τ (y) is defined as:

τ (y) = (c(y) -c 0 ) : ε(y) + p(y) (35)
Taking now the Fourier transform, defined, for any Y-periodic function F (y), by:

F (ξ) = F(F (y)) = 1 V ∫ V F (x) exp(iξ.x)dx (36)
the elasticity problem reads:

             iσ(ξ).ξ + g(ξ) = 0 σ(ξ) = c 0 : ε(ξ) + τ (ξ) ε(ξ) = i 2 (ξ ⊗ u(ξ) + u(ξ) ⊗ ξ) (37)
Let us now introduce θ(ξ) defined by:

θ(ξ) = i ∥ξ∥ 4 [ ξ ⊗ ξg(ξ).ξ -(g(ξ) ⊗ ξ + ξ ⊗ g(ξ))∥ξ∥ 2 ] (38)
and such that g(ξ).ξ = iθ(ξ). The solution is:

ε(ξ) = -Γ 0 (ξ) : (τ (ξ) + θ(ξ)) (39)
Or equivalently, in the real space:

ε(y) = -Γ 0 (y) * (τ (y) + θ(y)) ( 40
)
where Γ 0 (y) represents the Green and the symbol " * " denotes the convolution product. Finally, replacing τ (ξ) by its expression given by ( 35), the operator for the strain and solution of (A.1) complies with the following integral equation:

ε(y) = -Γ 0 (y) * [ (c(y) -c 0 ) : ε(y) + p(y) + θ(y) ] (41)

The FFT based algorithm

Similarly to Moulinec and Suquet 1994, we propose to search the solution of Lippmann-Schwinger equation ( 41) by means of the following iterative scheme:

ε i+1 (ξ) = -Γ 0 (ξ) * [ (c(y) -c 0 ) : ε i (y) + p(y) + θ(y) ] (42)
the scheme starting from ε 1 = 0. A simplification of the above recurrence relation is possible. Since the tensors ε i (y) is compatible and < ε i (y) > V = 0 at each iteration i, it follows that:

Γ 0 (ξ) * [ c 0 : ε i (y) ] = ε i (y) (43)
Consequently:

ε i+1 (y) = ε i (y) -Γ 0 (ξ) * [ c(y) : ε i (y) + p(y) + θ(y) ] (44)
The following iterative scheme is then employed:

ε i (y) = F -1 (ε i (ξ)) ε i (y) = c(y) : ε i (y) + p(y) σ i (ξ) = F(σ i (y)) convergence test ε i+1 (ξ) = ε i (ξ) -Γ 0 (ξ) : (σ i (ξ) + θ(ξ)) (45)
The exact Fourier transform is thereafter approximated by the finite discrete Fourier transform which is computed by using the FFT algorithm. The iterative scheme is stopped when the stress field complies with the local equilibrium. The following convergence test has been used in our computations:

∥P (ξ) : σ i (ξ) + θ(ξ)∥ ∥σ i (ξ)∥ < h ( 46
)
where ∥ • ∥ denotes the Frobenius norm and P (ξ) is the fourth order tensor defined by

P (ξ) = I -Q(ξ) and Q ijpq (ξ) = (k ip k jq + k iq k jp )/2 and k ij = δ ij -ξ i ξ j /∥ξ∥ 2 .
Tensors P (ξ) and Q(ξ) are two projectors; for any secondorder tensor a, Q(ξ) : a is interpreted as the projection of a along the plane orthogonal to the wave-vector ξ whereas P (ξ) : a is the out of plane projection of a. Following Moulinec and Suquet 1994, the reference medium, giving the convergence of the FFT-based iterative scheme, is chosen as:

c 0 = 1 2 [ max y c(y) + min y c(y) ] (47)

Illustration of the method

As an illustration purpose we compute the higher-order elastic coefficients in the case of a composite made up of long fibers with a circular cross section. The sections of the fibers are assumed to be arranged along a periodic squared lattice (see figure 3). For the calculations, a squared unit cell has been considered and a representation with a grid containing 128 × 128 points has been used for effecting the Fourier transform and its inverse. The width of the cell is denoted by h and the radius of the cylinder is denoted by a. The elastic moduli of the matrix are λ 1 , µ 1 whereas those of the fibers are denoted by λ 2 , µ 2 . For all applications, the values of the elastic moduli are the ones already used in the previous section. Plane strain conditions are considered; it follows that only the components E ij,k with i, j, k = 1, 2 are used. Due to the symmetry of the unit cell, the relation giving the components of the hyperstress T as functions of the components of the macroscopic gradient of strain ∇E is:

                   T 11 T 21 T 62 T 22 T 12 T 61                    = h 2                    C 1111 C 1121 C 1162 0 0 0 C 1121 C 2121 C 2162 0 0 0 C 1162 C 2162 C 6262 0 0 0 0 0 0 C 1111 C 1121 C 1162 0 0 0 C 1121 C 2121 C 2162 0 0 0 C 1162 C 2162 C 6262                                       E 1,1 E 2,1 E 6,2 E 2,2 E 1,2 E 6,1                    (48)
where the compact notation for the components of the hyperstress and straingradient are used. The fourth-order elasticity tensor can be obtained from usual homogenization procedures and its coefficients are not reported here. Relation (48) involves only 6 elastic coefficients. The variation of these coefficients with ratio a/h are shown on figure 4. It can be observed that the gradient elastic coefficients are again null when a = 0, corresponding to the case of a unit cell having the constant elastic moduli λ 1 and µ 1 ., whereas these coefficients are maximal when the ratio a/h is maximal and equal to 1/2, the composite being still heterogeneous in this last case. The changes with concentration of theses coefficients are monotonous for C 1111 and C 1121 while they are characterized by a minimum between a/h equal to 0.3 and 0.5 for all other coefficients.

Conclusion

In this paper we have proposed a simple and efficient method for the determination of the gradient elastic properties of periodic composite materials. The approach can be summarized as follows: (i) the elastic solution is expanded along a series expansion for which higher order terms are kept (in addition to those generally considered in the usual homogenization framework) (ii) the elastic energy of the equivalent homogeneous material is computed by means of a generalized Hill-Mandel lemma, (iii) the constitutive relations are derived from the "state law" associated with this macroscopic potential. The averaging rule giving the higher-order elasticity tensors have been derived. The components of these higher-order elasticity tensors are obtained by solving higher-order homogenization problems for the unit cell of the periodic medium. Closed-form solutions for these higher-order elasticity tensors are provided in the particular case of a stratified composite. For arbitrary microstructures, a computational method has been developed in the last part of the paper. The method of resolution uses an iterative scheme and the exact expression of the Green's operator in the Fourier space. The relevance of this approach has been illustrated in the case of a periodic array of elastic fibers embedded in an elastic matrix. This work opens the way to many new applications and extensions not only in the domain of linear but also in the domain of non linear homogenization. For instance, the determination of gradient properties of elastic or plastic polycrystalline materials is an array of research of high importance. The extension of our approach in this context will be the scope of a future work.

B Solution for a stratified composite

For simplicity, the variable y 1 is denoted thereafter by y. The solution at the first order involves the components of the third order tensor X (1) (y) which are: in which H(ψ) and D(ψ) are given, for any ψ, by:

H(ψ) = [ 1 -τ ψ a + τ ψ b ] -1 D(ψ) = τ (1 -τ )(ψ a -ψ b ) (B.4)
Where ψ a and ψ b are respectively the values of ψ in layer a and in layer b. The solution at the second order involves the components of X 2 (y) which are:

Fig. 1 .

 1 Fig. 1. The periodic unit cell for a stratified composite.

Fig. 2 .

 2 Fig. 2. Variations of the component C 2222 , C 1212 , C 6262 and C 1222 of the six th -order elasticity tensor for the stratified medium with the volume fraction τ .

Fig. 3 .

 3 Fig. 3. Regular array of cylindrical inclusions aligned along direction Ox 3 .

Fig. 4 .

 4 Fig. 4. Variations of the component C 1121 , C 2121 , C 1121 , C 6262 , C 1162 , C 2162 of the six th -order gradient elasticity tensor as a function of a/h for the fiber-reinforced composite.

X 1 111

 1 (y) = kf (y) X 1 122 (y) = X 1 133 (y) = k[D(λ)/D(λ + 2µ)]f (y) X 1 212 (y) = X 1 313 (y) = mf (y) (B.1)where function f (y) is given by: h-(1 -τ )/2]/(1 -τ ) in layer a -[y/h + τ /2]/τ in layer b (B.2)coefficients k and m are given by:k = H(λ + 2µ)D(1/(λ + 2µ)) m = H(µ)D(1/µ) (B.3)

  161 C 112 . . . C 162 C 113 . . . C 163

										
	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .	     
										6×18

C 611 . . . C 661 C 612 . . . C 662 C 613 . . . C 663

  1161 C 1112 . . . C 1162 C 1113 . . . C 1163 C 6161 C 6112 . . . C 6162 C 6113 . . . C 6163

	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .
			C 1212 . . . C 1262 C 1213 . . . C 1263
				. . .	. . .	. . .	. . .	. . .
					C 6262 C 6213 . . . C 6263
		Sym				C 1313 . . . C 1363
							. . .	. . .

A Hierarchy of problems

The hierarchy of linear elastic problems, leading to the identification of tensors X α (y) for α = 1, 2, 3, ..., can be put into the form:

In which the coordinate x is omitted for simplicity. Expressions of g α (y) and p α (y) are: order 0 :

order 2 :

Where β(y) is given by:

When the body force is of the type: f (x, y) = ρ(y)f 0 (x) , ρ(y) being the density of the material constituents. The solutions of these linear elastic problems for α = 0, 1, 2, etc. are:

where tensors X 1 (y), X 2 (y) X 3 (y) are localization tensors which do not depend on macroscopic variables E ,∇E and ∇ 2 E, but depend only on the microstructure, as in the classical homogenization framework.

where function F (y) is defined by: