Jeremy Bleyer
email: jeremy.bleyer@enpc.fr

Multiphase continuum models for fiber-reinforced materials

Keywords: fiber-reinforced materials, homogenization, generalized continuum, scale effects, crack bridging

This contribution addresses the formulation of a generalized continuum model called multiphase model aimed at describing more accurately the mechanical behavior of fiber-reinforced materials. Improving on the classical macroscopic description of heterogeneous materials by an effective homogeneous Cauchy medium, such models rely on the superposition of several continua (or phases) possessing their own kinematics at the macroscopic level and being in mutual interaction (in the same spirit of deformable porous media). Up to now, they have only been formulated based on phenomenological assumptions and the identification of the corresponding constitutive parameters remained unclear. The aim of this paper is three-fold. First, a homogenization procedure is described, enabling to derive constitutive parameters from the resolution of a generalized auxiliary problem on a classical heterogeneous microstructure. Second, analytical and numerical derivation of these properties is performed in various cases. Finally, illustrative applications on boundary-value problems assess the validity of the homogenization procedure and illustrate the relevance of such generalized models which are able to capture scale effects and to model crack bridging and delaminated configurations at the macroscopic level. It is also shown to encompass results of shear-lag models for analyzing stress transfers in fiber/matrix composites.

Introduction

Fiber-reinforced materials, consisting of short or long linear inclusions embedded in a surrounding matrix, offer extremely interesting mechanical performances compared to more conventional materials. Optimal properties are achieved through careful design of many parameters such as fiber volume fraction, elastic modulus contrast, fiber orientation, fiber aspect ratio, ductility or brittleness of each phase, quality of the interface, and many more.

Obviously, this offers a large range of mechanical behaviors which must be well understood for each particular industrial application. Fiber-reinforced materials are now present across many scales, ranging for instance from carbon nanotubes to fiber-reinforced concretes or even piled-raft foundations in civil engineering structures. Besides, environmental aspects also drive the development of bio-composites (natural fibers or bio-derived matrix) and many biological materials such as bones or nacre are based on a complex heterogeneous and hierarchical structure from which interesting insights can be gained.

For engineering applications, the composite properties are generally described in an effective macroscopic way. Homogenization theory is a formidable tool to perform this up-scaling procedure [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF][START_REF] Sánchez-Palencia | Non-homogeneous media and vibration theory[END_REF] and has been used to establish classical results for the specific case of reinforced solids [START_REF] Hashin | Analysis of Properties of Fiber Composites With Anisotropic Constituents[END_REF][START_REF] Hashin | The elastic moduli of fiber-reinforced materials[END_REF][START_REF] Hill | Elastic properties of reinforced solids: some theoretical principles[END_REF]. In this classical framework, the macroscopic effective behavior is that of a standard Cauchy continuum associated with a macroscopic stiffness tensor C hom . An important result of a two-scale asymptotic expansion with respect to the scale factor is that the lowest-order displacement u 0 is a rigid-body motion at the microscopic scale, depending upon the macroscopic space variable only. The effective constitutive relation then involves the macroscopic stiffness C hom obtained as the solution of an auxiliary (or corrector) problem defined on the heterogeneous unit cell or RVE.

However, this simple picture fails to describe many situations, as pointed out for instance in [START_REF] Sanchez-Palencia | Homogenization in mechanics. A survey of solved and open problems[END_REF] and associated references. First, the previous homogenized behavior pertains to the limit case in which → 0. In practice, heterogeneous structures correspond to fixed but nonzero values of and therefore exhibit so-called scale effects when is not sufficiently small i.e. when scale separation is moderate. Second, a well-known difficulty in deriving higherorder models able to capture such scale effects is related to what happens near the boundary of the heterogeneous structure. The presence of heterogeneities creates the appearance of boundary layers [START_REF] Dumontet | Study of a boundary layer problem in elastic composite materials[END_REF][START_REF] Sanchez-Palencia | Homogenization in mechanics. A survey of solved and open problems[END_REF] which do not affect the zero-order terms far from the boundary but have a non-negligible impact on higher-order terms. Thirdly, the previous asymptotic expansion assumes no other "small" parameter than the scale factor , in particular, it does assume that heterogeneous moduli do not scale with . In the presence of large contrasts of stiffness properties, e.g. C1 = O( 2 )C 2 for a biphasic material as considered in [START_REF] Boutin | Generalized inner bending continua for linear fiber reinforced materials[END_REF], the zero-order asymptotic behavior differs from the classical one. In this case, u 0 is not a rigid-body motion in the whole unit cell anymore but only in the rigid phase, the remaining part exhibiting a relative displacement which depends on the microscopic variable. The resulting effective behavior then becomes that of a two-phase medium in which the relative displacement satisfies a Darcy-like law. In the limit case of regions with completely debonded phases, the effective behavior has also been shown to be one of a two-phase medium including two independent macroscopic displacement fields corresponding to the individual motion of the matrix and of the debonded fibers [START_REF] Berrehili | The homogenized behavior of unidirectional fiber-reinforced composite materials in the case of debonded fibers[END_REF]. Finally, slenderness of heterogeneities also introduce an additional small parameter for which higher-order terms may contribute to the macroscopic behavior [START_REF] Boutin | Generalized inner bending continua for linear fiber reinforced materials[END_REF]. Unfortunately, the problematics of poor scale separation 1 , boundary layer effects, large stiffness contrasts, weak interfaces and geometrical slenderness are all frequently encountered in fiber-reinforced materials.

Treatment of such situations may rely on different strategies such as multi-scale computations or investigation of higher-order terms in the previous asymptotic expansions leading to generalized continua [START_REF] Boutin | Microstructural effects in elastic composites[END_REF]. One main interest of the latter is to introduce a length-scale dependence at the macroscopic level in order to capture scale or size effects characterizing the specific heterogeneous microstructure. Generalized continua can involve either additional degrees of freedom (Cosserat, micromorphic, stress-gradient, etc.) or higher-order derivatives (strain gradient or constrained micromorphic models) [START_REF] Altenbach | Generalized Continua-from the Theory to Engineering Applications[END_REF][START_REF] Forest | Homogenization methods and the mechanics of generalized continua[END_REF]. The wide range of possible theory makes it extremely difficult to assess their precise individual domain of validity or relevant field of application. Since all theories introduce, at least indirectly, an internal length scale, different size effects can be obtained. Regarding the specific case of fiber-reinforced materials, this aspect can be quite complex to assess as evidenced by the simple example considered in [START_REF] De Buhan | Elastic, Plastic and Yield Design of Reinforced Structures[END_REF][START_REF] De Buhan | Multiphase approach as a generalized homogenization procedure for modelling the macroscopic behaviour of soils reinforced by linear inclusions[END_REF] of a multilayered 2D domain exhibiting stiffening size effects (i.e. when the structure scale L decreases for a fixed microstructure length s) when loaded in shear and softening size effects when loaded in tension/compression in the transverse direction. As a result, since strain gradient theories are known to lead to stiffening size effects whereas other models (e.g. stress gradient) lead to softening size effects [START_REF] Tran | Modélisation à plusieurs échelles d'un milieu continu hétérogène aléatoire[END_REF], choosing a specific generalized continuum model for fiber-reinforced media is not obvious and depends on the precise application.

Mathematical asymptotic analyses have been conducted in the specific case of a linear inclusion embedded in a matrix, either in a (scalar) conductivity framework [START_REF] Bellieud | Homogenization of elliptic problems in a fiber reinforced structure. Non local effects[END_REF][START_REF] Briane | Nonlocal effects in two-dimensional conductivity[END_REF][START_REF] Cherednichenko | Non-local homogenized limits for composite media with highly anisotropic periodic fibres[END_REF] or for three-dimensional elasticity [START_REF] Bellieud | Homogenization of an elastic material reinforced by very stiff or heavy fibers. Non-local effects. Memory effects[END_REF][START_REF] Briane | Homogenization of two-dimensional elasticity problems with very stiff coefficients[END_REF][START_REF] Pideri | A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium[END_REF]. The main result of these works is that, depending on the precise scaling of the fiber/inclusion material contrast (conductivity or elastic modulus) and of the fiber volume fraction with respect to the scale factor , different types of homogenized behaviors can be obtained. In particular, in the case of elastic fiberreinforced composites 2 , these models can be:

• a standard single-phase Cauchy continuum • a single-phase continuum including second-gradient effects (due to the bending of the fibers)

• a two-phase continuum without bending effects with one displacement variable for each phase (fiber and matrix)

• a two-phase continuum including bending effects Similar results involving effective two-phase continua have also been obtained in the context of porous media with double porosity [START_REF] Arbogast | Derivation of the double porosity model of single phase flow via homogenization theory[END_REF] or for an elastic porous solid filled with a composites, it may not be the case when considering their behavior around small defects or cracks for instance. Besides, scale separation may sometimes be poor in civil engineering reinforced structures such as piled-raft foundations. 2 For scalar conductivity problems, no second-gradient effect but only single or two-phase continua can be obtained as mentioned in [START_REF] Pideri | A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium[END_REF]. low-viscosity fluid [START_REF] Sanchez-Palencia | Homogenization in mechanics. A survey of solved and open problems[END_REF]. Experimental evidences of such generalized behaviors in a dynamic setting can be found for instance in [START_REF] Boutin | Experimental evidence of the highgradient behaviour of fiber reinforced materials[END_REF][START_REF] Hegemier | A continuum mixture theory of wave propagation in laminated and fiber reinforced composites[END_REF]. These results provide an explanation of the different scale effects mentioned earlier for the multilayered domain. The shear-loading situation induces bending of the fibers, the energy of which has a non negligible contribution for slender and stiff inclusions and which can be described by a continuum involving second-gradient effects. Conversely, the softening size effect triggers a relative displacement between the fiber and the matrix near the lateral free boundary which can be accounted by a generalized continuum considering two separate displacements for each phase.

Adopting a multi-component or a multi-phase modeling of a continuum system has already been proposed in the context of the theory of mixtures, see for instance [START_REF] Atkin | Continuum theories of mixtures: basic theory and historical development[END_REF][START_REF] Bedford | Theories of immiscible and structured mixtures[END_REF] for a review. Such a point of view can apply to fluid mixtures, solid and fluid for poroelastic media [START_REF] Bowen | Incompressible porous media models by use of the theory of mixtures[END_REF][START_REF] Coussy | From mixture theory to Biot's approach for porous media[END_REF] or two solids for heterogeneous materials such as laminated composites [START_REF] Bedford | Theories of immiscible and structured mixtures[END_REF][START_REF] Bedford | A multi-continuum theory for composite elastic materials[END_REF][START_REF] Mcniven | A mixture theory for elastic laminated composites[END_REF]. Following the same idea in the case of fiber-reinforced materials, so-called multiphase models have been developed initially in [START_REF] Sudret | Modélisation multiphasique des ouvrages renforcés par inclusions[END_REF][START_REF] Sudret | Modélisation multiphasique de matériaux renforcés par inclusions linéaires[END_REF] without bending effects and in [START_REF] De Buhan | Micropolar model for materials reinforced by linear inclusions[END_REF] with bending effects by postulating a specific constitutive equation for the matrix and the fiber phase. Extensions of these models have then concerned elastoplasticity, damage and dynamic effects, see [START_REF] De Buhan | Elastic, Plastic and Yield Design of Reinforced Structures[END_REF] for a recent review. Similar models have also been obtained in [START_REF] Boutin | Generalized inner bending continua for linear fiber reinforced materials[END_REF][START_REF] Soubestre | Non-local dynamic behavior of linear fiber reinforced materials[END_REF] based on asymptotic expansion methods. Finally, let us mention that in the dynamic case, particular dispersion relations and memory effects due to time non-locality are obtained [START_REF] Boutin | Experimental evidence of the highgradient behaviour of fiber reinforced materials[END_REF][START_REF] Nguyen | Analyse sismique des ouvrages renforcés par inclusions rigides à l'aide d'une modélisation multiphasique[END_REF][START_REF] Soubestre | Non-local dynamic behavior of linear fiber reinforced materials[END_REF].

In this work, we advocate the use of multiphase models to account for all these previously mentioned difficult aspects related to scale effects in fiber-reinforced composites. One key feature of such models is that when all phases possess the same kinematics the standard Cauchy continuum model is retrieved. It will be seen that this happens when the scale factor goes to zero. Another important additional feature is that second-order effects may be considered, for instance when one phase is described by a Cosserat continuum instead of a Cauchy continuum, to include shear and bending deformation of the reinforcements. This aspect will not be considered in the present manuscript and will constitute an important future perspective. As a result, only mutually interacting Cauchy phases will be considered, enabling to account only for higher-order effects induced by scale separation and boundary layers. However, one must keep in mind that the capabilities of these models go far beyond that.

In all the previously mentioned works on mixture theory or multiphase models, constitutive relations were often postulated or analytically derived in specific situations only. In particular, it was often assumed that the stress of one phase depends upon the strain of this particular phase only [START_REF] Bedford | Theories of immiscible and structured mixtures[END_REF][START_REF] Bedford | A multi-continuum theory for composite elastic materials[END_REF][START_REF] De Buhan | Elastic, Plastic and Yield Design of Reinforced Structures[END_REF][START_REF] Sudret | Modélisation multiphasique des ouvrages renforcés par inclusions[END_REF], an unnecessary hypothesis from which we will later depart. Further development of these models, although extremely simple in their implementation, has been therefore hindered by a lack of rationalized up-scaling schemes, even in a linear elastic setting. In this work, we aim at deriving a general micromechanical homogenization procedure linking a heterogeneous standard Cauchy medium at the microscopic scale to a multiphase continuum model at the macroscopic scale. Although illustrated here only for laminated or infinitely long fiber-reinforced media, this procedure can be applied to any microstructure (short-fiber reinforced composites for instance) and will serve as a rational basis for deriving more complex (non-linear in particular) effective behaviors in the fiber-reinforced material reinforcement phase matrix phase relative displacement Figure 1: Fiber-reinforced material (N = 1 reinforcing phase) represented by a multiphase continuum with a two-phase kinematics (U 1 , U 2 ) at any material point, interaction is described through the relative displacement V = U 2 -U 1 between both phases. future. Section 2 will thus first recall the governing equations of the multiphase model in a linear elastic setting and the homogenized procedure for deriving its effective properties will be presented in section 3. Explicit derivation of these effective properties will be discussed in section 4. The other main goal of this paper is to illustrate the modeling capabilities of such models and validate the previous homogenization procedure on some boundary-value problems, involving scale effects near free boundaries or modeling of crack bridging and partial fiber/matrix delamination. This will be the purpose of section 5. A link with shear-lag models will also be discussed in Appendix C.

Notations: Throughout this paper, the following notations will be used. Italic plain symbols correspond to scalar quantities. Vectors and second-order tensors are denoted by italic boldface symbols, e.g. u for the displacement and ε for the strain. Fourth-order tensors are denoted by blackboard bold letters, e.g. C for the elasticity tensor. Double dot product between tensors is given by a : b = a ij b ji and the transpose operator is denoted by ( )

T . The symmetric gradient operator is denoted by ∇ s and the divergence operator by div. When different spatial variables are present (e.g. x, y) a subscript will indicate upon which variable a differential operator acts, e.g. ∇ s

x or div y . Average of a quantity over a domain A will be noted as = 1 |A| A ( )dΩ whereas partial averages over a subset A i ⊆ A will be noted as

i = 1 |A i | A i ( )dΩ.

Governing equations of the multiphase model

As mentioned in the introduction, different authors proposed similar models in various frameworks, we will first start from the formulation proposed initially in [START_REF] Sudret | Modélisation multiphasique des ouvrages renforcés par inclusions[END_REF][START_REF] Sudret | Modélisation multiphasique de matériaux renforcés par inclusions linéaires[END_REF] excluding bending effects. It aims at treating the case of a composite material made of a matrix reinforced by N arrays of differently oriented linear inclusions, each family characterized by a unit vector e i with i = 2, . . . , N + 1. The set of virtual motions consists of N + 1 virtual velocity vectors: U 1 for the matrix and U i with i = 2, . . . , N + 1 for each inclusion family (Figure 1). Assuming that each inclusion family behaves in tension/compression only and postulating an appropriate form of the internal, external and inertial forces virtual work, the use of the virtual work principle enables to obtain the following equilibrium equations (quasi-static case):

div σ 1 + ρ 1 F + N +1 i=2 I i = 0 (1) div(σ i e i ⊗ e i ) + ρ i F -I i = 0 ∀i = 2, . . . , N + 1 (2)
where:

• σ 1 is the partial stress in the matrix phase and σ i the uniaxial partial stress in direction e i in each reinforcement family

• ρ 1 F (resp. ρ i F
) is an external body force exerted on the corresponding phase

• I i is an interaction force density exerted by phase i on the other phases (note that the interaction force exerted by all reinforcements on the matrix phase is equal to N +1 i=2 I i due to the action/reaction principle) These equilibrium equations are to be completed by the following traction boundary conditions:

σ 1 • n = t 1 on ∂Ω T (3)
σ i (e i • n)e i = t i on ∂Ω T ∀i = 2, . . . , N + 1 (4)
on some part ∂Ω T of the domain boundary of exterior unit normal n. It is important to highlight that the classical equilibrium and boundary conditions of a standard Cauchy medium are retrieved when summing the contribution of each phases:

div Σ + ρF = 0 (5) Σ • n = t on ∂Ω T (6)
where we introduced the total mass density ρ, the total traction t and the total stress Σ:

ρ = N +1 j=1 ρ j (7) t = N +1 j=1 t j (8)
Σ = σ 1 + N +1 i=2 σ i e i ⊗ e i (9)
The previous equilibrium equations correspond to the following deformation work density:

w def = σ 1 : ε 1 + N +1 i=2 (σ i e i ⊗ e i ) : ε i + I i • (U i -U 1 ) (10)
where ε j = ∇ s U j for j = 1, . . . , N +1 is the linearized strain tensor of each phase kinematics.

In the previous expression, it clearly appears that these quantities are work-conjugate to the partial stresses of each phase whereas the work-conjugate quantity of the interaction force I i is the relative displacement V i = U i -U 1 between the reinforcement phase i and the matrix phase. Again, it is important to mention that in the case when each phase has the same kinematics U = U 1 = U i , the deformation work density reduces to its classical expression w def = Σ : ε with ε = ε 1 = ε i for all i = 2, . . . , N + 1.

Postulating the existence of a convex free energy density ψ function of the phase strains ε j and of the relative displacements V i between phase i and phase 1, the constitutive equations of the multiphase model are derived from:

σ 1 = ∂ψ ∂ε 1 (ε 1 , ε i , V i ) ( 11
)
σ i = ∂ψ ∂ i (ε 1 , ε i , V i ) (12)
I i = ∂ψ ∂V i (ε 1 , ε i , V i ) ( 13
)
where i = e i • ε i • e i .

A few works have further postulated that the free energy density ψ can be decomposed into individual contributions of each phases, depending only upon the corresponding phase strain, and into a contribution due to the interaction between phases depending upon the relative displacements only:

ψ(ε 1 , ε i , V i ) = ψ 1 e (ε 1 ) + N +1 i=2 ψ i e (ε i ) + ψ i I (V i ) (14)
For instance, the following phenomenological linear elastic behavior is considered in [START_REF] De Buhan | Elastic, Plastic and Yield Design of Reinforced Structures[END_REF][START_REF] Sudret | Modélisation multiphasique des ouvrages renforcés par inclusions[END_REF]:

ψ 1 e (ε 1 ) = 1 2 ε 1 : C 1 : ε 1 ⇒ σ 1 = C 1 : ε 1 ( 15
)
ψ i e (ε i ) = 1 2 α i ( i ) 2 ⇒ σ i = α i i ( 16
)
ψ i I (V i ) = 1 2 V i • c i I • V i ⇒ I i = c i I • V i (17)
where C 1 is the matrix phase elasticity tensor, α i = φ i E i is the axial stiffness density per unit transverse area of the fiber family i with φ i being its volume fraction, E i its Young modulus and c i I is an interaction stiffness tensor between phase i and phase 1. Although the previous constitutive equations seem reasonable, their specific form have not been justified from a micromechanical approach. In particular, the precise value of the interaction stiffness remains unclear and needed to be identified through fitting the multiphase model solutions with finite-element simulations [START_REF] De Buhan | Multiphase approach as a generalized homogenization procedure for modelling the macroscopic behaviour of soils reinforced by linear inclusions[END_REF].

The multiphase model which has just been recalled has been initially developed in order to model the reinforcement of a given material by thin one-dimensional inclusions, the constitutive material of which is also much stiffer than the surrounding matrix in general. It will be shown in section 4.2 that the previous constitutive equations can indeed be justified in this very specific case.

In the present work, we will consider a more general formulation of the multiphase model, encompassing the previous one, in which the partial stress in the reinforcement phases will not be restricted to be uniaxial but will admit a general three-dimensional tensorial form (see [START_REF] Ben Hassine | A macroscopic model for materials reinforced by flexible membranes[END_REF] for a formulation using 2D stresses for modeling flexible membranes). The set of governing equations are given here in the case of only n = N + 1 = 2 phases:

• Generalized stresses: partial stress σ j in both phases and interaction force I (exerted by phase 2 over phase 1)

• Equilibrium equations:

div σ 1 + ρ 1 F + I = 0 (18) div σ 2 + ρ 2 F -I = 0 (19)
• Traction boundary conditions:

σ 1 • n = t 1 on ∂Ω T (20) σ 2 • n = t 2 on ∂Ω T (21)
• Kinematics: displacement vector U j in both phases

• Generalized strains: infinitesimal strain in each phase ε j = ∇ s U j and relative displacement between both phases V = U 2 -U 1

• Deformation work density:

w def = σ 1 : ε 1 + σ 2 : ε 2 + I • V (22)
• Constitutive equations based on a convex free energy density ψ:

σ 1 = ∂ψ ∂ε 1 (ε 1 , ε 2 , V ) (23)
σ 2 = ∂ψ ∂ε 2 (ε 1 , ε 2 , V ) (24)
I = ∂ψ ∂V (ε 1 , ε 2 , V ) (25)
In this context, the obtained governing equations share the same structure of those describing multiphase porous media [START_REF] Coussy | Poromechanics[END_REF][START_REF] Coussy | From mixture theory to Biot's approach for porous media[END_REF]. In particular, using a suitable averaging procedure, the balance equations of each phase at the macroscopic level can be obtained from quantities expressed at the microscopic level as follows:

div x φ f σ(y, x) f + ρ f F (x) + f int (x) = 0 (26a) div x φ s σ(y, x) s + ρ s F (x) -f int (x) = 0 (26b) with f int (x) = 1 |ω| Γ f /s σ(y, x) • n f dS (26c)
where x (resp. y) represents the macroscopic (resp. microscopic) space variable, ω is an elementary representative volume located at x, Γ f /s represents the fluid/solid interface, φ f , φ s are their relative volume fraction and n f denotes the outer unit normal of the fluid elementary volume. The generalized forces of the multiphase model therefore have a firm microscopic interpretation since the partial stresses σ i can, indeed, be interpreted as the phase average stress of constituent i whereas the multiphase interaction force I can be interpreted as the resulting contact force f int over the common interface exerted by phase 2 (solid) over phase 1 (fluid).

A generalized homogenization procedure for deriving multiphase continuum models

In this section, we present a generalized homogenization procedure aimed at deriving the material parameters characterizing the multiphase model. We will restrict here to the linear elastic setting and to the case of a biphasic material although the procedure can be easily generalized to nonlinear constitutive behaviors and n-phase materials. We will first build an auxiliary problem from a stress-based approach using the minimum complementary energy principle. However, this point-of-view is not well suited to a finite-element implementation and the associated displacement-based problem will also be derived. Effective constitutive equations are then obtained using standard micromechanics arguments and a formal twoscale asymptotic expansion linking microscopic and macroscopic quantities is provided.

An auxiliary problem using the minimum complementary energy principle

Let us consider the unit cell A of a periodic heterogeneous biphasic material, its restriction to phase i is noted A i . Both phases are assumed to obey a linear elastic behavior (with phasewise-uniform stiffness tensor C i ) and to be perfectly bonded along their common interface Γ.

We recall that the derivation of standard effective stiffness properties C hom can be performed from the resolution of the following auxiliary problem with periodic boundary con-ditions: Find (u, σ) such that

div σ = 0 in A (27a) σ(y) = C i : ε(y) in A i (27b) ε(y) = ∇ s U (y) = E + ∇ s u(y) in A (27c) U (y) = E • y + u(y) in A (27d) σ • n A-antiperiodic (27e) u(y) A-periodic (27f)
where E is the given macroscopic strain and in which the first equation is to be understood in the sense of distributions (i.e. with σ • n continuous).

The solution to the previous problem can be equivalently obtained by minimizing the complementary energy of the unit cell as follows:

Ψ * (Σ) = min σ 1 |A| A 1 2 σ : C -1 (y) : σ dΩ s.t. div σ = 0 in A σ • n A-antiperiodic σ = Σ (28)
where Σ is a controlled loading parameter of dimension 6. Optimality conditions of this constrained optimization problem can indeed be shown to yield equations ( 27), the macroscopic strain E appearing, in particular, as the Lagrange multiplier associated with the last constraint. Due to the problem linearity, the corresponding optimal objective function is the macroscopic complementary energy and admits the following quadratic form:

Ψ * (Σ) = 1 2 Σ : (C hom ) -1 : Σ.
As regards the multiphase model effective properties, we propose the following generalization of the previous optimization problem. Instead of controlling only the average stress Σ, we will control independently the partial stress in each phase using a loading parameter of dimension 12 characterized by Σ 1 and Σ 2 . In addition, we will also independently control the average traction exerted by phase 2 on phase 1 over the interface Γ using a three-dimensional parameter I. However, if σ(y) is divergence-free at the local scale, then I must necessary be zero due to the divergence theorem. As a result, a fictitious body force needs to be introduced in both phases in order to control the average traction. This body force is chosen as:

b(y) = ζ(y)I where ζ(y) = -1/φ 1 in phase 1 1/φ 2 in phase 2 (29)
It has the important property of being self-equilibrated in average, that is b = 0. The proposed problem therefore reads as:

Ψ * (Σ 1 , Σ 2 , I) = min σ 1 |A| A 1 2 σ : C -1 (y) : σ dΩ s.t. div σ + ζ(y)I = 0 in A σ • n A-antiperiodic φ 1 σ 1 = Σ 1 φ 2 σ 2 = Σ 2 (30)
Using the divergence theorem on both phases and the antiperiodicity on the unit cell boundary, the local equilibrium equation necessarily induces that:

I = 1 |A| Γ σ(y) • n 1→2 dS (31)
so that the constraint linking I to the average traction [START_REF] Dumontet | Study of a boundary layer problem in elastic composite materials[END_REF] does not need to be included explicitly in the previous list of constraints. As before, the optimal objective function corresponds to the 15-dimensional macroscopic complementary energy of the multiphase model.

Owing to the linearity of [START_REF] Dormieux | Microporomechanics[END_REF], it can already be anticipated that it will depend quadratically upon its parameters (Σ 1 , Σ 2 , I).

Although we do not pretend to give a rigorous justification of this auxiliary problem, this choice is motivated by the fact that the controlled macroscopic loading parameters are precisely the generalized stresses of the multiphase model. In addition, due to the partial balance equations of each phase (26) (ignoring macroscopic body forces), the chosen link with microscopic quantities is such that the multiphase model equilibrium equations will be necessarily satisfied at the macroscopic level. More details regarding this matter will be discussed in section 3.6.

Associated displacement-based problem

Although ( 30) is a well-posed problem, its stress-based form is not well-suited to a finiteelement implementation. A displacement-based formulation can be derived on the basis of the corresponding optimality conditions. Introducing the following Lagrangian:

L(σ, u, v, E i ) = 1 |A| A 1 2 σ : C -1 (y) : σ dΩ + 1 |A| A (div σ + ζ(y)I) • u dΩ (32) + 1 |A| ∂A + (t + + t -) • v dS + 2 i=1 (Σ i -φ i σ i ) : E i
where the Lagrange multipliers associated with the constraints of (30) are respectively u, which is defined on the whole unit cell A, v, defined on half of the unit cell boundary ∂A + , and E i which are constant tensors. In the Lagrangian, the antiperiodicity constraint has been written by defining t + (y

+ ) = σ • n(y + ) with y + ∈ ∂A + and t -(y + ) = σ • n(y -)
where y -∈ ∂A -is the corresponding point with respect to y + on the opposite boundary.

The previous Lagrangian depends implicitly on the fixed parameters Σ i and I. Using the divergence theorem and after rearranging some terms, the Lagrangian can also be expressed as:

L(σ, u, v, E i ) = 1 |A| A 1 2 σ : C -1 (y) : σ -σ : (∇ s u + χ 1 E 1 + χ 2 E 2 ) dΩ (33)
+ 1 |A| ∂A + (t + + t -) • v dS + ∂A (σ • n) • u dS + 2 i=1 Σ i : E i + I • ( u 2 -u 1 )
where χ i are the characteristic functions of phase i.

The primal-dual solution (σ, u, v, E i ) to the optimization problem ( 30) is given by the following saddle-point problem:

L(σ, u, v, E i ) = max u ,v ,E i L(σ, u , v , E i ) = max u ,v ,E i min σ L(σ , u , v , E i ) (34) = min σ max u ,v ,E i L(σ , u , v , E i ) = min σ L(σ , u, v, E i )
Its optimality conditions with respect to σ are thus given by:

C -1 (y) : σ -(∇ s u + χ 1 E 1 + χ 2 E 2 ) = 0 in A (35) u(y + ) = u(y -) on ∂A (36)
which therefore results in the following boundary-value auxiliary problem: Find (u, σ) such that

div σ + ζ(y)I = 0 in A (37a) σ(y) = C i : (E i + ∇ s u) in A i (37b) σ • n A-antiperiodic (37c) u(y) A-periodic (37d)
In this problem, the controlled loading parameters are therefore (E 1 , E 2 , I). Compared to the standard auxiliary problem ( 27), (37) admits a similar form but differs by the presence of the self-equilibrated body forces and the local strain in each phase being given by

ε el = E i + ∇ s u instead of ε = E + ∇ s u.
Although there is no mention of any macroscopic displacement U (y), it can however be rewritten as: with

div σ + ζ(y)I = 0 in A (38a) σ(y) = C i : (ε(y) -(E -E i )) in A i (38b) U (y) = E • y + u(y) in A (38c) σ • n A-antiperiodic (38d) u(y) A-periodic (38e)
E = φ 1 E 1 + φ 2 E 2
and where η i = E -E i can be seen as an eigenstrain of zero average over the unit cell: η = 0. Hence, instead of controlling problem ( 37) with E 1 and E 2 , one can equivalently control problem [START_REF] Hegemier | A continuum mixture theory of wave propagation in laminated and fiber reinforced composites[END_REF] with the macroscopic strain E and some self-

equilibrated eigenstrain ∆E = φ 1 η 1 = -φ 2 η 2 = φ 1 φ 2 (E 2 -E 1
) as independent parameters (Figure 2). If E 1 = E 2 (or ∆E = 0) and I = 0, one obviously recovers [START_REF] De Buhan | Elastic, Plastic and Yield Design of Reinforced Structures[END_REF].

Let us point out that, when E i = 0, the previous problems are reminiscent of those considered for deriving macroscopic permeabilities when homogenizing a Stokes flow at the pore scale towards Darcy equation at the macroscopic scale [START_REF] Auriault | Etude du comportement macroscopique d'un milieu poreux saturé déformable[END_REF][START_REF] Bignonnet | Caractérisation expérimentale et modélisation micro-mécanique de la perméabilité et la résistance de roches argileuses[END_REF][START_REF] Bignonnet | Upper bounds on the permeability of random porous media[END_REF][START_REF] Boutin | Study of permeability by periodic and self-consistent homogenisation[END_REF]. In such a situation, the forcing term of the auxiliary problem is a uniform pressure gradient acting upon the fluid phase while the solid phase is assumed to rigid. In order to extend the definition of stresses into the solid phase, an equilibrated opposite body force is also considered inside the solid part. Our considered auxiliary problem is thus very similar, at least formally, with the difference that both phases can deform in response to this self-equilibrated body force.

Finally, homogenization procedures for generalized continua usually require establishing a macro-homogeneity condition related to a generalized Hill-Mandel lemma. In our case, the presence of fictitious body forces somehow complicates the analysis since they appear as external forces in the potential energy of the unit cell which is to be related with the macroscopic virtual work of deformation. At this macroscopic scale, these body forces do not appear as external forces but as generalized internal forces. Upon defining the set SA(I) of statically admissible stress fields and the set KA of kinematically admissible periodic fluctuations as follows:

SA(I) =    σ s.t. div σ + ζ(y)I = 0 in A σ • n continuous in A σ • n antiperiodic on ∂A    (39)
KA = {u continuous, piecewise differentiable in A and periodic on ∂A} [START_REF] Logg | Automated solution of differential equations by the finite element method: The FEniCS book[END_REF] and considering any stress field σ ∈ SA(I) and any periodic fluctuation u ∈ KA, one easily shows that:

σ : ε el = σ : (E + ∇ s u ) = φ 1 σ 1 : E 1 + φ 2 σ 2 : E 2 + I • [[u ]] (41) = Σ 1 : E 1 + Σ 2 : E 2 + I • V where ε el = ε -E + E , E = χ 1 E 1 + χ 2 E 2 and [[ ]] = 2 - 1 .
This extended Hill-Mandel lemma shows that the average work of deformation of σ in the local "elastic" strain ε el is equal to the macroscopic work of deformation of the multiphase model characterized by the partial stresses Σ i = φ i σ i , dual to the generalized strains E i , and the interaction force I, dual to the average relative displacement

V = u 2 -u 1 .
Using the set (E, ∆E) of strain parameters, the previous expression can also be rewritten as:

σ : ε el = σ : E + ( σ 2 -σ 1 ) : ∆E + I • [[u ]] (42) = Σ : E + ∆Σ : ∆E + I • V
in which Σ can be interpreted as the classical macroscopic stress, associated with the macroscopic strain E, and the average stress mismatch ∆Σ = σ 2 -σ 1 = [[σ ]] between phase 2 and phase 1, associated with ∆E. Again, setting ∆E and I to zero, the classical homogenization framework is recovered.

The extended Hill-Mandel lemma (41) also provides bounding theorems since σ : ε ≤ ψ(ε ) + ψ * (σ ) for any (σ , ε ) with ψ and ψ * being the local strain energy and stress energy densities. Hence, ∀(σ , u ) ∈ SA(I) × KA:

σ : ε el = Σ 1 : E 1 + Σ 2 : E 2 + I • V ≤ ψ(E + ∇ s u ) + ψ * (σ ) (43)
in which the equality is attained for the solution pair (u, σ) of problem (37).

Resolution of the auxiliary problem

In this subsection, we consider the auxiliary problem (37) controlled by (E 1 , E 2 , I). The solution can be defined up to a rigid body motion on the unit cell. To avoid this indeterminacy, we will look for displacement solutions of zero average on the unit cell. Owing to its linear character, the solution is a linear combination of the different loading parameters:

u(y) = a 1 (y) : E 1 + a 2 (y) : E 2 + d(y) • I ( 44
)
where a 1 , a 2 , d are displacement concentration tensors relative to the periodic fluctuation u, respectively corresponding to the solution of elementary auxiliary problems in which (E 1 ∈ {e i ⊗ e j for i, j ∈ {1, 2, 3}}, E 2 = 0, I = 0), (E 1 = 0, E 2 ∈ {e i ⊗ e j for i, j ∈ {1, 2, 3}}, I = 0) and (E 1 = 0, E 2 = 0, I ∈ {e i for i ∈ {1, 2, 3}}). The local elastic ε el (y) and total ε(y) strains are then also expressed as linear combinations of the loading parameters:

E (y) + ∇ s u(y) = ε el (y) = A 1 (y) : E 1 + A 2 (y) : E 2 + D (y) • I (45) E + ∇ s u(y) = ε(y) = A 1 (y) : E 1 + A 2 (y) : E 2 + D (y) • I (46)
where:

A i (y) = χ i (y)I + ∇a i (y) i = 1, 2 (47)
A i (y) = φ i I + ∇a i (y) i = 1, 2 (48)
D (y) = ∇d(y) (49)
I is the 4th-order identity tensor and ∇ is such that ( ∇X ) ijkl = (X ikl,j + X jkl,i )/2 for any 3rd-order tensor X and ( ∇Y ) ijk = (Y ik,j + Y jk,i )/2 for any 2nd-order tensor Y . Note that tensors A i (y), A i (y) are closely related to the eigenstrain influence functions introduced in [START_REF] Dvorak | On Transformation Strains and Uniform Fields in Multiphase Elastic Media[END_REF]. Let us also remark that if E 1 = E 2 = E and I = 0, we recover the standard homogenization framework, hence, A(y) = A 1 (y) + A 2 (y) = A 1 (y) + A 2 (y) corresponds to the standard strain concentration tensor.

From [START_REF] Nguyen | Analyse sismique des ouvrages renforcés par inclusions rigides à l'aide d'une modélisation multiphasique[END_REF], the local stress is then given by:

σ(y) = C i : A 1 (y) : E 1 + C i : A 2 (y) : E 2 + C i : D(y) • I ∀y ∈ A i (50)
From the definition of the partial stresses Σ i and of the relative displacement V , we obtain the macroscopic constitutive equations:

Σ 1 = φ 1 C 1 : A 11 : E 1 + φ 1 C 1 : A 12 : E 2 + φ 1 C 1 : D 1 • I (51a) Σ 2 = φ 2 C 2 : A 21 : E 1 + φ 2 C 2 : A 22 : E 2 + φ 2 C 2 : D 2 • I (51b) V = a 21 -a 11 : E 1 + a 22 -a 12 : E 2 + d 2 -d 1 • I (51c) with A ij = A j (y) i and a ij = a j (y) i .
Summing up the first two relations yields the macroscopic total stress:

Σ = Σ 1 + Σ 2 = C(y) : A 1 (y) : E 1 + C(y) : A 2 (y) : E 2 + C(y) : D (y) • I (52)
and, again, if E 1 = E 2 = E and I = 0, it reduces to the standard macroscopic behavior:

Σ = C(y) : (A 1 (y) + A 2 (y)) : E = C(y) : A(y) : E = C hom : E (53)
where κ = [[d]] -1 is the interaction stiffness. Obviously, the previous first two relations can be inverted in order to express the generalized strains E 1 , E 2 as functions of both partial stresses Σ 1 , Σ 2 as follows:

E 1 = F 11 : Σ 1 + F 12 : Σ 2 (58a) E 2 = (F 12 ) T : Σ 1 + F 22 : Σ 2 (58b) I = κ • V (58c)
with the following relations between the partial stiffness tensors D ij and the partial compliance tensors F ij :

F 11 = D 11 -D 12 : (D 22 ) -1 : (D 12 ) T -1 (59)
F 22 = D 22 -(D 12 ) T : (D 11 ) -1 : D 12 -1 (60)
F 12 = -F 11 : D 12 : (D 22 ) -1 = -(D 11 ) -1 : (D 12 ) T : F 22 (61)

Formal justification through asymptotic expansion and discussion

Let us now interpret the former auxiliary problem and the link between macroscopic and microscopic quantities by means of a formal two-scale asymptotic expansion [START_REF] Sánchez-Palencia | Non-homogeneous media and vibration theory[END_REF]. The x space variable will denote the macroscopic scale (slow variable) and y = x/ the microscopic space (fast) variable. Mechanical fields will be represented as a series in terms of the scale factor . In the following, we will use standard results of spatial derivatives of total and partial averages in a periodic setting [START_REF] Dormieux | Microporomechanics[END_REF].

Let us consider the following first-order expansion for the stress field:

σ(x, y) = σ (0) (x, y) + σ (1) (x, y) = C i : A(y) : E + (C i : ∆A(y) : ∆E + C i : D (y) • I) (62)
with ∆A(y) = A 2 (y)/φ 2 -A 1 (y)/φ 1 and where the supscripts in parenthesis denote the order in the expansion. This expansion can be directly derived from [START_REF] Soubestre | Non-local dynamic behavior of linear fiber reinforced materials[END_REF] in which it has been assumed that the partial generalized strains are equal to the macroscopic strain at the lowest-order plus a first-order correction depending on ∆E i.e.

E i = E + (-1) i ∆E/φ i .
Similarly, the interaction force real amplitude has been assumed of order 1 in . The three macroscopic variables E, ∆E, I are assumed to be functions of the macroscopic variable x.

The previous expansion is such that:

• The total average is equal to the macroscopic stress:

σ = σ (0) + σ (1) = σ (0) = C hom : E = Σ (63)
where the fact that σ (1) = 0 can be obtained from ( 47)-( 49) and periodicity.

• The partial stresses are given by:

Σ i = φ i σ i = φ i σ (0) i + φ i σ (1) i = (D i1 + D i2 ) : E + D i2 φ 2 - D i1 φ 1 : ∆E (64)
that is they are given at the lowest order by the localization obtained from classical homogenization in addition to a first-order correction depending on the relative strain ∆E.

• Expressing the local balance equation at each order (assuming body forces of order 0) gives:

-1 (div y σ (0) ) + 0 (div x σ (0) + div y σ (1) + ρF (0) ) + 1 (div x σ (1) ) = 0 (65)

Since σ (0) is the solution to ( 38) with E 1 = E 2 = E and I = 0, the term of order -1 is zero by construction. From the periodicity of σ (0) , σ (1) and ( 63), the total average of (65) yields the macroscopic equilibrium equation for the total stress:

div x σ (0) + ρF (0) + div x σ (1) = div x σ (0) + ρF (0) + div x σ (1) = div x Σ + ρF (0) = 0 (66)
As a result, if the total stress satisfies the macroscopic balance equation, its microscopic counterpart is satisfied on average on the unit cell.

• Taking now the average on each phase, one has:

φ i div x σ (0) i + φ i div y σ (1) i + φ i ρF (0) i + φ i div x σ (1) i = 0 div x (φ i σ (0) i + φ i σ (1) i ) -(-1) i I + φ i ρF (0) i = 0 div x Σ i -(-1) i I + φ i ρF (0) i = 0 ( 67
)
where we used that σ (1) is the solution to [START_REF] Hegemier | A continuum mixture theory of wave propagation in laminated and fiber reinforced composites[END_REF] with non-zero interaction force I i.e. div y σ (1) = -ζ(y)I. As a result, if the partial stresses satisfy the macroscopic balance equations for each phase ( 18)- [START_REF] Boutin | Generalized inner bending continua for linear fiber reinforced materials[END_REF] in which the partial body forces are now given a clear microscopic interpretation from (67), the microscopic balance equation is satisfied on average on each phase.

Similarly, let us now consider the following expansion for the displacement field:

U (x, y) = U (0) (x) + U (1) (x) + u (1) (x, y) + 2 U (2) (x) + u (2) (x, y) (68)
u (1) (x, y) = a (y) : E(x) (69)
u (2) (x, y) = ∆a (y) : ∆E(x) + d(y) • I(x) (70)
with a (y) = a 1 (y) + a 2 (y) and ∆a (y) = a 2 (y)/φ 2 -a 1 (y)/φ 1 . In the proposed expansion u (1) (x, y) + u (2) (x, y) corresponds exactly to the displacement solution to the auxiliary problem (38) using the same hypothesis regarding the order of magnitudes of the macroscopic variables as for the previous stress expansion, the total displacement being the addition of this periodic fluctuation to a macroscopic displacement expanded up to the second-order. The proposed expansion is therefore such that:

• The average displacement is given by:

U (x) = U (0) (x) + U (1) (x) + 2 U (2) (x) (71)
due to the fact that the periodic fluctuation is of zero-average.

• The average displacements on each phase are given by:

U i (x) = U (x) + a i : E(x) + 2 ∆a i : ∆E(x) + d i • I(x) (72)
their relative difference being given by:

U 2 -U 1 = [[a ]] : E(x) + 2 [[∆a ]] : ∆E(x) + [[d]] • I(x) = V (x) (73)
showing that the average displacement on each phase is equal to the macroscopic displacement at the lowest-order plus a first-order correction depending on V

• The local strain is given by:

ε(x, y) = ∇ s x U + -1 ∇ s y U = 0 ∇ s x U (0) + ∇ s y u (1) + 1 ∇ s x U (1) + ∇ s x u (1) + ∇ s y u (2) + 2 ∇ s x U (2) + ∇ s x u (2) (74) = ∇ s x U + ∇ s y (a : E) + ∇ s y (∆a : ∆E) + D • I + ∇ s x u (1) + 2 ∇ s x u (2)
Its total average is then ε = ∇ s x U = ε because of periodicity and zero-average of u (1) , u (2) . Comparing now this strain with that obtained from the stress expansion (62) yields:

ε -(C i ) -1 : σ = ∇ s x U -E + ∇ s x u (1) -ζ(y)∆E + 2 ∇ s x u (2) (75)
On the one hand, the zero-order term is null if the macroscopic strain is identified with the macroscopic gradient of the average displacement E = ∇ s x U , as expected from the standard homogenization approach. On the other hand, it can be seen that nor the first-order neither the second-order terms can be null in general. However, both terms are of zero-average over the unit cell so that both strain fields are equal on average on the unit cell. Besides, taking the difference of both phase-averages, one obtains:

[[ε -(C i ) -1 : σ]] = ∇ s
x [[u (1) ]] + [[u (2) )]] -∆E/(φ

1 φ 2 ) = (∇ s x V -∆E/(φ 1 φ 2 )) (76)
As a result if ∆E(x) = φ 1 φ 2 ∇ s x V (x), the previous right-hand side vanishes. As a result, since

E 2 -E 1 = ∆E/(φ 1 φ 2 ) = ∇ s x V = ∇ s x U 2 -∇ s x U 1 and since E = φ 1 E 1 + φ 2 E 2 = ∇ s x U = φ 1 ∇ s x U 1 + φ 2 ∇ s
x U , we finally obtain that:

E 1 (x) = ∇ s x U 1 (x), E 2 (x) = ∇ s x U 2 (x) (77)
As a conclusion, it can be shown that the solutions to the proposed auxiliary problem can be associated with microscopic stress and displacement fields (62)-(68) such that local equilibrium is verified on total average and on phase-average and such that the stress-strain constitutive equation is satisfied locally at order 0 and in total and phase-average provided that:

• the macroscopic displacement of both phases are identified with the average displacement on the corresponding phase:

U i = U i
• the generalized strains E 1 , E 2 are given by the compatibility conditions (77)

• the generalized stresses Σ 1 , Σ 2 , I satisfy the equilibrium equations ( 18)-( 19)

• the generalized stresses and strains are linked by the macroscopic constitutive equations [START_REF] Tran | Modélisation à plusieurs échelles d'un milieu continu hétérogène aléatoire[END_REF] This derivation shows that the multiphase model can be thought of as a first-order correction to the standard Cauchy continuum obtained at order 0, in which the relative displacement, the strain mismatch and the interaction force control the first-order term. One key feature is that the higher-order terms of the balance equation and strain compatibility equations are not verified exactly, contrary to other higher-order approaches such as [START_REF] Boutin | Microstructural effects in elastic composites[END_REF][START_REF] Smyshlyaev | On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media[END_REF] yielding strain gradient models. Instead, these equations are satisfied on average on each phase, yielding a simpler macroscopic model. Although the previous derivation is purely formal, it is by no means a proof of convergence, it can be anticipated that the model will be able to better describe what happens near the boundary, since it is driven by average fields in each phase and not by average fields on the whole unit cell only. The validity of the proposed model will then be assessed by means of the examples of section 5.

Remark: In problem [START_REF] Dormieux | Microporomechanics[END_REF], the fictitious body force enabling to control the interaction force I has been chosen as b(y) = ζ(y)I with ζ being phase-wise uniform. However, any function ζ(y) such that ζ = 0 and φ i ζ i = (-1) i would also fulfill [START_REF] Dumontet | Study of a boundary layer problem in elastic composite materials[END_REF]. The whole procedure would therefore be very similar, the main change being that the dual quantity associated with I would now be given by V = ζU , reducing to the difference of phase-average with the retained choice. As a result, the macroscopic displacements U i and associated strains E i would then have the following micromechanical interpretation:

U 1 = U -φ 2 ζU , U 2 = U + φ 1 ζU (78)
E 1 = ∇ s x U -φ 2 ∇ s x ζU , E 2 = ∇ s x U + φ 1 ∇ s x ζU (79)
A whole range of multiphase models can then be obtained from different choices of ζ, yielding most certainly different values for the associated interaction stiffness. Exploring the differences between such choices is out-of-the scope of the present paper, our initial choice being the most simple and having the most satisfactory microscopic interpretation of the macroscopic displacements and strains.

Derivation of effective properties

The homogenization procedure presented in section (3) can be easily extended to the situation of a multiphase model consisting of n > 2 phases. In general, the determination of the various effective moduli needs to resort to a numerical resolution of the auxiliary problem (37) (or [START_REF] Hegemier | A continuum mixture theory of wave propagation in laminated and fiber reinforced composites[END_REF]), by means, for instance, of a finite-element discretization. It is to be noted that this poses no specific difficulty since taking into account periodic boundary conditions, eigenstrains and body forces can be straightforwardly done in finite-element codes. The main difference compared to the standard homogenization framework comes from a larger number of elementary solicitations to be solved which induces a negligible additional cost in a linear elastic setting.

However, many situations can benefit from a direct evaluation or estimation of the effective moduli as a function of the constituents properties and the microstructure morphology. We give here different results extending classical continuum micromechanics approaches to the multiphase continuum framework, enabling to derive either an exact evaluation, an estimate or bounds on the effective moduli. As regards the results obtained in this section, a few of them have already been established by other authors in some specific situations such as homogenization of fiber-reinforced composites or porous materials. Our proposed homogenization procedure is therefore shown to extend or encompass a wide range of results arising in different fields.

Partial stiffness tensors of a two-phase material

In the case of a two-phase (n = 2) material, no estimation of the concentration tensors A ij is needed if the standard overall homogenized stiffness C hom is known in advance (e.g. from a numerical resolution, experimental results or also from an estimation based on standard micromechanical schemes) [START_REF] Dvorak | On uniform fields in heterogeneous media[END_REF][START_REF] Zaoui | Continuum micromechanics: survey[END_REF].

Indeed, using the fact that ∇ s u = 0, it can be deduced from (47) that:

A 1 = φ 1 A 11 + φ 2 A 21 = φ 1 I, A 2 = φ 1 A 12 + φ 2 A 22 = φ 2 I (80)
Using now relations (54) (in particular the symmetry relations of D ij ) and (56), it can be shown that:

D 11 = φ 1 C 1 -C 1 : [[C]] -1 : ∆C : [[C]] -1 : C 1 (81a) D 22 = φ 2 C 2 -C 2 : [[C]] -1 : ∆C : [[C]] -1 : C 2 (81b)
D 12 = C 1 : [[C]] -1 : ∆C : [[C]] -1 : C 2 (81c) with [[C]] = C 2 -C 1 and ∆C = C -C hom .
We can also easily check that

D 11 + D 22 + D 12 + (D 21 ) T = C hom .
4.2. Partial stiffness tensors for a two-phase material with stiff inclusions in small volume fraction The phenomenological multiphase continuum model used in [START_REF] De Buhan | Elastic, Plastic and Yield Design of Reinforced Structures[END_REF][START_REF] Sudret | Modélisation multiphasique des ouvrages renforcés par inclusions[END_REF][START_REF] Sudret | Modélisation multiphasique de matériaux renforcés par inclusions linéaires[END_REF] has been developed for the case of a matrix (phase 1) reinforced by a periodic array of linear inclusions (phase 2) which are assumed to be much stiffer than the matrix (C 2 C 1 ). At the same time, the volume fraction of such inclusions is also supposed to be very small φ 2 1. More precisely, it is assumed that φ 2 C 2 admits a finite limit C 0 when φ 2 → 0 so that inclusions still have an effect with vanishing volume fraction. Considering this particular situation, the previous effective moduli (81a)-(81c) admit a simpler form. We have indeed that:

[[C]] ∼ C 2 (82) C → C 1 + C 0 (83)
Due to the last relation and since C hom C , the homogenized stiffness C hom will also tend to a finite value denoted by C hom,0 and which depends on the microstructure. Using all these relations, we have:

D 11 → C 1 (84)
D 22 → C hom,0 -C 1 (85)
D 12 → 0 (86)
In this limit, the stress-strain behavior of one phase is therefore decoupled from the other and the matrix phase partial behavior is exactly that of its constitutive material.

In the particular case of (transverse) isotropic linear inclusions oriented along direction e x , it can be shown, using for instance the results of [START_REF] Hashin | The elastic moduli of fiber-reinforced materials[END_REF], that:

lim φ 2 →0 C hom = C hom,0 = C 1 + E 0 e x ⊗ e x ⊗ e x ⊗ e x (87)
where E 0 = lim

φ 2 →0 φ 2 E 2 (88)
with E 2 being the (axial) Young's modulus of the inclusion material. In this specific case:

D 22 → E 0 e x ⊗ e x ⊗ e x ⊗ e x (89)
thereby justifying a posteriori the phenomenological multiphase constitutive relations ( 15)-( 16) in the case of stiff reinforcements.

Bounds on the effective properties

From [START_REF] Nairn | On the use of shear-lag methods for analysis of stress transfer in unidirectional composites[END_REF], one can easily deduce the following bounds on the multiphase model effective properties if the controlled loading parameters are (E 1 , E 2 , I):

φ 1 σ 1 : E 1 + φ 2 σ 2 : E 2 -ψ * (σ ) ≤ Ψ(E 1 , E 2 , I) ≤ ψ(E + ∇ s u ) -I • [[u ]] (90)
for any σ ∈ SA(I) and any u ∈ KA.

The generalization of Voigt's upper bound to the multiphase effective moduli is obtained by considering a zero periodic fluctuation u = 0 in the upper bound estimate (90), that is:

Ψ(E 1 , E 2 , I) ≤ ψ(E ) = 1 2 E 1 : (φ 1 C 1 ) : E 1 + 1 2 E 2 : (φ 2 C 2 ) : E 2 (91)
This upper bound thus corresponds to the situation in which both phases are perfectly bonded at the macroscopic level (I being indeterminate but with V = 0 ⇔ κ → ∞) and in which the partial stiffness tensor are given by4 D ii = φ i C i and D ij = 0 for i = j. In this situation, the multiphase model reduces to a standard one-phase Cauchy medium with effective properties C .

A generalized Reuss lower bound can also be derived by considering a uniform stress field σ (y) = Σ, which is statically admissible with I = 0, in the lower bound estimate (90), that is:

Ψ(E 1 , E 2 , 0) ≥ Σ : (φ 1 E 1 + φ 2 E 2 ) - 1 2 Σ : C -1 : Σ (92)
This lower bound holding true for any choice of Σ, the right-hand side can be replaced by its maximum over all Σ, yielding:

Ψ(E 1 , E 2 , 0) ≥ max Σ {Σ : (φ 1 E 1 + φ 2 E 2 ) - 1 2 Σ : C -1 : Σ} (93) ≥ 1 2 (φ 1 E 1 + φ 2 E 2 ) : C -1 -1 : (φ 1 E 1 + φ 2 E 2 )
This leads to an effective behavior in which both phases are completely debonded (I = 0 ⇔ κ = 0) and in which the partial stiffness tensors are given by

D ij = φ i φ j C -1 -1 .
As a result, simple bounds considering uniform strain or stress states do not produce interesting results regarding the interaction stiffness since its behavior degenerates to either fully bonded or fully debonded at the macroscopic level. It is interesting to point out that this is reminiscent of the difficulty of deriving simple bounds when estimating the effective permeability of a porous medium [START_REF] Bignonnet | Caractérisation expérimentale et modélisation micro-mécanique de la perméabilité et la résistance de roches argileuses[END_REF][START_REF] Bignonnet | Upper bounds on the permeability of random porous media[END_REF][START_REF] Boutin | Study of permeability by periodic and self-consistent homogenisation[END_REF]. Obviously, obtaining more accurate bounds, using for instance the Hashin-Shtrikman framework, would be of great value but is out of the scope of the present paper.

Effective properties of a multi-layered medium

We consider here the case of a multi-layered medium consisting of alternating isotropic (Lamé coefficients λ i , µ i ) layers of thicknesses t 1 and t 2 in the y-direction (see Figure 3). We will restrict here to the plane-strain case in the (x, y)-plane. The unit cell will be characterized by a square domain of dimensions s × s (actually, the length in the x-direction will play no role in the subsequent calculations) with s = t 1 + t 2 and with η = t 2 /s being the volume fraction of phase 2 (later identified with the reinforcement phase). We also introduce the non-dimensional coordinates (y 1 , y 2 ) = (x/s, y/s).

Derivation of partial stiffness moduli

As previously mentioned, the partial stiffness moduli D ij can be expressed directly from the knowledge of the homogenized stiffness tensor C hom in the case of a biphasic material. The expression of C hom is available, for instance, in [START_REF] De Buhan | Elastic, Plastic and Yield Design of Reinforced Structures[END_REF]. Instead, we will solve directly the corresponding auxiliary problems since their resolution can be easily generalized to more than two phases for which no simple expression in terms of C hom is available.

Let us consider problem (P I ) in which phase 1 is subjected to generalized strains

E 1 = E xx E xy E xy E yy (x,y)
, E 2 = 0 and with no body forces (I = 0). The displacement solution is searched in the form of U (y) = φ 1 E 1 • y + u(y 2 )e x + v(y 2 )e y where u(y 2 ), v(y 2 ) are piecewise-linear periodic fluctuations in both directions such that u (y 2 ) (resp. v (y 2 )) is equal to a given constant A i (resp. B i ) in phase i. The periodicity condition requires that A = B = 0 so that:

u(y 2 ) =    A 2 y 2 for 0 ≤ y 2 ≤ η (in phase 2) A 2 η 1 -η (1 -y 2 ) for η ≤ y 2 ≤ 1 (in phase 1) (94) v(y 2 ) =    B 2 y 2 for 0 ≤ y 2 ≤ η (in phase 2) B 2 η 1 -η (1 -y 2 ) for η ≤ y 2 ≤ 1 (in phase 1) (95)
Such a choice yields piecewise uniform stress states in each phases:

σ xx =    M 1 E xx + λ 1 E yy -B 2 η 1 -η in phase 1 λ 2 B 2 in phase 2 (96)
σ yy =    M 1 E yy -B 2 η 1 -η + λ 1 E xx in phase 1 M 2 B 2 in phase 2
(97)

σ xy =    µ 1 2E xy -A 2 η 1 -η in phase 1 µ 2 A 2 in phase 2 (98)
where M i = λ i + 2µ i is the oedometric modulus of phase i. Expressing traction continuity at the interface y 2 = η fixes the values of A 2 and B 2 and thus fully determines the solution. Partial stresses Σ 1 and Σ 2 in both phases can then be completely expressed as a function of the macroscopic generalized strain E 1 components as follows:

   Σ 1 xx Σ 1 yy Σ 1 xy    = (1 -η)        M 1 -η λ 2 1 M * λ 1 -η λ 1 M 1 M * 0 λ 1 -η λ 1 M 1 M * (1 -η) M 1 M 2 M * 0 0 0 (1 -η) µ 1 µ 2 µ *           E xx E yy 2E xy    (99)    Σ 2 xx Σ 2 yy Σ 2 xy    = η(1 -η)        λ 1 λ 2 M * λ 2 M 1 M * 0 λ 1 M 2 M * M 1 M 2 M * 0 0 0 µ 1 µ 2 µ *           E xx E yy 2E xy    ( 100
)
where X * = φ 2 X 1 + φ 1 X 2 for any X. These two equations respectively represent the plane-strain components of tensors D 11 and D 21 = (D 12 ) T . A similar derivation considering problem (P II ) with ,y) and I = 0 yields the corresponding expressions for D 22 :

E 1 = 0, E 2 = E xx E xy E xy E yy (x
   Σ 2 xx Σ 2 yy Σ 2 xy    = η        M 2 -(1 -η) λ 2 2 M * λ 2 -(1 -η) λ 2 M 2 M * 0 λ 2 -(1 -η) λ 2 M 2 M * η M 1 M 2 M * 0 0 0 η µ 1 µ 2 µ *           E xx E yy 2E xy    (101)
Note that the obtained results also agree with [START_REF] Mcniven | A mixture theory for elastic laminated composites[END_REF].

Let us finally remark that adopting a similar reasoning for a generalized loading of the form E 1 = E 0 0 1 (x,y) and E 2 = E 0 0 2 (x,y) where 1 and 2 are left free such that σ yy (y) = 0, one finds:

F 11 xxxx = 1/(φ 1 E * 1 ) (102)
F 22 xxxx = 1/(φ 2 E * 2 ) ( 103
)
F 12 xxxx = 0 ( 104
)
where

E * i = E i 1 -ν i
1 + ν i with E i being the Young's modulus of phase i and ν i its Poisson ratio.

Derivation of the interaction stiffness

Because of the unit-cell central symmetry, there is no coupling between the partial stresses and the interaction force in the multiphase constitutive equations (57). The only remaining parameters are therefore those related to the interaction stiffness κ which is of the form κ = κ x e x ⊗ e x + κ y e y ⊗ e y .

Problem (P III

) where E 1 = E 2 = 0 and I = I x e x + I y e y is now considered (see Figure 3). As before the displacement field is sought of the form U (y) = u(y 2 )e x + v(y 2 )e y with u and v being periodic. Expressing local equilibrium, displacement and traction continuities, one can easily show that the solution of (P III ) (up to a rigid-body motion) is given by:

u(y 2 ) =        s 2 I x 2(1 -η)µ 1 (y 2 2 -(1 + η)y 2 + η) in phase 1 s 2 I x 2ηµ 2 y 2 (η -y 2 ) in phase 2 (105) v(y 2 ) =        s 2 I y 2(1 -η)M 1 (y 2 2 -(1 + η)y 2 + η) in phase 1 s 2 I y 2ηM 2 y 2 (η -y 2 ) in phase 2 (106)
Computation of the average relative displacement between both phases yields:

V = 1 η η 0 U (y 2 )dy 2 - 1 1 -η 1 η U (y 2 )dy 2 = s 2 µ * 12µ 1 µ 2 I x e x + s 2 M * 12M 1 M 2 I y e y ( 107
)
from which we deduce that:

κ x = 12 s 2 1/µ , κ y = 12 s 2 1/M (108)
Let us remark that these expressions are consistent with the derivations made in [START_REF] Bedford | Toward A Diffusing Continuum Theory of Composite Materials[END_REF][START_REF] Mcniven | A mixture theory for elastic laminated composites[END_REF][START_REF] Stern | Wave propagation in elastic laminates using a multi-continuum theory[END_REF].

Effective properties of unidirectional fiber-reinforced materials

In this section, the derivation of the multiphase model constitutive parameters is discussed in the case of unidirectional fiber-reinforced composites aligned along the z-direction (see Figure 4a). Again, it is possible to numerically solve the auxiliary problem using a finite-element discretization of the real fiber-matrix microstructure. However, many analytical results are available when adopting, for instance, a simplified representation of the composite microstructure in the form of a composite cylinder assemblage (CCA) model in which a fiber cylinder of radius a is embedded in a matrix cylinder of radius b (Figure 4b). In [START_REF] Hegemier | A continuum mixture theory of wave propagation in laminated and fiber reinforced composites[END_REF], a procedure based on stress and displacement expansions with r is proposed but does not lead to the correct expressions of the macroscopic (and hence the partial) stiffness and the interaction stiffness.

Derivation of partial stiffness moduli

As regards the partial stiffness moduli D ij , it has already been mentioned that they can be derived from the sole knowledge of the standard homogenized moduli C hom of the corresponding composite. Classical results from Hashin and Rosen [START_REF] Hashin | Analysis of Properties of Fiber Composites With Anisotropic Constituents[END_REF][START_REF] Hashin | The elastic moduli of fiber-reinforced materials[END_REF] provide either exact expressions or bounds on the effective moduli in the case of isotropic or transverse isotropic constituents.

Some particular loading cases of the auxiliary problem offer simple exact solutions based on the CCA model. For instance, consider the problem for which E 2 = Ee z ⊗ e z and E 1 = 0 in the case of isotropic constituents. Building a piecewise-linear periodic fluctuation displacement field such that ∇ s u = A i (e x ⊗ e x + e y ⊗ e y ) in phase i with φ 1 A 1 + φ 2 A 2 = 0 to satisfy the periodicity requirements leads to the following uniform stresses in each phases:

σ zz = 2λ 1 A 1 in phase 1 (λ 2 + 2µ 2 )E + 2λ 2 A 2 in phase 2
(109)

σ xx = σ yy = 2(λ 1 + µ 1 )A 1 in phase 1 λ 2 E + 2(λ 2 + µ 2 )A 2 in phase 2 (110)
Continuity of the transverse stresses further imposes that 2(λ 1 +µ 1 )A 1 = λ 2 E +2(λ 2 +µ 2 )A 2 , so that:

A 1 = φ 2 λ 2 2 λ + µ * E, A 2 = -φ 1 λ 2 2 λ + µ * E (111)
We then deduce that:

D 22 zzzz = φ 2 λ 2 + 2µ 2 - φ 1 λ 2 λ + µ * (112) D 12 zzzz = φ 1 φ 2 λ 1 λ 2 λ + µ * (113) D 12 zzxx = D 12 zzyy = D 22 xxzz = D 22 yyzz = φ 1 φ 2 λ 2 (λ 1 + µ 1 ) λ + µ * (114)
Exchanging the role of 1 and 2 we also easily deduce the corresponding expressions for D 11 zzzz , D 11 xxzz and D 12 xxzz . Extension of these expressions for transverse isotropic constituents is also straightforward.

Similarly, let us now consider the case where

E 2 =   2 0 0 0 2 0 0 0 E   (x,y,z) and E 1 =   1 0 0 0 1 0 0 0 0   (x,y,z)
where 1 and 2 are left free such that Σ xx = Σ yy = 0. Considering the same periodic fluctuation field, one now has:

σ zz = 2λ 1 (A 1 + 1 ) in phase 1 (λ 2 + 2µ 2 )E + 2λ 2 (A 2 + 2 ) in phase 2 ( 115
)
σ xx = σ yy = 2(λ 1 + µ 1 )(A 1 + 1 ) in phase 1 λ 2 E + 2(λ 2 + µ 2 )(A 2 + 2 ) in phase 2 (116)
Again, because of the continuity condition on the transverse stresses, 2(λ

1 + µ 1 )(A 1 + 1 ) = λ 2 E + 2(λ 2 + µ 2 )(A 2 + 2
) which should both be zero to satisfy the macroscopic condition

Σ xx = Σ yy = 0. One finds A 1 = -1 and φ 1 1 + φ 2 2 = E xx = E yy = - λ 2 2(λ 2 + µ 2 ) φ 2 E. This results in Σ 1 = 0 and Σ 2 = φ 2 λ 2 + 2µ 2 - λ 2 2 λ 2 + µ 2 Ee z ⊗ e z = φ 2 E 2 Ee z ⊗ e z where E 2
is the Young's modulus of phase 2. Obviously, a similar result is obtained when exchanging the role of phase 1 and 2. One can then deduce that:

F 11 zzzz = 1/(φ 1 E 1 ) (117)
F 22 zzzz = 1/(φ 2 E 2 ) ( 118
)
F 12 zzzz = 0 (119)
Note that these results are also easily extended to the case of transverse isotropy for each phases.

Interaction stiffness

Assuming a transverse isotropic behavior of axis e for both constituents, the interaction stiffness is necessary of the following form:

κ = κ L e z ⊗ e z + κ T (e x ⊗ e x + e y ⊗ e y ) ( 120
)
Longitudinal interaction stiffness κ L . The longitudinal interaction stiffness is obtained by solving an auxiliary problem with self-equilibrated uniform body forces along the cylinder axis:

div σ + I/φ 2 e z = 0 for 0 < r < a (121) div σ -I/φ 1 e z = 0 for a < r < b (122)
in addition to the traction and displacement continuity at the interface r = a and the constitutive equations and where φ 2 = a 2 /b 2 = 1 -φ 1 . This problem is solved by considering purely longitudinal displacements u = u i (r)e z in both phases associated with shear stresses

σ i rz = µ i du i dr
where µ i denotes the longitudinal shear modulus of phase i. The displacements are then solutions to:

d 2 u 2 dr 2 + 1 r du 2 dr = - I µ 2 φ 2 (123) d 2 u 1 dr 2 + 1 r du 1 dr = I µ 1 φ 1 (124)
which can be solved as:

u 2 (r) = - I 4µ 2 φ 2 r 2 + A (125) u 1 (r) = I 4µ 1 φ 1 (r 2 -b 2 ) + B ln(r/b) + C (126)
The constants A, B and C are then obtained by solving the system of equations corresponding to the following conditions:

u 1 (a) = u 2 (a) (127) σ 1 rz (a) = σ 2 rz (a) (128) u 1 (b) = 0 (129)
Note that the last condition enables to remove the rigid body motion, it can be equivalently replaced by another condition fixing this indeterminacy. One finds C = 0 and:

A = b 2 µ 1 φ 1 -µ 2 (φ 1 + ln(φ 2 )) 4µ 1 µ 2 φ 1 (130) B = - b 2 I 2µ 1 φ 1 (131)
Finally, the average displacement jump between both phases is computed as:

V = 2 a 2 u 2 (r)rdr - 2 b 2 -a 2 b a u 1 (r)rdr (132)
The corresponding longitudinal interaction stiffness κ L = I/V is given, after computations, by:

κ L = 8µ 1 µ 2 φ 1 b 2 (µ 1 φ 1 -µ 2 (2 + φ 1 + 2(ln φ 2 )/φ 1 )) (133)
In the limit of rigid inclusions µ 2 µ 1 the previous expression reduces to:

κ L → κ rigid L = - 8µ 1 φ 1 b 2 (2 + φ 1 + 2(ln φ 2 )/φ 1 ) (134)
It is worth noting that this expression gives a larger value (up to 50% for large reinforcement volume fractions) of the interaction stiffness compared to Sudret's estimate [START_REF] Sudret | Modélisation multiphasique des ouvrages renforcés par inclusions[END_REF]:

κ Sudret L = - 4µ 1 b 2 (1 + (ln φ 2 )/φ 1 ) (135)
derived from the solution of a rigid pull-out displacement of the fiber. Interestingly, expression (134) also coincides for small fiber volume fractions (φ 2 → 0) with the permeability estimate derived by Boutin [START_REF] Boutin | Study of permeability by periodic and self-consistent homogenisation[END_REF] using a static approach for a flow in a fluid-solid cylindrical inclusion:

κ Boutin L = - 8µ 1 b 2 (φ 1 (3 -φ 2 ) + 2 ln φ 2 ) ( 136
)
Transverse interaction stiffness κ T . The transverse interaction stiffness κ T is now obtained from a similar auxiliary problem in which uniform distributed forces are oriented along e x : div σ + I/φ 2 e x = 0 for 0 < r < a (137) div σ -I/φ 1 e x = 0 for a < r < b

We look for displacements of the following form u = U i (r) cos θe r + W i (r) sin θe θ in each phase. As a consequence, the plane components of the stress tensor are of the form σ rr = Σ rr (r) cos θ, σ θθ = Σ θθ (r) cos θ and σ rθ = Σ rθ (r) sin θ in both phases.

It can be shown that the following displacements will induce stresses in equilibrium with the prescribed body forces:

U 1 (r) = A + Br 2 + Cr -2 + A + D ln r (139) W 1 (r) = -A - I φ 1 (2µ 1 -k 1 ) r 2 + γ 1 Br 2 + Cr -2 - k 1 k 1 + 2µ 1 D -D ln(r) (140)
U 2 (r) = A + B r 2 (141)
W 2 (r) = -A + I φ 2 (2µ 2 -k 2 ) r 2 + γ 2 B r 2 (142)
with γ i = (3k i + 2µ i )/(2µ i -k i ), µ i denoting now the transverse shear modulus and k i being the transverse bulk modulus used in the replacement scheme of [START_REF] Hashin | Analysis of Properties of Fiber Composites With Anisotropic Constituents[END_REF].

The six unknown constants are then determined by solving the linear system associated with the displacement and traction continuities at the fiber/matrix interface r = a (4 conditions), one condition fixing the rigid body motion (e.g. U 2 (0) = 0) and the antiperiodicity conditions for the traction at the outer boundary. The latter imposes that Σ rr (b) = Σ rθ (b) = 0 which seems to lead to an overdetermined system. However, because the body forces are self-equilibrated, the divergence theorem shows that both conditions are in fact linearly dependent so that only one is necessary to solve the problem.

The transverse interaction stiffness is then determined after computing the horizontal displacement average jump. The analytical expression is unfortunately too long to be reproduced here. If the matrix phase is assumed to be isotropic and incompressible and in the limit of an infinitely stiff reinforcement, the transverse interaction stiffness reduces to:

κ rigid T = 8φ 2 1 (1 + φ 2 2 )µ 1 b 2 (φ 2 2 -1 -(1 + φ 2 2 ) ln(φ 2 )) (143)
which, once again, has the same behavior when φ 2 → 0 as the permeability estimate derived by Boutin [START_REF] Boutin | Study of permeability by periodic and self-consistent homogenisation[END_REF]:

κ Boutin T = 8(1 + φ 2 )µ 1 b 2 (φ 2 2 -1 -(1 + φ 2 2 ) ln(φ 2 )) (144)

Comparison with a hexagonal unit cell model

The previous analytical estimates based on a CCA model are now compared with respect to a periodic hexagonal microstructure (see material properties in Table 1). Due to the microstructure invariance along the fiber axis direction, the auxiliary problem is solved in generalized plane strain conditions along with periodic boundary conditions, see Figure 4c. The spacing between fibers is equal to 2b, that is the outer diameter of the CCA model, and the fiber radius a is such that the fiber volume fraction φ 2 is equal for both models.

The validity of the CCA model has already been established in previous work but we reproduce in Figure 5a the homogenized longitudinal and transverse shear moduli computed from the resolution of the auxiliary problem (with ∆E = I = 0) and from the lower bounds derived in [START_REF] Hashin | Analysis of Properties of Fiber Composites With Anisotropic Constituents[END_REF]. The values for both models are in very good agreement for a wide range of Matrix (epoxy) Fibers (graphite)

E L = 345 GPa E T = 9.66 GPa E = 3.45 GPa ν L = 0.2 ν = 0.35 ν T = 0.3 µ L = 2.07 GPa µ T = 3.

GPa

Table 1: Material properties of a graphite/epoxy composite (graphite fibers are transverse isotropic with respect to their longitudinal (L) direction, epoxy matrix is isotropic), taken from [START_REF] Hashin | Analysis of Properties of Fiber Composites With Anisotropic Constituents[END_REF].

(a) Longitudinal and transverse shear moduli (b) Longitudinal and transverse interaction stiffness volume fractions, this is also the case for other homogenized moduli. As a consequence, the agreement of both models would also be very good when comparing directly the different coefficients of the D ij partial stiffness tensors.

In addition, the longitudinal κ L and transverse κ T interaction stiffness derived previously are also compared with their numerical estimates using a hexagonal unit cell in Figure 5b. It can be observed that the CCA model and the hexagonal unit cell give very similar values at low volume fractions. Both models differ by a few percents for higher volume fraction, especially regarding the longitudinal interaction stiffness. However, it can be observed that the analytical estimate (133) is much more accurate than the estimates assuming that the inclusion is infinitely stiff (134) or the one from Sudret (135).

multiphase continuum

Illustrative applications of boundary-value problems

In the subsequent subsection, a numerical implementation of the multiphase model has been performed using the FEniCS finite-element library [START_REF] Alnaes | The FEniCS project version 1.5[END_REF][START_REF] Logg | Automated solution of differential equations by the finite element method: The FEniCS book[END_REF]. One key advantage of such models regarding their numerical implementation compared to other generalized continuum media is that no specific finite-element needs to be developed since it is just sufficient to use a two-field formulation with standard displacement interpolation for both fields. In particular, no higher regularity is needed (as in strain gradient models for instance) and no locking behavior is observed.

Reinforced block in compression

For this first example, we revisit a problem already considered in [START_REF] De Buhan | Elastic, Plastic and Yield Design of Reinforced Structures[END_REF][START_REF] Sudret | Modélisation multiphasique des ouvrages renforcés par inclusions[END_REF] which illustrates the capability of the multiphase model to capture size-dependency due to boundary-layer effect. A rectangular block of dimensions L × H in 2D plane strain conditions and reinforced by a regular distribution of N horizontal layers is considered (Fig. 6-left), isotropic material properties of both phases are recalled in Table 2. A downwards vertical displacement of intensity δ is imposed by means of a rigid plate in smooth contact with the block top boundary. The lower boundary is also in smooth contact with a fixed plate while the lateral boundaries are stress-free. Finite-element (FE) computations on the heterogeneous problem have been performed for various numbers of reinforcing layers and will serve as a reference solution.

This problem is also modeled in the multiphase framework by considering identical boundary conditions for both phases (see Fig. 6-right) with the effective properties derived in section 4.4. It has been shown in the previous references that the vertical displacement is the same for both phases and corresponds to a homogeneous vertical strain: U m y = U r y = -δy/H. Conversely, the horizontal displacement is a priori different for both Matrix (phase 1) Reinforcement (phase 2) Volume fraction 1 -η = 0.9 η = 0.1 Young's modulus phases and depends only upon the x-variable. The detailed analytical solution is provided in Appendix B.3. A key feature of this solution is that both phase displacements can be written as a zero-order contribution corresponding to the standard homogenization solution U hom (x) = ν hom, * (δ/H)x plus a first-order (with respect to the scale factor = s/L) correction which enables to satisfy the lateral stress-free boundary condition for both phases. More particularly, this boundary-layer correction is controlled by a characteristic length scale (see (B.12)) depending on the partial stiffnesses and the horizontal interaction stiffness κ x in such a way that it is proportional to the scale factor. For the retained numerical values, it is approximately ≈ 0.75s. The particular structure of this solution justifies a posteriori the scaling used in the two-scale asymptotic expansion of section 3.6.

E 1 = 10 MPa E 2 = 1000 MPa Poisson ratio ν 1 = 0.45 ν 2 = 0.3
Comparison of the obtained analytical solution and the reference FE solution is reported in Figure 7 for various values of the scale factor. Due to the macroscopic invariance along the vertical direction, each symbol of the heterogeneous computations corresponds to a partial average over the corresponding phase in the vertical direction for a given horizontal location. First, both solutions are found to be in perfect agreement. Second, one can observe the influence of the scale factor = 1/N on the boundary-layer effect which manifests itself by a relative displacement between both phases and a drop of the partial stress to zero near the free boundary. Although the characteristic length is ≈ 0.75s, the reinforcement partial stress differs from its constant value in the middle of the sample on a distance corresponding to a few heterogeneities spacing. For sufficiently small scale factors, the central value corresponds to the homogenized solution.

The solution computed in [START_REF] Sudret | Modélisation multiphasique des ouvrages renforcés par inclusions[END_REF] has also been represented in Figure 7. It clearly differs by a few percent from the reference FE computations and the present solution, despite exhibiting the correct trend. This discrepancy is due to two different reasons: first, the phenomenological constitutive equations used by Sudret are valid only for infinitesimally small volume fraction and large stiffness of the reinforcing material; second, the computations relied on an estimate of the longitudinal interaction stiffness κ Sudret x = 8µ 1 /(1 -η)/s 2 based on a rigid pull-out loading of the reinforcing layer which turns out to be 50% smaller than the value (108) obtained from the proposed homogenization procedure when µ 2 µ 1 .

Size-effect of stress concentration near a circular hole

It has been previously demonstrated that, contrary to a standard homogenized medium, the multiphase model can account for size effects due to relative displacements between the fibers and the matrix in the vicinity of free-edges. This example explores the implication on the prediction of stress concentration near a circular hole in an infinite plate. To this end, a quarter of a square plate perforated by a circular hole of radius R in its bottomleft corner is loaded by a uniform uniaxial tension very far from the hole. The plate, in plane-strain conditions, of a fiber-reinforced composite (properties of Table 1) with a fiber volume fraction φ 2 = 10% and with fibers oriented along the x-direction. The far-field uniaxial tension σ ∞ is either aligned with the fibers (longitudinal tension) or perpendicular to it (transverse tension). In the following, the problem will be solved thanks to a FE multiphase model with properties derived in section 4.5 with the CCA model. Different scale factors = s/R will be considered in which s denotes the outer matrix diameter of the CCA model and will be compared to the standard homogenized medium solution. Let us mention that the far-field partial stress boundary conditions have been imposed as σ m nn = φ 1 σ ∞ and σ r nn = φ 2 σ ∞ (with n = x or y depending on the loading case) although the result is insensitive to the specific linear combination as long as σ m nn + σ r nn = σ ∞ . Figure 8a represents the evolution of the horizontal stress near the hole boundary on the vertical line (x = 0, y = R+d) in the case of a longitudinal tension loading. The stress values are very close to the homogenized solution for distances d larger the heterogeneities length scale s. Close inspection of the stress fields near the hole boundary reveals that the multiphase model predicts a larger stress, all the more larger for increasing values of s. The stress field then decreases sharply before reaching the homogenized value from below.

On Figure 8b, the vertical stress near the hole boundary on the horizontal line (x = R + d, y = 0) is represented in the case of a transverse tension loading. In this case, the picture is quite different since both models differ on a larger distance (roughly 2s) but, more interestingly, the vertical stress converges to a well-defined value at the hole boundary, independently of the scale factor. With decreasing scale factors, the stress variation indeed converges to the homogenized solution but not in a uniform manner. More precisely, the stress concentration factor (SCF) predicted by the multiphase model is here approximately 22% smaller than the homogenized value (2.52 against 3.25). Let us mention that a similar behavior is obtained when investigating separately the matrix and fiber partial stresses. Besides, we obtained very similar results when increasing the fiber volume fraction to 50% (SCF = 2.70 for the multiphase model and SCF = 3.28 for the homogenized continuum).

It turns out that the discrepancy between both models occurs on a much larger area when investigating the stress values along the vertical line x = R. The partial vertical stresses along this line have been represented in Figure 9a. It can be seen that both solutions differ on a much larger distance, roughly given by 10s. In Figure 9b, the relative gap ∆ σ = |σ multiphase yy -σ homog. yy |/σ ∞ is used as a measure of the discrepancy between both solutions. Regions where ∆ σ > 5% been represented for three different scale factors, revealing strong anisotropy originating at the hole boundary and extending along the vertical direction on scales much larger than the microscopic scale s. Finally, this picture illustrates that both models differ in critical regions in which delamination or cracking are likely to occur. The correct evaluation of the stress levels in these regions is of paramount importance in order to predict initiation of such phenomena, a task for which multiphase models may prove to be very useful.

Modeling crack bridging and delamination

This last example serves to illustrate the efficiency of multiphase models to simulate complex mechanisms such as crack bridging or delamination in fiber-reinforced materials. Let us consider a rectangular domain of dimension 2L × H consisting of a multilayered medium oriented along direction e x . A vertical crack of length H/2 is introduced in the matrix phase only at middle of the sample (see Figure 10-left). The reinforcing phase is assumed to be perfectly sound and can still transfer loads across this surface. This situation, typical of crack bridging in fiber-reinforced materials cannot be modeled using a standard homogenized medium, except by introducing some ad hoc cohesive forces. On the contrary, the multiphase model can account for such a situation by prescribing different boundary conditions for the matrix and the reinforcing phase along the crack surface. More precisely, one half of the sample is modeled in the multiphase framework (see Figure 10-right) by prescribing zero horizontal displacements for the reinforcing phase along the symmetry plane (x = 0) boundary, whereas zero horizontal displacements for the matrix phase are prescribed only along the uncracked part of the plane (x = 0, y ≥ H/2); on the cracked part, the matrix phase is stress-free. The loading consists of a horizontal displacement (δ = 1) prescribed on the left part to both phases, the vertical displacement remaining free on this boundary. Once again, the solution of the multiphase model will be compared to a fully heterogeneous computation with N = 8, 16 or 32 layers across the sample height, corresponding to a scale factor of s/H = 1/N . Dimensions are fixed to L = 2, H = 1 and material properties are those of Table 2.

Comparing the horizontal displacement on the vertical line x = 0 in the matrix phase for both models shows a very good agreement for all values of the scale factor (see Figure 11a). In particular, the crack-opening of the matrix phase relative to the fixed reinforcements is well reproduced by the multiphase model and decreases (roughly linearly) with the scale factor, eventually converging to zero for a crack bridged by an infinite number of reinforcements. It should be mentioned that the discrete symbols relative to the heterogeneous model correspond to an average displacement over each matrix layer.

In Figure 11b, the interfacial shear stress σ xy on the horizontal line corresponding to the top interface relative to the reinforcement located near y = H/4 is represented (no average for the heterogeneous model in this case). On the same plot, we also represented a reconstructed interfacial shear stress from the multiphase model. Using relation [START_REF] Dumontet | Study of a boundary layer problem in elastic composite materials[END_REF] an infinitesimal reinforcement length, we can indeed relate the horizontal interaction force I x to the top (resp. bottom) interfacial shear stress σ + xy (resp. σ - xy ) as follows:

I x (x) = 1 s (σ - xy (x) -σ + xy (x)) (145)
Using the fact that top and bottom interfacial shear stresses are opposite σ - xy = -σ + xy , we obtain the estimate σ + xy = -sI x /2. It can be observed that this reconstruction is of extremely good quality for all scale factors, even regarding the maximal value of the shear stress at x = 0. The dashed lines originating from the maximal shear stress intercept the horizontal axis at the scale-dependent distance and characterize the length-scale over which the shear stress decreases to zero. Quite interestingly, the maximal shear stress depends only weakly on the scale factor. Once again, this example illustrates that the multiphase model can be used to predict initiation of delamination in the critical region of the crack surface.

As a follow-up to this problem, we also investigated the case where each interface bridging the crack has delaminated over a distance d = 0.1. This problem, with N = 16 layers, has been modeled, on the one hand, by explicitly including a discontinuity in the heterogeneous finite-element mesh and, on the other hand, by forcing the interaction stiffness κ, and thus the interaction force I, to be zero in the region (0 ≤ x ≤ d, 0 ≤ y ≤ H/2) for the multiphase model (see Figure 12). Figure 13a represents each phase horizontal displacement on the vertical line located halfway through the delaminated zone (x = d/2 = 0.05), once again each symbol corresponds to the total displacement averaged over the corresponding phase layer. On Figure 13b, the horizontal partial stresses σ i xx in the sixth layer have been represented along a horizontal line. As before, both stresses and displacements are very well reproduced by the multiphase model, even in this complex case involving a matrix crack and a partially debonded interface as evidenced in Figure 14 by the shear stress field in the vicinity of the debonded region. This figure also indicates some bending of the bottom and top debonded layers which is not accounted for in this version of the multiphase model, possibly explaining the slight discrepancy observed in Figure 13a.

Conclusions and perspectives

Multiphase models have been initially proposed as a generalized continuum capable of accounting for scale effects in media reinforced by stiff linear inclusions in small volume fraction. The reinforced medium is represented as the superposition of two mutually interacting continuous media, each of them possessing its own kinematics along with phenomenological constitutive relations.

In this contribution, a more general formulation has been considered which is not limited to stiff inclusions in small volume fraction. Reminiscent of mixture models for porous media, the macroscopic equilibrium equations have been linked to their microscopic counterpart. A homogenization procedure has then been developed so as to control exactly the averaged quantities involved in this phase-average balance equations. The corresponding auxiliary problem is parametrized by the partial stresses of each phase and an interaction force which can be interpreted as the average traction on the interface separating both phases. In a more convenient manner, this problem can be equivalently formulated by considering phase-wise uniform self-equilibrated eigenstrains and body forces in addition to the standard macroscopic strain. The obtained macroscopic constitutive equations involve, in particular, partial stiffness tensors D ij which can be expressed on the basis of the macroscopic stiffness for biphasic materials and an interaction stiffness tensor κ which can be obtained from the resolution of a permeability-like problem.

Analytical expressions have been derived in the case of multilayered and fiber-reinforced media, recovering, in a unified way, independent results of the literature. Finally, some boundary-value problems have been considered in the multiphase framework. The solution of a multilayered block in compression is perfectly well reproduced with the previously derived effective properties, differing from earlier works which considered the same problem. More complex scale effects induced by boundary layers near free surfaces have been investigated by considering stress concentrations near circular notches. The importance of the discrepancy between the multiphase model and a standard homogenized continuum depends on the loading conditions and the considered region. Non-negligible differences can be reached on regions much larger than the typical heterogeneity size which also happen to be those prone to fracture and delamination. Regarding these two phenomena, it has also been shown that the multiphase approach can model very accurately complex situations such as crack bridging or interface debonding which is very hard to do using traditional continuum models. Investigation of some specific one-dimensional problems also reveals that the multiphase framework generalize in some way shear-lag models and provide a rational basis for deriving their stress-transfer β parameter.

Inspection of the form of multiphase analytical solutions reveals that the phase displacements consist of a zero-order macroscopic displacement which is identical for both phases and is equal to the standard homogenization solution. A first-order correction is governed by the relative displacement between both phases, the amplitude of which is maximal near free boundaries or boundaries with fixed displacement conditions for one phase only. Although the homogenization procedure has been built in a periodic setting, which is not valid near the boundaries, we believe that the multiphase framework provides more degrees of freedom for matching what happens near the boundary with the inner expansion, thereby improving on the standard Cauchy continuum. In addition, the proposed two-scale asymptotic expansion of section (3.6) shows that the multiphase auxiliary problem is consistent with a microscopic description of mechanical fields in terms of partial averages on each phases. In particular, the macroscopic equilibrium equations are equivalent to satisfying the microscopic balance equations on average on each phases, and not only on average on the unit cell as for standard homogenization.

Although interesting results have been obtained regarding stress concentration, bridged cracks or debonded interfaces, it should be acknowledged that scale effects induced by boundary layers are not very strong in a purely elastic setting since the first-order contribution rapidly drops to zero when considering macroscopic quantities such as an apparent modulus for instance. However, situations involving bending of the reinforcing phase can induce much stronger scale-effects, depending on the phase contrast and aspect ratio. A generalization of the proposed procedure to such a situation is therefore a natural perspective. It will however require considering phases with a Cosserat or even micromorphic kinematics, for which proper homogenization procedures still seem lacking at this moment. Besides, the proposed approach can also be applied, in principle, to short fiber-reinforced composites. It seems that, in this case, the effect of relative displacement between both phases is much smaller than for long reinforcements whereas bending effects seem to play an important role. The applicability of the proposed model for such a situation will therefore be the subject of a future work.

When considering material non-linearities such as plasticity or fracture and damage, scale effects may also become much more important. As evidenced by the considered examples, multiphase models offer an extremely natural way to model independently the damaging process of each phase and possible debonding of their common interface. Such aspects are obviously of paramount importance for the assessment of the bearing capacity of many fiber-reinforced materials.

Finally, other topics would be worthy of investigation such as elastodynamics or the extension to geometrical non-linearities. except that A and B are now solutions to the linear system corresponding to boundary conditions at x = L, namely: Recalling that σ m = φ 1 σ 1 , the obtained expression are exactly those of the shear-lag model with β = 1/ . Upon replacing κ by the corresponding longitudinal interaction stiffness of the multilayered medium (108), one has:

1/ 2 = 12 E * φ 1 φ 2 E * 1 E * 2 (t 1 + t 2 ) 2 1/µ = 12 (t 1 E * 1 ) -1 + (t 2 E * 2 ) -1 (t 1 /µ 1 ) + (t 2 /µ 2 ) = β 2 (C.11)
that is the multiphase solution corresponds exactly to the best shear-lag solution with (C.4).

Appendix C.3. Other comparisons

An alternative to the previous microcracked problem corresponds to the case where a prescribed displacement U is imposed to the uncracked phase. In this case, the solution derived in section Appendix B.2.2 reduces here to the same expressions as (C.10) with:

σ 0 E * = U/L 1 + φ 1 E * 1 φ 2 E * 2 L
tanh(L/ ) (C.12)

the denominator of which also appears in the derivation of the effective modulus of a microcracked composites through shear-lag models. Finally, when considering shear lag models for fiber-matrix concentric cylinders, the β parameter appearing in [START_REF] Mccartney | Analytical Models of Stress Transfer in Unidirectional Composites and Cross-Ply Laminates, and Their Application to the Prediction of Matrix/Transverse Cracking[END_REF][START_REF] Nairn | On the use of shear-lag methods for analysis of stress transfer in unidirectional composites[END_REF] is exactly the one corresponding to the longitudinal interaction stiffness derived in (133).

Appendix C.4. Conclusions

As a result, the proposed multiphase models can be seen as a generalization of the shearlag models in the sense that the former degenerate to the latter when investigating purely uniaxial solutions. The proposed homogenization procedure also provide a systematic way for deriving optimal values for the shear-lag β parameter. By considering 3D kinematics and equilibrium equations, multiphase models are therefore able to generalize shear-lag type predictions for complex loadings and damage modes.

Figure 2 :

 2 Figure 2: Auxiliary problem defined on the unit cell consisting of a macroscopic strain E, self-equilibrated eigenstrains η i and body forces ζ(y)I.

Figure 3 :

 3 Figure 3: Multi-layered medium and associated unit cell (left), auxiliary problem corresponding to the longitudinal interaction stiffness κ x identification (right).

  (a) Fiber-reinforced composite (b) Composite cylinder assem-(CCA) model (c) Hexagonal unit cell

Figure 4 :

 4 Figure 4: Unidirectional fiber-reinforced composite and associated microstructure models

Figure 5 :

 5 Figure 5: Comparison of the analytical predictions based on the CCA model and finite-element computations on an hexagonal periodic unit-cell for varying fiber volume fraction

Figure 6 :

 6 Figure 6: Reinforced multilayered block in plane-strain compression: heterogeneous problem (left) and equivalent multiphase model with associated boundary conditions (right)

Figure 7 :

 7 Figure 7: Comparison of matrix and reinforcement displacements (left) and partial stresses the reinforcement (right) for different scale factors = s/L = 1/N obtained with the original multiphase model of [52] (dashed line), the present multiphase model (solid line) and FE computations on the heterogeneous structure (symbols) with N layers.

  (a) Longitudinal tension, stresses along x = 0, y = R + d. Inset: zoom near the hole edge (b) Transverse tension, stresses along x = R + d, y = 0. Inset: zoom near the hole edge

Figure 8 :

 8 Figure 8: Total stress concentration profile near a circular hole for tension loading parallel or perpendicular to the fiber direction for varying values of the scale factor s/R = 0.01, 0.02, 0.05, 0.1 and for the standard homogenized Cauchy medium (black dashed line).

  (a) Vertical partial stresses along x = R, y > 0, matrix stress in thick lines, fiber stress in thin lines. (b) Regions where ∆ σ > 5% for different scale factors, the corresponding microscopic length s is indicated for comparison.

Figure 9 :

 9 Figure 9: Stress concentration along the vertical direction near the hole edge in transverse tension.

Figure 10 :

 10 Figure 10: Left: Tension of a multilayered medium with a matrix crack bridged by reinforcements. Right: Boundary conditions for the associated multiphase model.

Figure 11 :

 11 Figure 11: Matrix displacements along the symmetry plane x = 0 (left) and interfacial shear stress σ xy along a horizontal line near y = H/4 for both models (heterogeneous model : filled symbols, multiphase model : solid lines)

Figure 12 :Figure 13 :Figure 14 :

 121314 Figure 12: Matrix crack bridged by reinforcements with delaminated interfaces: heterogeneous model (left) and multiphase model with zero interaction stiffness in the delaminated zone (right)

(α 11 + 1 E * σ 0 1 - 2 E * σ 0 1 + φ 1 E * 1 φ 2 E * 2 cosh

 1111211122 u r (L) = AL + α 11 + α 21 α hom sinh(L/ )B = U (B.31) σ m (L) = (α 11 + α 12 )A -det α α hom cosh(L/ ) B = 0 (B.32)Resolution of this system gives the following expressions for the displacement and stress fields:u m (x) = U 1 + α L tanh(L/ ) x L -α α 12 + α 22 α 11 + α 21 L sinh(x/ ) cosh(L/ ) (B.33) u r (x) = U 1 + α L tanh(L/ ) x L + α L sinh(x/ ) cosh(L/ ) (B.34) σ m (x) = (α 11 + α 12 ) U/L 1 + α L tanh(L/ ) 1 -cosh(x/ ) cosh(L/ ) + τ m (B.35) σ r (x) = (α 21 + α 22 ) U/L 1 + α L tanh(L/ ) 1 + α 11 + α 12 α 21 + α 22 cosh(x/ ) cosh(L/ ) + τ r (B.36) where α = (α 11 + α 12 )(α 11 + α 21 )/ det α. Appendix B.3. Solution for the multilayered block in compression In the context of the multilayered block problem (section 5.1), u m and u r represent the horizontal component of the partial displacements u i = U i x and σ m = σ m xx and σ r = σ r xx are the uniaxial stresses in the horizontal direction. Since ε 1 yy = ε 2 yy = -δ/H, the partial stresses are given by (B.1)-(B.2) with α ij = D ij xxxx (in this case α 12 = α 21 ) and:τ m = -δ H (γ 11 + γ 12 ) (B.37) τ r = -δ H (γ 21 + γ 22 ) (B.38)whereγ ij = D ijxxyy (note that γ 12 = γ 21 ). The boundary conditions are those of subsection Appendix B.2.1 with t m = t r = 0. Equations (B.21)-(B.22) reduce here to: t w = γ 11 + γ 12 + γ 21 + γ 22 α hom α 12 )(γ 21 + γ 22 ) -(α 21 + α 22 )(γ 11 + γ 12 ) det α (B.40)The first result shows that w(x) = ν hom, * δ H x indeed corresponds to the displacement of the standard homogenized Cauchy medium for which ν hom, * = ν hom /(1 -ν hom ) is the t j = σ 0 /(φ 2 E * 2 ) and expressions of section Appendix B.2.1 particularize to:

Table 2 :

 2 Material properties for the plane-strain multilayered block.

If the microscopic scale is in general well separated with respect to the structural scale in fiber-reinforced

Note that the antisymmetric terms associated with [[a i ]] are due to the fact that the constitutive relations are expressed partially in terms of stress (Σ 1 , Σ 2 ) and of generalized strain (V ). Inverting the last relation to express (Σ 1 , Σ 2 , I) as a function of (E 1 , E 2 , V ) yields a fully symmetric expression.

Note that these estimates are consistent with replacing C hom by its Voigt upper bound C in (81).

It corresponds to the Young's modulus E in the case of an isotropic phase in plane stress state and to E(1 -ν)/(1 + ν) in plane strain conditions.

Homogenized constitutive equations

Let us now introduce the partial stiffness tensors D ij = φ i C i : A ij . Making use of the Hill-Mandel lemma and the Maxwell-Betti reciprocity theorem, the following relations hold true (see Appendix A):

The homogenized constitutive equations therefore admit the following symmetric 3 expression:

From [START_REF] Sudret | Modélisation multiphasique de matériaux renforcés par inclusions linéaires[END_REF], we have the following important result:

C hom = D 11 + D 22 + D 12 + (D 12 ) T (56)

Centro-symmetric case

In the important case of a unit-cell presenting a center of symmetry at its origin, both third-order tensors [[a 1 ]] and [[a 2 ]] vanish. In the remainder of this paper, we will always consider this centro-symmetric case.

The first two constitutive equations regarding partial stresses then decouple from the last one related to the interaction force so that the system of constitutive equations can then be rewritten as:

Appendix A. Symmetry properties of the homogenized constitutive equations

In this section we establish some relations between the constitutive coefficients appearing in the homogenized behavior [START_REF] Stern | Wave propagation in elastic laminates using a multi-continuum theory[END_REF].

We consider the following solution pairs of a particular auxiliary problem (37): (u I , σ I ) associated with problem (P I ) in which E 2 = 0 and I = 0, (u II , σ II ) associated with problem (P II ) in which E 1 = 0 and I = 0 and (u III , σ III ) associated with problem (P III ) in which E 1 = E 2 = 0. For each pair J = I, II, III, we also consider the associated local elastic strain ε el,J .

Recalling that the generalized Hill-Mandel lemma [START_REF] Mccartney | Analytical Models of Stress Transfer in Unidirectional Composites and Cross-Ply Laminates, and Their Application to the Prediction of Matrix/Transverse Cracking[END_REF]:

applies for any pair (u, σ) ∈ SA(I) × KA, the following relations hold true:

Applying now the Maxwell-Betti reciprocity theorem to problems (P I ) and (P II ), we have:

so that (A.5) and (A.6) are in fact equal. We obtain the same result for (A.7) and (A.8) on the one hand and (A.9) and (A.10) on the other hand. Using now the expressions of the localized displacements [START_REF] Nairn | On the use of planar shear-lag methods for stress-transfer analysis of multilayered composites[END_REF] and stresses [START_REF] Soubestre | Non-local dynamic behavior of linear fiber reinforced materials[END_REF] in terms of their loading parameters in both sides of the previous relations, we obtain:

Since these relations hold for any value of the loading parameters, we obtain the relations of (54).

Appendix B. One-dimensional solutions of multiphase continua

We consider here one-dimensional solutions of the equations describing multiphase continua depending on one space variable x in the absence of macroscopic body forces.

Appendix B.1. General structure

Let us assume that the partial uniaxial stresses in both phases are of the following form:

where u m (x), u r (x) are some displacement component for both phases, α ij are generalized stiffness coefficients and τ m , τ r are given constants. The corresponding uniaxial partial equilibrium equations are given by:

where

) is the uniaxial interaction force and κ the corresponding interaction stiffness. Injecting all constitutive relations into the equilibrium equations yields the following coupled second-order differential equations:

These equations can be conveniently rewritten as a linear system of differential equations as follows:

where the prime symbol now denotes the derivative with respect to x. By performing a generalized eigenvalue decomposition of matrices α and κ, one can easily show that the previous system can be written as the following two uncoupled differential equations:

where:

From the last two relations, we deduce the inverse relations:

as well as the partial stresses as functions of w(x) and j(x):

This decomposition makes it clear, at least in the one-dimensional case, that the displacement solutions of multiphase continua can also be expressed as the superposition of w(x) corresponding to the standard displacement of the associated Cauchy continuum since, if u m (x) = u r (x) (j(x) = 0), both fields are then given by w(x), and a part j(x) corresponding to the contribution of the relative displacement between both phases. This last contribution is governed by a modified Helmholtz equation (B.10) characterized by a characteristic length scale which is related to both the partial stiffnesses α ij and the interaction stiffness κ (B.12). It can be seen in the following examples that the amplitude of the relative displacement j(x) scales with /L (L being a structural characteristic length) that is with the heterogeneities scale factor = s/L since ∝ 1/ √ κ and κ ∝ 1/s 2 . Hence, the relative displacement contribution is a first-order correction to the standard homogenized displacement.

Appendix B.2. General solutions for different types of boundary conditions

We now give the expression of the previous system general solution for different types of boundary conditions. The matrix (resp. reinforcement) phase corresponds to phase 1 (resp.

2).

Appendix B.2.1. Fixed displacements at x = 0 and imposed tractions at x = L

We assume here that u m (0) = u r (0) = 0 and that some tractions t m , t r are imposed at x = L. The solution will be more conveniently expressed as a function of the derivatives t w = w (L) and t j = j (L) at x = L. Hence, we will reformulate the stress boundary conditions at x = L using these derivatives by inverting relations (B.13) and (B.14) so that:

that is:

The solution of (B.9)-(B.10) can then be expressed as:

which can be further recombined in terms of partial displacements:

and of partial stresses:

Fixed displacement at x = L for one phase only

We assume here that the displacement at x = 0 is fixed for both phases and imposed to some value U at x = L for one phase only, say the fiber phase u r (L) = U , whereas the other phase is assumed to be stress-free at x = L : σ m (L) = 0. The global displacement w(x) and the jump displacement j(x) are of the same form as before:

macroscopic pseudo-Poisson ratio in plane-strain conditions. The partial stresses in both phases take here the following simple form:

where we recover that total horizontal stress Σ xx (x) = σ m (x) + σ r (x) is identically null for this problem. Finally, it can be shown that the prefactor t j det α α hom corresponds exactly to the partial stress that one would obtain when relocalizing the solution of the standard homogenized Cauchy medium in each phase. However, this stress level is reached only in the central part of the sample and if L.

Appendix C. A link with shear-lag models Shear-lag models have been developed in order to analyze stress states in composite materials and, more particularly, the stress transfer between a fiber embedded in a matrix. Various models have been built using different simplifying hypotheses with respect to the equations of planar or axisymmetric elasticity. We refer to Nairn [START_REF] Nairn | On the use of shear-lag methods for analysis of stress transfer in unidirectional composites[END_REF][START_REF] Nairn | On the use of planar shear-lag methods for stress-transfer analysis of multilayered composites[END_REF] for many references on the subject. One key feature of such models is, for instance, to characterize the heterogeneous axial stress occurring in the case of transverse cracking of the fibers. When subjected to a macroscopic tensile stress, the fiber axial stress is zero at both cracked ends and increases towards the center of a cracked fiber, eventually reaching a homogeneous value corresponding to an uniaxial strain equal to that of the matrix. This occurs if the fiber is sufficiently long compared to the stress transfer length corresponding to the intermediate region on which the axial stress increases from zero to its homogeneous value. As regards the interfacial shear stress, it is maximal at both fiber ends (with opposite sign) and decreases to zero over this stress transfer region. Note that, in the case of matrix cracking with bridging fibers, a similar behaviour is obtained when exchanging the role played by the matrix and the fibers.

Appendix C.1. Shear-lag models for microcracked laminates

Considering the case of a plane strain bilayered laminate under uniaxial stress σ 0 with one layer (phase one) being damaged by the presence of microcracks separated by a distance 2L, shear-lag models predict that the average axial stress in this phase is given by [START_REF] Nairn | On the use of planar shear-lag methods for stress-transfer analysis of multilayered composites[END_REF]:

where E * 1 is the longitudinal modulus of the corresponding phase 5 and where the shear-lag parameter β can be interpreted as the inverse of a stress-transfer length and depends only on the constituent moduli and ply thicknesses t i . Depending on the precise hypotheses of the underlying shear-lag model, various expressions for β have been proposed such as [START_REF] Nairn | On the use of planar shear-lag methods for stress-transfer analysis of multilayered composites[END_REF]:

for [START_REF] Mccartney | Analytical Models of Stress Transfer in Unidirectional Composites and Cross-Ply Laminates, and Their Application to the Prediction of Matrix/Transverse Cracking[END_REF] (C.4)

Nairn derived a more general expression using shape functions which encompasses all these formulas depending on a specific choice of these shape functions. However, a comparison with finite element computations showed that the formula (C.4) proposed in [START_REF] Mccartney | Analytical Models of Stress Transfer in Unidirectional Composites and Cross-Ply Laminates, and Their Application to the Prediction of Matrix/Transverse Cracking[END_REF] is almost exact, yielding an error less than a percent.

Appendix C.2. Solution with the multiphase model

Let us consider here the constitutive equations of the mutiphase model for uniaxial tension in the longitudinal direction. From the constitutive equations expressed in terms of partial compliances, we have: From the general form of the one-dimensional solutions (B.12) with α 11 = φ 1 E * 1 , α 22 = φ 2 E * 2 and α 12 = 0, the characteristic length scale is here given by:

where E * = α hom = φ 1 E * 1 + φ 2 E * 2 . Considering that the matrix phase is cracked at both ends z = ±L and that the fiber phase carries a total tension σ 0 , the following boundary conditions hold true: u m (0) = u r (0) = 0, σ m zz (L) = t m = 0 and σ r zz (L) = t r = σ 0 . One then has t w = σ 0 / E * and