Noël Lahellec
email: lahellec@lma.cnrs-mrs.fr

Pierre Suquet
email: suquet@lma.cnrs-mrs.fr

On the effective behavior of nonlinear inelastic composites: I. Incremental variational principles

Keywords: Composite, homogenization, nonlinear viscoelasticity, elasto-viscoplasticity, variational method, random micro-structure, Hashin-Shtrikman estimates

A new method for determining the overall behavior of composite materials comprised of nonlinear inelastic constituents is presented. Upon use of an implicit timediscretization scheme, the evolution equations describing the constitutive behavior of the phases can be reduced to the minimization of an incremental energy function. This minimization problem is rigorously equivalent to a nonlinear thermoelastic problem with a transformation strain which is a nonuniform field (not even uniform within the phases). In this first part of the study the variational technique of Ponte Castañeda is used to approximate the nonuniform eigenstrains by piecewise uniform eigenstrains and to linearize the nonlinear thermoelastic problem. The resulting problem is amenable to simpler calculations and analytical results for appropriate microstructures can be obtained. The accuracy of the proposed scheme is assessed by comparison of the method with exact results.

Introduction

This study is devoted to the overall response of nonlinear composites comprised of phases which, when deformed, have a partly reversible and partly irreversible behavior. This is the case of nonlinear viscoelastic constituents, where both elastic and viscous effects are always present and coupled, and of elasto-viscoplastic or elasto-plastic constituents for which there exists a domain in stress space where the material behaves in a purely elastic manner, and more generally of most engineering materials [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF].

Reversible effects are associated with a free-energy density which depends only of the state variables of the material. It consists of all the energy available in a given state of the system (at a given material point). The driving forces (or thermodynamic forces) responsible for the evolution of the system, characterized by changes in the state variables, are obtained by derivation of the free-energy with respect to the state variables (Gibbs relations).

On the other hand, irreversibility is characterized by dissipation along an evolution path of the system. Changes in the state variables are governed by the aforementioned driving forces. A convenient (and commonly used) formulation of the corresponding evolution equations for dissipative materials is based on the assumption that they can be derived from a dissipation potential (see [START_REF] Rice | On the structure of stress-strain relations for time-dependent plastic deformation in metals[END_REF][START_REF] Halphen | Sur les matériaux standard généralisés[END_REF]or Germain et al., 1983, for a review). The present study is conducted in this setting of constitutive relations which can be deduced from two potentials, the free-energy function and the dissipation potential.

The problem of predicting the effective response of heterogeneous materials comprised of elasto-plastic or elasto-viscoplastic phases has been mainly addressed in the literature by numerical means (see [START_REF] Christman | An experimental and numerical study of deformation in metal-ceramic composites[END_REF][START_REF] Bao | Particle reinforcement of ductile matrices against plastic flow and creep[END_REF][START_REF] Fish | Computational plasticity and viscoplasticity for composite materials and structures[END_REF][START_REF] Gonzalez | Numerical simulation of elastoplastic deformation of composites: evolution of stress microfields and implications for homogenization models[END_REF], to cite only a few). However numerical simulations are not the objective of the present study but will rather be used to assess the accuracy of the theoretical methods hereafter proposed.

Regarding theoretical or analytical predictions, which are our main objective here, several directions of research have been explored and it is worth mentioning how they differ from the present study.

(1) Elasto-plastic composites are mostly considered in the literature in the framework of a deformation theory in which the actual incremental constitutive relations for the phases are replaced by nonlinear elastic relations. Within this approximation, classical micromechanical schemes, such as the secant method, have been applied (see [START_REF] Chu | Plastic behavior of composites and porous media under isotropic stress[END_REF]Hashin, 1971, Berveiller andZaoui, 1979, among others), as well as more recent homogenization techniques for nonlinear composites with one potential (see for instance [START_REF] Li | The effect of particle shape and stiffness on the constitutive behavior of metal-matrix composites[END_REF][START_REF] Suquet | Effective properties of nonlinear composites[END_REF][START_REF] Bardella | An extension of the secant method for the homogenization of the nonlinear behavior of composite materials[END_REF]and Gonzalez et al., 2004, for particle-reinforced composites). These recent homogenization methods are based on the existence of a single potential and make use of the variational characterization of the local fields and of the effective potential. They include on the one hand rigorous bounds obtained by variational methods, such as the generalization by [START_REF] Talbot | Variational principles for inhomogeneous nonlinear media[END_REF] of the Hashin-Shtrikman technique, the variational method of Ponte Castañeda (1992) and its subsequent interpretation as a modified secant method by [START_REF] Suquet | Overall properties of nonlinear composites : a modified secant moduli theory and its link with Ponte Castañeda's nonlinear variational procedure[END_REF], [START_REF] Suquet | Effective properties of nonlinear composites[END_REF], and other equivalent variational schemes based on convexity inequalities [START_REF] Suquet | Overall potentials and extremal surfaces of power law or ideally plastic materials[END_REF][START_REF] Olson | Bounding the effective yield behavior of mixtures[END_REF]. They also include estimates, which are exact to second-order in the contrast of the phases, obtained by taking a Taylor expansion of the potential [START_REF] Ponte Castañeda | Exact second-order estimates for the effective me-chanical properties of nonlinear composite materials[END_REF][START_REF] Ponte Castañeda | Second-order homogenization estimates for nonlinear composites incorporating field fluctuations. I -Theory[END_REF]. The notion of a linear comparison composite is central in these studies. However, the approximation on which the deformation theory of plasticity is based is legitimate only when the loading is locally proportional (i.e. at all material points in the volume element), the axes of principal stresses remaining fixed and the local stress being monotonically increased with no unloading. This situation is very rarely met in composite materials. The approximation of constant principal directions of stress at each material point at the local level is a reasonable approximation (but by no means an exact result) when the overall loading is itself proportional but it becomes unrealistic in other circumstances when the loading path involves a rotation of the principal axis of macroscopic stress. Similarly, the assumption of monotone loading (at each local material point) is acceptable (but not exact) as long as the macroscopic stress is applied monotonically, but is completely unrealistic under general loading conditions with loading and unloading sequences. These limitations have motivated the present study where the approximation of a deformation theory is not made and where the incremental form of the actual constitutive relations is taken into account.

(2) The Transformation Field Analysis (TFA) proposed by [START_REF] Dvorak | Transformation field analysis of inelastic composite materials[END_REF] and implemented by different authors (Dvorak et al., 1994, Fish and[START_REF] Fish | Computational plasticity and viscoplasticity for composite materials and structures[END_REF] among others) offers an elegant direction of research, in which the incremental form of the constitutive relations is preserved. In this approach, the plastic strain (or more generally the transformation strain which can be due to various effects such as thermal effects, plasticity, phase transformation, damage....) is assumed to be uniform within each individual phase in the composite. However, as is now well-recognized [START_REF] Suquet | Effective properties of nonlinear composites[END_REF][START_REF] Chaboche | Towards a micromechan-ics based inelastic and damage modeling of composites[END_REF], this model leads to excessively stiff predictions (which can be softened by subdividing each phase into subdomains, at the expense of increasing the complexity of the model). This excessive stiffness is due to the assumed uniformity of the transformation strain in each phase which is unrealistic in practice (except in very particular situations such as laminates or nonlinear inclusions in a linear matrix). As will be shown in appendix A, the TFA can be obtained as a special case of the present theory by enforcing the fields of internal variables to be uniform per phase, which is a very crude approximation from which the present study is free.

(3) Linearization of the constitutive relations for nonlinear viscoelastic heterogeneous materials have been proposed by [START_REF] Li | A secant-viscosity approach to the time-dependent creep of an elastic-viscoplastic composite[END_REF], [START_REF] Masson | Self-consistent estimates for the rate-dependent elastoplastic behaviour of polycrystalline materials[END_REF], [START_REF] Brenner | A "quasi-elastic" affine formulation for the homogenised behaviour of nonlinear viscoelastic polycrystals and composites[END_REF], [START_REF] Lévesque | A micromechanical model for nonlinear viscoelastic particle reinforced polymeric composite materials-undamaged state[END_REF] and [START_REF] Brenner | Improved affine estimates for nonlinear viscoelastic composites[END_REF] (among others) to reduce the problem to a simpler one for a linear viscoelastic comparison composite. The present study is also based on an appropriate linearization of the constitutive relations and the two parts of the paper make use of two different types of linearization. But an essential difference with [START_REF] Masson | Self-consistent estimates for the rate-dependent elastoplastic behaviour of polycrystalline materials[END_REF], [START_REF] Brenner | A "quasi-elastic" affine formulation for the homogenised behaviour of nonlinear viscoelastic polycrystals and composites[END_REF] and [START_REF] Brenner | Improved affine estimates for nonlinear viscoelastic composites[END_REF] is that in these previous studies, the stress history is computed at each time step from the initial time t = 0 (and not from the previous time step) to the present time. In other words, assume that the loading is imposed on a time interval [0, T ] discretized into successive time-steps t n , t n+1 ... and that one is interested in the response of the composite at t n+1 , knowing the overall stress and strain at time t n . According to these previous studies, one would have to determine an "equivalent" linear viscoelastic composite at time t n+1 and compute the whole history of the stress and strain from 0 to t n+1 without using the stress and strain determined at time t n (see for instance [START_REF] Lévesque | A micromechanical model for nonlinear viscoelastic particle reinforced polymeric composite materials-undamaged state[END_REF]. By contrast, the present study is a true step-by-step procedure and is radically different from the previous ones.

Here, the strain and internal variables in the composite are determined at time t n+1 using the strain and internal variables at time t n and the loading conditions at time t n+1 . (4) The special case of linear viscoelastic composites, which is a particular case of the class of materials considered here for which the dissipation potential is a quadratic function, deserves a special mention. There exists an abundant literature on the subject, mostly based on the correspondence principle (see for instance [START_REF] Hashin | Complex moduli of viscoelastic composites-I General theory and application to particulate composites[END_REF]. In this approach, the equations governing the local state and the effective properties of linear viscoelastic composites are, after application of the Laplace transform, completely parallel to that of linear elastic composites with complex moduli. Theoretical results and predictive schemes initially developed for elastic composites can be therefore extended to viscoelastic composites by the correspondence principle. One of these theoretical results is that short memory effects in the individual constituents give rise, after homogenization, to long memory effects in the composite [START_REF] Sanchez-Hubert | Sur certains problèmes physiques d'homogénéisation donnant lieu à des phénomènes de relaxation[END_REF][START_REF] Sanchez-Hubert | Sur certains problèmes physiques d'homogénéisation donnant lieu à des phénomènes de relaxation[END_REF][START_REF] Suquet | Elements of Homogenization for Inelastic Solid Mechanics[END_REF]. As a consequence of these long memory effects, the overall stress at time t depends on the whole history of the strain between times 0 and t. The exact effective constitutive relations require the storage at time t of an infinite number of informations corresponding to the overall strains at all previous times. This result is consistent with a more general result of [START_REF] Suquet | Elements of Homogenization for Inelastic Solid Mechanics[END_REF] which shows that the exact effective relations of elasto-viscoplastic composites requires, in general, the knowledge at time t of an infinite number of internal variables. This theoretical result is of limited practical use since a realistic effective constitutive relation should require only the knowledge of a finite (and hopefully limited) number of informations. But this exact result has the merit of highlighting the complexity of the homogenization procedure when the local constituents exhibit both elastic and dissipative effects.

Reducing the size of necessary information cannot be achieved exactly in general, but only within certain approximations. Therefore the question which we are left with is not that of finding an exact theory (too ambitiou)s, but that of finding a set of approximations leading to accurate results and to an easy implementation.

The paper is organized as follows. The constitutive relations of the individual constituents are presented in section 2.1. An incremental variational principle for a time-discretized version of the resulting evolution equations is derived in section 2.2 at the local level and in section 3.2 at the level of the composite.

Variational incremental principles for dissipative, or elasto-plastic systems, have been known for some time but have received for some time only little attention. One of the first occurrence of such a variational principle in a form close to the one that will be used here, can be found in [START_REF] Mialon | Eléments d'analyse et de résolution numérique des relations de l'élasto-plasticité[END_REF] for rateindependent systems (other variational approaches can be found in Simo and Hughes, 1998, chapter 4 and the references herein). [START_REF] Mialon | Eléments d'analyse et de résolution numérique des relations de l'élasto-plasticité[END_REF] established a variational principle for bodies comprised of generalized standard materials and used it to derive several computational algorithms for elasto-plasticity with (or without) hardening. These algorithms (but not the variational principles themselves) were subsequently applied to homogenization problems in plasticity by [START_REF] Marigo | Plasticité et homogénéisation : un exemple de prévision des charges limites d'une structure périodiquement hétérogène[END_REF]. Recently (and independently) variational principles similar in spirit to that of Mialon, but much more general in the sense that they apply also to rate-dependent elasto-viscoplastic materials at finite strains, have been derived by [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF].A renewed interest for incremental variational principles followed this paper, which not only gave new directions for the time integration of rate-dependent or rate-independent dissipative structures, but also opened new perspectives in the understanding of the formation of microstructures in elastoplastic materials [START_REF] Ortiz | Nonconvex energy minimization and dislocation structures in ductile single crystals[END_REF]. The importance of these more recent variational principles in homogenization of nonlinear materials has also been recognized by [START_REF] Miehe | Strain-driven homogenization of inelastic micro-structures and composites based on an incremental variational formulation[END_REF] and [START_REF] Miehe | Homogenization of inelastic materials at finite strains based on incremental variational principles. Application to the texture analysis of polycrystals[END_REF] who made several computational applications of these principles to finite strain plasticity of composites or polycrystals with evolving microstructures.

In the present study, these variational principles are used to derive approximate but almost analytical schemes for the prediction of the global response and first and second-order statistics of the local fields in nonlinear dissipative composites. The method is general and applies both to nonlinear viscoelastic constituents, where elastic and viscous effects are always coupled, and elastoviscoplastic constituents (where there exists a domain in stress space inside which the material behaves in a purely elastic manner). For brevity, the the-ory is illustrated here only in the case of nonlinear viscoelastic materials, but applications to elasto-viscoplastic materials do exist and will be published separately.

The exact and time-discretized effective constitutive relations are recalled in section 3. The time-discretization allows us to reduce the problem with two potentials to the minimization of a single potential for a nonlinear thermoelastic body subjected to nonuniform eigenstrains (not even uniform within each individual phase). The proposed approximation consists in replacing the nonuniform eigenstrains by piecewise uniform (within in each individual phase) eigenstrains and in linearizing the thermoelastic composite. This is done by a procedure close to the modified secant method of [START_REF] Suquet | Overall properties of nonlinear composites : a modified secant moduli theory and its link with Ponte Castañeda's nonlinear variational procedure[END_REF] (equivalent to the variational method of Ponte Castañeda (1992)). The resulting thermoelastic problem is therefore a classical problem for a N-phase composite. The accuracy of the proposed scheme is assessed in section 5 by comparing the predictions of the scheme with exact results for two or threedimensional microstructures.

2 Individual constituents

Generalized standard materials

The composite materials considered in this study are comprised of individual constituents exhibiting a dissipative behavior which can be modelled in the more general framework of constitutive relations deriving from two thermodynamic potentials. The main ingredients for such a model are threefold : a) A finite number of internal variables α which, in addition to the observable strain ε, describe the internal state of the material. Attention is limited here to infinitesimal strains. b) A free-energy function w(ε, α) which gives, for each possible state (ε, α)

of the system (material) the energy available to trigger its evolution. The driving forces associated with the state variables are :

σ = ∂w ∂ε (ε, α), A = - ∂w ∂α (ε, α).
(1) c) A dissipation potential ϕ( α) which relates the evolution of the internal variables α to the driving forces :

A = ∂ϕ ∂ α ( α), or equivalently α = ∂ϕ * ∂A (A), (2)
where ϕ * is the Legendre transform of ϕ. When the two potentials w and ϕ are convex functions of their arguments (ε, α) and α respectively), the material governed by (1) and ( 2) is said to be a generalized standard material [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF]Nguyen, 1975, Germain et al., 1983).

Upon elimination of A between (1) and ( 2), the constitutive relations of the materials under consideration can be re-written as a system of two coupled equations, one of them being a differential equation in time :

σ = ∂w ∂ε (ε, α), ∂w ∂α (ε, α) + ∂ϕ ∂ α ( α) = 0. ( 3
)
Most classical nonlinear viscoelastic or elasto-viscoplastic models can be formulated in this general framework by appropriate choices of the potentials w and ϕ (Germain et al., 1983, Lemaitre and[START_REF] Lemaitre | Mechanics of Solid Materials[END_REF].

Incremental potential

The present section presents the extension to rate-dependent systems of the variational principle of [START_REF] Mialon | Eléments d'analyse et de résolution numérique des relations de l'élasto-plasticité[END_REF] and [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF] in the simplified framework of infinitesimal strains and of a fully implicit timeintegration scheme.

The time derivative in (3) can be approximated by a difference quotient after use of an implicit Euler-scheme. The time interval of study [0, T ] is discretized into time intervals t 0 = 0, t 1 , ...., t n , t n+1 , ..., t N = T . The time step t n+1 -t n is denoted by ∆t (its dependence on n is omitted for simplicity). The timederivative of a function f at time t n+1 is replaced by the difference quotient (f (t n+1 ) -f (t n ))/∆t and the constitutive relations are written at the end of the time step. This time-discretization procedure applied to (3) leads to the discretized system :

σ n+1 = ∂w ∂ε (ε n+1 , α n+1 ), ∂w ∂α (ε n+1 , α n+1 ) + ∂ϕ ∂ α α n+1 -α n ∆t = 0. (4)
Assuming that the fields σ n , ε n and α n are known at time t n , the unknowns σ n+1 , ε n+1 and α n+1 at time t n+1 solve this discretized version of the system (3). It is readily seen by taking the derivative of J defined in (5) with respect to α that the second equation in ( 4) is the Euler-Lagrange equations for the following variational problem (with solution α n+1 ):

Inf α J(ε n+1 , α), J(ε, α) = w(ε, α) + ∆t ϕ α -α n ∆t . (5)
Then, defining the following potential :

w ∆ (ε) = Inf α J(ε, α), (6)
we obtain the following remarkable result which gives the stress as the derivative of a single potential with respect to the strain :

σ n+1 = ∂w ∆ ∂ε (ε n+1 ). ( 7
)
J and w ∆ are respectively called the incremental potential and the condensed incremental potential of the constituent under consideration.

In order to prove that relation ( 7) follows from ( 6) and ( 4), we note that the derivative of w ∆ with respect to ε reads as :

∂w ∆ ∂ε (ε) = ∂J ∂ε (ε, α) + ∂J ∂α (ε, α) : ∂α ∂ε , (8)
where α(ε) denotes the solution of the infimum problem in ( 6). The last term in (8) vanishes by virtue of the stationarity of J with respect to α, and we are left with :

∂w ∆ ∂ε (ε) = ∂J ∂ε (ε, α) = ∂w ∂ε (ε, α). (9)
This relation, applied with ε = ε n+1 gives (7).

Remark 1 : Under the assumption that w and ϕ are convex functions of their arguments, the incremental potential J is a convex function of (α, ε) and the condensed incremental potential w ∆ is a convex function of ε. The latter property follows from a general result [START_REF] Ekeland | Convex Analysis and Variational Problems[END_REF], for instance)

Let Φ(u, p) be a convex function on V × Y (where V and Y are two Banach spaces) with values in R ∪ {+∞} , and let h be defined on Y as

h(p) = Inf u∈V Φ(u, p).
Then h is a convex function.

In the present context, (u, p) = (α, ε), Φ = J and h = w ∆ .

Furthermore, when w and ϕ are twice-differentiable functions, w ∆ is also twicedifferentiable and its second derivative can be expressed after a straightforward calculation as :

∂ 2 w ∆ ∂ε 2 (ε) = L εε -L εα : L αα + 1 ∆t L v -1 : L αε , (10)
where :

L εε = ∂ 2 w ∂ε 2 , L εα = ∂ 2 w ∂ε∂α , L αε = ∂ 2 w ∂α∂ε , L αα = ∂ 2 w ∂α 2 et L v = ∂ 2 ϕ ∂ α2 .
The convexity of w ∆ can alternatively be deduced from (10). Indeed, the convexity of w and ϕ imply (all inequalities should be understood in the sense of positivity for symmetric quadratic forms) :

L εε -L εα : (L αα ) -1 : L αε ≥ 0, L v ≥ 0,
and therefore

L εε -L εα : L αα + 1 ∆t L v -1 : L αε ≥ L εε -L εα : (L αα ) -1 : L αε ≥ 0.
3 Composite materials. Effective behavior

Local problem

A representative volume element (r.v.e.) V of the composite is comprised of N phases occupying domains V (r) with characteristic functions χ (r) and volume fraction c (r) . Each individual phase is governed by the differential equations ( 3) with potentials w (r) and ϕ (r) . The free-energy w and the dissipation potential ϕ at position x are given by:

w(x, ε, α) = N r=1 χ (r) (x)w (r) (ε, α), ϕ(x, α) = N r=1 χ (r) (x)ϕ (r) ( α). ( 11
)
The r.v.e. V is subjected to a path of macroscopic strain E(t) and the local problem which is solved by the local fields σ(x, t), ε(x, t) and α(x, t) reads as :

σ = ∂w ∂ε (ε, α), ∂w ∂α (ε, α) + ∂ϕ ∂ α ( α) = 0 for (x, t) ∈ V × [0, T ], div σ = 0 for (x, t) ∈ V × [0, T ], ε(t) = E(t) + boundary conditions on ∂V.                 
(12) All fields σ, ε, α depend on x and t. The bracket . denotes spatial averaging over V . For definiteness, periodic boundary conditions are imposed on the boundary of V .

The homogenized or effective response of the composite along the path of prescribed strain {E(t), t ∈ [0, T ]} is the history of average stress {Σ(t), t ∈ [0, T ]} where Σ(t) = σ(x, t) .

Effective incremental potential

Upon discretization of the time interval [0, T ], the discretized version of the local problem ( 12) can be written with the help of (7) as :

σ n+1 = ∂w ∆ ∂ε (ε n+1 ), div σ n+1 = 0, in V, ε n+1 = E n+1 + boundary conditions on ∂V.        ( 13
)
where the incremental potentials J and w ∆ read as :

J(ε, α, x) = N r=1 w (r) (ε, α) + ∆t ϕ (r) α -α n (x) ∆t χ (r) (x), w ∆ (ε, x) = Inf α J(ε, α, x).          (14)
The average stress Σ n+1 = σ n+1 satisfies :

Σ n+1 = ∂ w ∆ ∂E (E n+1 ) , (15)
where the effective energy w∆ has the variational characterization :

w∆

(E n+1 ) = Inf ε =E n+1 w ∆ (ε) = Inf ε =E n+1 Inf α J(ε, α) . ( 16
)
To prove the validity of (15), we compute the derivative of w∆ with respect to E :

∂ w∆ ∂E (E n+1 ) = ∂J ∂ε : ∂ε ∂E + ∂J ∂α : ∂α ∂E . (17)
The second term vanishes (stationarity of J with respect to α) and the first term gives, thanks to the relation σ = ∂J /∂ε and Hill's lemma :

σ n+1 : ∂ε ∂E = σ n+1 : ∂ε ∂E = σ n+1 = Σ n+1 .
Thanks to this variational characterization, the homogenization of the evolution problem ( 12) is reduced to the variational problem ( 16). The latter problem amounts to finding the effective energy w∆ of a composite material with one potential w ∆ . Two important things are worth noticing about this condensed potential :

1. The condensed potential w ∆ is not explicitly known and is certainly nonquadratic. 2. The condensed potential w ∆ depends on x not only through the characteristic functions χ (r) (x) but also through the fields α(x) and α n (x). In other words w ∆ (ε) cannot be put under the familiar form (11).

The first point makes the variational problem ( 16) strongly reminiscent of the problem of nonlinear composites composed of N different phases governed by a single potential. However the second point shows that the problem at hand is more complicated than the latter one. Even the full incremental potential J does not depend only on the phase r but also depends on the internal variable α n at time t n which may be completely inhomogeneous, even within each single phase.

Orientation for the rest of the study

Two objectives will be pursued in the rest of this study :

1. First, in order to make contact with recent work on nonlinear composite materials with a single nonquadratic potential, one has to approximate the potentials w ∆ or J which depend on the position through α n (x), by potentials in the form (11). This is more easily achieved with J than with w ∆ . The approximation consists in replacing in the expression ( 14) the field α n (x) in phase r by an effective internal variable α (r) n which remains to be defined. The resulting problem is then close to that for a N-phase nonlinear composite, but not completely identical since there remains an additional minimization over the field α.

A first natural guess for the effective internal variable α (r) n could be to consider the fields α n (x) and α(x) to be piecewise uniform within each individual phase. As shown in appendix A, this choice leads precisely to the Transformation Field Analysis (TFA) of [START_REF] Dvorak | Transformation field analysis of inelastic composite materials[END_REF]. However, it has been observed in the literature [START_REF] Suquet | Effective properties of nonlinear composites[END_REF][START_REF] Chaboche | Towards a micromechan-ics based inelastic and damage modeling of composites[END_REF] that the TFA may lead to inaccurate results in certain circumstances. Therefore the present approach goes beyond piecewise uniform fields of internal variables and accounts for the first and the second moments of these variables in each individual phases. 2. Second, a strategy has to be chosen in order to approximate J in the variational problem ( 16) by a quadratic function. As is known from the study of nonlinear composites with only one potential, several procedures with different degree of complexity and accuracy are available to this effect. Two of them, which have proven to be flexible and efficient in the latter context, have been explored in the present study. The first approach, presented in section 4 derives from the variational method of Ponte Castañeda (1992) and its variant called the modified secant method (Suquet, 1995, Ponte Castañeda andSuquet, 1998). It makes use of an isotropic linear comparison composite. In the present context the linear comparison composite is comprised of isotropic, linearly viscoelastic phases.

The second approach, presented in the second part of this study, is based on a generalization, due to [START_REF] Lahellec | Nonlinear composites: a linearization procedure, exact to second-order in contrast and for which the strain-energy and the affine formulations coincide[END_REF], of the second-order method of Ponte Castañeda (1996). It makes use of an anisotropic linear comparison composite and is presented in the second part of this study.

As will be seen the two above steps, finding an effective internal variable for each phase and linearizing the constitutive relations, are not completely independent.

Working assumptions

For simplicity, several working assumptions will be made in the following. They are assumed to hold for each individual constituent in the composite.

H1:

The internal state variable α is a symmetric, traceless, second-order tensorial variable which can be seen as the inelastic strain.

H2: The free energy w (r) is a quadratic function of ε and α.

For further simplicity it will be assumed that w (r) can be written as :

w (r) (ε, α) = 1 2 (ε -α) : L (r) : (ε -α), (18)
where the fourth-order tensor L (r) has minor and major symmetries :

L (r) ijkh = L (r) jikh = L (r) khij .
H3: The dissipation potential ϕ, which is a convex function of α, is assumed to be a function of the sole invariant αeq :

ϕ (r) ( α) = φ (r) ( αeq ) where αeq = 2 3 α : α 1/2 , ( 19
)
It is further assumed that φ (r) can be written as

φ (r) ( αeq ) = f (r) ( α2 eq ),
where f (r) is a concave function of its scalar argument.

The secant viscosity η sct of the constituent is defined as :

∂ϕ (r) ∂ α ( α) = 2η (r) sct ( αeq ) α, (20)
and can be expressed in terms of the functions φ (r) and f (r) by the relations :

η (r) sct ( αeq ) = 1 3 φ ′(r) ( αeq ) αeq = 2 3 f ′(r) ( α2 eq ). ( 21
)
Remark 2 : Assumption H3 is met by standard potentials. For the sake of simplicity we shall restrict our attention to simple potentials with no account of hardening effects. Consider for instance the power-law potential:

ϕ (r) ( α) = σ (r) 0 ε0 (m + 1) αeq ε0 m+1 , with f (r) (x) = σ (r) 0 (m + 1) εm 0 (x) m+1 2 . ( 22
)
The associated constitutive relations describe a nonlinear viscoelastic behaviour and read as:

ε = M (r) : σ + 3 2 ε0 σ eq σ (r) 0 n s σ eq , M (r) = L (r) -1 , n = 1 m . ( 23
)
The limiting case m = 0 (or equivalently n → +∞) corresponds to rateindependent elastoplasticity with yield stress σ 0 .

Similarly, a simple constitutive model for elasto-viscoplasticity is obtained when the potential ϕ is chosen in the form :

ϕ (r) ( α) = R (r) αeq + σ (r) 0 ε0 (m + 1) αeq ε0 m+1 . ( 24
)
The corresponding constitutive relations read as :

ε = M (r) : σ + 3 2 ε0 (σ eq -R (r) ) + σ (r) 0 n s σ eq , (25)
R (r) is the threshold below which the material behaves elastically and beyond which the strain-rate contains a viscous part.

As can be seen, the two potentials ( 22) and ( 24) have very similar forms although the second one has an elastic domain. Both potentials can be handled in almost the same manner (with the present variational method), except for the numerical implementation where, for the second potential, it is more convenient to work with the inverse of the secant viscosity rather than with the secant viscosity itself. For brevity, the theory will be illustrated in section 5 only for nonlinear viscoelasticity.

Remark 3 : Free-energy potentials more general than (18) should be considered to account for hardening effects (isotropic and kinematic hardening). This is possible within the very same framework but is left for future work.

We would like to concentrate here on the simplest forms of the potentials.

4 Estimation of the effective incremental potential by a variational method

In this section the determination of the effective incremental potential w∆ ( 16) is reduced to a simpler problem involving only a N-phase linear comparison composite by means of a method inspired by the variational procedure of Ponte Castañeda (1992).

Variational procedure

The idea of the variational method is to add and subtract to the original potential J which is difficult to homogenize, an energy J 0 (that for the linear comparison composite) which is more amenable to homogenization, while the difference J -J 0 can still be estimated. Unlike the incremental potential J, the energy J 0 of the linear comparison composite is chosen to be quadratic and piecewise uniform,

J 0 (x, ε, α) = N r=1 J (r) 0 (ε, α)χ (r) (x),
where

J (r) 0 (ε, α) = w (r) (ε, α) + η (r) 0 ∆t (α -α (r) n ) : (α -α (r) n ). (26)
In this expression η (r) 0

and α (r) n are uniform in phase r and will be chosen appropriately in the sequel. Note that J 0 is the incremental potential for a linear viscoelastic comparison composite.

Let ∆J be the difference between J and J 0 :

∆J(α) = N r=1   ∆t ϕ (r) α -α n ∆t - η (r) 0 ∆t (α -α (r) n ) : (α -α (r) n )   χ (r) (x). (27)
Then ( 16) can be re-written as:

w ∆ (E) = Inf ε =E Inf α J 0 (ε, α) + ∆J(α) ≤ Inf ε =E Inf α J 0 (ε, α) + Sup α ∆J(α) . (28)
The last expression, where the supremum of ∆J over α is taken, provides a rigorous upper bound for w ∆ . This upper bound can be too stiff in certain circumstances as reported in [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: a time-integration approach[END_REF].

To improve on this prediction, it is worth noting that Ponte Castañeda and [START_REF] Ponte Castañeda | Variational second-order estimates for nonlinear composites[END_REF] and Ponte Castañeda ( 2002) have observed, in a different but similar context, that sharper estimates can be obtained by requiring only stationarity of the "error function" (here ∆J) with respect to its argument (here α) :

w ∆ (E) ≈ Inf ε =E Inf α J 0 (ε, α) + Stat α ∆J(α) . ( 29
)
We emphasize that the resulting expression has no upper bound character and is only an estimate (hopefully accurate) of w ∆ .

Linear viscoelastic constituents

To show how the method works independently of the nonlinearity associated with a general potential ϕ, let us first begin with the situation where the individual constituents are linearly viscoelastic, following the lines of [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: a time-integration approach[END_REF] where additional results may be found. In this case, the dissipation potential ϕ is a quadratic function :

ϕ (r) ( α) = η (r) α : α, (30)
and the stationarity condition, Stat α ∆J(α), becomes :

Stat α N r=1   η (r) ∆t (α -α n ) : (α -α n ) - η (r) 0 ∆t (α -α (r) n ) : (α -α (r) n )   χ (r) (x).
The solution α of the above stationarity problem satisfies in phase r :

2η (r) α -α n ∆t = 2η (r) 0 α -α (r) n ∆t , (31)
and therefore can be expressed in phase r as :

α(x) = α n (x) -θ (r) α (r) n 1 -θ (r) , (32)
where θ (r) = η (r) 0 /η (r) . With this relation, the last term in (29) can be evaluated and the resulting estimate for w ∆ reads :

w ∆ (E) ≈ Inf ε =E Inf α J 0 (ε, α) + N r=1 c (r) η (r) θ (r) ∆t(θ (r) -1) (α n -α (r) n ) : (α n -α (r) n ) r (33)
where . r denotes the spatial average over phase r and where θ (r) and α (r) n remain to be determined.

4.

2.1 Determination of θ (r) and α (r) n .

The estimate (33) can be optimized with respect to θ (r) and α (r) n . Stationarity of the right-hand side of (33) with respect to these variables yields

∂ ∂θ (r) J (r) 0 (ε, α) + η (r) θ (r) ∆t(θ (r) -1) (α n -α (r) n ) : (α n -α (r) n ) r = 0,
which implies, by virtue of the expression (26) of J (r) 0

:

θ (r) = 1 ± (α n -α (r) n ) : (α n -α (r) n ) r (α -α (r) n ) : (α -α (r) n ) r . ( 34
)
The sign in the above expression is left undetermined at this stage. It turns out that taking the sign + corresponds to solving the problem (28) with a supremum. The sign -corresponds to an infimum and therefore to a rigorous lower bound for the effective condensed potential w ∆ .

Similarly, the stationarity condition over α (r) n reads :

∂ ∂α (r) n J (r) 0 (ε, α) + η (r) θ (r) ∆t(θ (r) -1) (α n -α (r) n ) : (α n -α (r) n ) r = 0 ⇒ α (r) n = α n r + (θ (r) -1) α r θ (r) , (35)
The two relations (34) and ( 35) show that the two unknowns θ (r) and α (r) n depend only on the first and on the second moment of α n (information available from time t n ) and on the first and the second moment of α over the phase r which remain to be determined.

Minimization of J 0 (ε, α).

Once θ (r) and α (r) n are known through the relations (34) and ( 35), the function J 0 is completely specified and can be minimized, first with respect to α and then, after averaging, with respect to ε. The field α, solution of the infimum problem in (33), can be expressed in phase r, after due account of the form (26) of J (r) 0 , as :

α(x) = 2η (r) θ (r) ∆t K + L (r) -1 : K : L (r) : ε(x) + 2η (r) θ (r) ∆t α (r) n , ( 36
)
where K is the fourth-order tensor associated with the projection over deviatoric tensors (we recall that α is traceless), L (r) is the fourth-order tensor of elastic moduli in phase r and ε is the strain field solution of the minimization problem of the incremental potential among all admissible strain fields compatible with an average strain E. Substituting (36) into the expression of

J (r) 0 yields Inf α J (r) 0 (ε, α) = 1 2 ε : L (r) 0 : ε + ρ (r) 0 : ε + f (r) 0 , (37)
where the tensors L (r) 0 and ρ (r) 0 are uniform in phase r and given by :

L (r) 0 = L (r) -L (r) : K : 2η (r) θ (r) ∆t K + L (r) -1 : K : L (r) , ρ (r) 0 = -L (r) : 2η (r) θ (r) ∆t K + L (r) -1 2η (r) θ (r) ∆t α (r) n f (r) 0 = η (r) θ (r) ∆t α (r) n : 2η (r) θ (r) ∆t K + L (r) -1 : L (r) : α (r) n .                        (38)
It is useful to define for further reference1 :

w (r) 0 (ε) = Inf α J (r) 0 (ε, α),
which is the free-energy of a linear thermoelastic phase. Then the relation (33) yields the following estimate for w ∆ :

w ∆ (E) = w 0 (E) + N r=1 c (r) η (r) θ (r) ∆t(θ (r) -1) (α n -α (r) n ) : (α n -α (r) n ) r , (39)
where w 0 is the effective energy of the auxiliary thermoelastic composite

w 0 (E) = Inf ε =E N r=1 c (r) w (r) 0 (ε) r . ( 40
)
The effective response of the composite at time t n+1 , as predicted by the present model, reads as :

Σ = ∂ w ∆ ∂E (E) = ∂ w 0 ∂E (E). ( 41
)
There are apparently more terms in the derivative of w ∆ than the sole derivative of w 0 but these additional terms cancel out by virtue of the stationarity conditions on θ (r) , α (r) n and α.

Remark 4 : Note an important consequence of the relation ( 41), whose validity is not limited to linear viscoelastic constituents: the macroscopic stress can be computed from the macroscopic strain, either by derivation of the condensed potential w ∆ , or by derivation of the effective energy w 0 of the thermoelastic composite. In other words the stress in the actual composite coincides with the stress in the thermoelastic composite.

The general nonlinear case

Let us come back now to the general case where ϕ is a nonquadratic potential.

In order to express the stationarity of the "gap function" ∆J given by ( 27) with respect to α, we note (as in [START_REF] Suquet | Overall properties of nonlinear composites : a modified secant moduli theory and its link with Ponte Castañeda's nonlinear variational procedure[END_REF] that the concavity of f (r) yields :

ϕ (r) α -α n ∆t r = f (r) (α -α n ) 2 eq ∆t 2 r ≤ f (r) (α -α n ) 2 eq ∆t 2 r . (42)
Therefore ∆J r is bounded from above by

∆t f (r) (α -α n ) 2 eq ∆t 2 r - η (r) 0 ∆t (α -α (r) n ) : (α -α (r) n ) r
Stationarity of the above expression with respect to α gives :

2η (r) α -α n ∆t = 2η (r) 0 α -α (r) n ∆t , (43)
where η (r) is the "effective secant viscosity" of phase r,

η (r) = η (r) sct   α -α n (r) ∆t   with α -α n (r) = (α -α n ) 2 eq 1/2 r , (44)
and where the function η (r)

sct is given by the relation (21).

Note that the optimality condition (43) coincides with the optimality condition (31) provided that the η (r) 's are chosen equal to the secant viscosities of the phases. The problem for general nonlinear phases is therefore reduced to that of section (4.2) for linear viscoelastic phases with the noticeable difference that the present η (r) 's depend on the field α which itself depends on the η (r) 's. Not surprisingly, the consistency condition for these two unknowns leads to a nonlinear problem which is solved iteratively.

Local fields in the individual phases

The relation (41) implies that the average stress in the dissipative composite at time t n+1 coincides with the average stress in the linear thermoelastic composite defined by ( 38) :

Σ = n r=1 c (r) L (r) 0 : ε r + ρ (r) 0 = σ , (45)
where ε and σ denote the local strain and stress fields in the auxiliary thermoelastic problem. It is therefore quite natural to approximate the actual stress and strain fields in the actual composite by the same fields in the linear thermoelastic composite.

The way in which the field of internal variables is approximated by the procedure requires a clarification. Two different variational problems have been solved for α leading to two (apparently) different expressions ( 32) and (36) for α. Which one is to be used for a proper approximation of α? For clarity, let α * denote the solution (32) of the stationarity problem Stat α ∆J(α) and α * * the solution (36) of the infimum problem Inf α J 0 (ε, α). Interestingly, it follows from the stationarity conditions for θ (r) and α (r) n , that the first and second moments over each phase of α * and α * * coincide. Indeed, it follows from the expression of ∆J and straightforward algebra that :

∂ ∂α (r) n J (r) 0 (ε, α * * ) + ∆J(α * ) r = 0 ⇒ α * * r = α * r , (46)
∂ ∂θ (r) J (r) 0 (ε, α * * ) + ∆J(α * ) r = 0 ⇒ α * * : α * * r = α * : α * r , (
47) Therefore, θ (r) and α (r) n , which depend only on the first and second moments of α as can be seen from the expressions (34) and ( 35), are defined unambiguously. Consequently there is no ambiguity on the linear thermoelastic composite. The relation (36) can be used to define the approximate field of internal variables α. There remains an ambiguity on the full field α(x), but its first and second moments are defined unambiguously and this is all what matters in order for the scheme to be implemented.

Remark 5 : For further use, an alternative expression of the secant viscosities can be derived by means of the relations ( 46) and ( 47). Indeed, it follows from (32) that

α -α n (r) = θ (r) |1 -θ (r) | α n -α (r) n (r)
.

Therefore

η (r) = η (r) sct   α -α n (r) ∆t   = θ (r) =1 η (r) sct θ (r) |1 -θ (r) |∆t α n -α (r) n (r) . ( 48
)
Remark 6 : θ (r) and α (r) n given by ( 34) and ( 35) depend on the first and second moment of the field of internal variables α, which can be expressed through (36) in terms of the first and second moments of the strain field ε.

Effective response of the composite

The procedure for determining the effective response of the composite can be summarized as follows.

1. At time t n , the first moment α n r and the second moment α n : α n r of the internal variables are known for each individual phase r. 35) and ( 34) are solved for θ (r) and α (r) n by a modification of the hybrid method of [START_REF] Powell | A hybrid method for nonlinear algebraic equation[END_REF]. In the present study use has been made of the MINPACK library (HYBRD1 routine).

The nonlinear equations (

At each step of this iterative procedure the residues of equations ( 35) and ( 34) are obtained by the following procedure :

2.1 The secant viscosities η (r) are known from (48) and the linear thermoelastic composite is known from (38).

2.2

The first and second moments of the strain in the thermoelastic composite are evaluated either by an exact (numerical) computation, or by means of a predictive scheme appropriate for the microstructure of the composite. In the latter case, they can be deduced from the relations ( 69), ( 70) and ( 62). These moments depend on the thermoelastic constants (38).

2.3

The first and second moment per phase of the unknown field of internal variable α at time t n+1 are computed by means of the relations (36) from the first and second moment of the strain field ε computed at step 3.

2.4

The residues are computed using equations ( 35) and (34).

3. Finally, after convergence is reached, the macroscopic stress Σ n+1 can be obtained as the stress in the thermoelastic composite (see remark 4) by taking the average of the microscopic stress field :

Σ n+1 = σ n+1 = N r=1 c (r) L (r) 0 : ε r + ρ (r) 0 . ( 49
)
Details about the rate of convergence of the algorithm are given in section 5.2.

An intermediate model for evaluating the approximation of the "effective internal variable"

As already noted in section 3.3, the reduction from the fully nonlinear incremental potential to the viscoelastic incremental potential (26) involves two separate approximations. The first one is the linearization corresponding to the introduction of the secant viscosity η (r) 0 . The second approximation consists in replacing the inhomogeneous field α n (x) by an "effective internal variable" (EIV) α (r) n in each individual phase. In order to evaluate the influence of each approximation separately, an intermediate model where linearization is the only approximation made, can be proposed. It consists in adding and subtracting to the original incremental potential the following function :

J int 0 (ε, α, x) = N r=1   w (r) (ε, α) + η (r) 0 ∆t (α -α n (x)) : (α -α n (x))   χ (r) (x).
(50) Note that the full field α n (x) enters the expression of this potential whereas it was approximated by α (r) n in (26).

A procedure analogous to that followed in section 4.3 can be developed with the above approximate potential. The stationarity condition ( 43) is replaced by

η (r) 0 = η (r) sct   α -α n (r) ∆t   with α -α n (r) = (α -α n ) 2 eq 1/2 r , (51)
and the expression (36) for the field of internal variables in phase r is replaced by :

α(x) = 2η (r) θ (r) ∆t K + L (r) -1 : K : L (r) : ε(x) + 2η (r) θ (r) ∆t α n (x) . ( 52
)
The free-energy of the thermoelastic composite can be expressed as :

w int 0 (ε, x) = Inf α J int 0 (ε, α, x),
where

J int 0 (ε, α, x) = N r=1 1 2 ε : L (r) 0 : ε + ρ int 0 (x) : ε + f int 0 (x) χ (r) (x), ( 53
)
where the elastic moduli L (r) 0 are uniform in each phase but where the thermal stress ρ int 0 and the initial energy f int 0 are fields given in phase r as :

ρ int 0 (x) = -L (r) : 2η (r) θ (r) ∆t K + L (r) -1 2η (r) θ (r) ∆t α n (x), f int 0 (x) = η (r) θ (r) ∆t α n (x) : 2η (r) θ (r) ∆t K + L (r) -1 : L (r) : α n (x).
Because of the spatial dependence of ρ int 0 , the thermoelastic problem

Inf ε =E w int 0 (ε)
is not a problem for N-phase thermoelastic composite and its solution cannot be expressed analytically. This problem has to be solved numerically with due account of the fact that α n is a field.

This model is therefore of limited practical interest, but it is useful in understanding the error introduced by the "effective internal variable" α (r) n in the "full" model summarized in section 4.5. It will be implemented in section 5.

Applications

In this section, the accuracy of the model presented in section 4 is assessed by comparing its prediction with the exact (computed numerically) response of specific classes of microstructures.

It is important to note that there are two levels of approximation involved in the practical implementation of the present model.

(1) A first approximation (which is in fact two-fold) is made in the relation (29) where the exact effective incremental potential is replaced by an approximate potential involving a secant viscosity η (r) 0

and an effective internal variable α (r) n for each phase. This set of approximations is central in the present study and is the one that we would like to evaluate.

(2) The second approximation lies in the evaluation of the effective energy of the thermoelastic composite (40). This energy can be either evaluated exactly (numerically)in which case there is (almost) no error introduced, or by means of an analytical model in which case an error is introduced if the analytical scheme is not exact for the composite microstructure under consideration.

We would like to concentrate on the first approximation and avoid any additional error which could be introduced at the stage of evaluating the effective energy of the thermoelastic composite. To this aim, two different approaches have been pursued. First, when there is no exact analytical scheme for the microstructure under consideration, the thermoelastic problem can be solved exactly, i.e numerically. This is the approach followed in section 5.1 where a simple three-dimensional microstructure is considered for which no analytical scheme predicting exactly the thermoelastic properties is available. Both the fully nonlinear problem and the thermoelastic problem (40) are solved by the Finite Element Method (FEM).

A second option to avoid approximations in the resolution of the thermoelastic problem, is to consider microstructures for which accurate predictive schemes are available. These microstructures are often complex and constructed iteratively or asymptotically. The microstructures considered in section 5.3 are fair approximations of composite cylinder assemblages for which one of the Hashin-Shtrikman bound can be considered to be an accurate result when the individual phases are linear and isotropic.

The focus here is on two-phase composites composed of inclusions dispersed in a surrounding matrix. The inclusion phase is identified as phase 1, whereas the continuous phase, called the matrix, is identified as phase 2.

5.1 Three-dimensional elastically-reinforced composites 5.1.1 Computational method, microstructure and loading Full-field computations are performed using the Finite Element Method. A simplified volume element consisting of a cylindrical block of matrix containing an inclusion at its center is considered as an approximation to a periodic distribution of inclusions. The applied loading is a uniaxial tensile stress :

Σ(t) = Σ 33 (t)e 3 ⊗ e 3 . ( 54
)
The problem is axisymmetric and solved using axisymmetric finite elements on the quarter of a vertical section of the unit cell (figure 1b). The boundary conditions are symmetry conditions (zero shear on all sides of the domain, zero normal displacement on z = 0 and r = 0 and uniform horizontal and vertical displacements α and β on the sides r = 1 and z = 1 of the domain. The vertical displacement β is increased at a constant rate β = Ė33 = 10 -2 s -1 , whereas the horizontal displacements α is adjusted iteratively within each time-step to ensure that the average stress in the unit cell is a uniaxial tension in the form (54).

The inclusions are linearly elastic, whereas the matrix is nonlinear viscoelastic with a power-law dissipative potential in the form ( 22). The following material data were used :

-Fibers :

E = 400 GPa, ν = 0.2. -Matrix : E = 70 GPa, ν = 0.3, σ 0 = 480 MPa, ε0 = 10 -2 s -1 .
The rate-sensitivity exponent m of the matrix was varied from 0.1 to 1.

Discussion of the results

The response of the microstructure of figure 1 under monotone uniaxial tension is shown in figure 2a, b and c for m = 1, 0.2 and 0.1 respectively. The overall axial stress is plotted as a function of the overall axial strain. Four sets of results are shown. The exact nonlinear simulations are shown as solid lines labelled "Exact(FEM)". The predictions of the intermediate model of section 4.6 where the thermoelastic problem is solved "exactly" using the FEM without making the approximation of an "effective internal variable" are shown as dashed lines and labelled "SEC(FEM)". The predictions of the full model summarized in section 4.5 where in addition to the secant linearization, the approximation of the field α n by an "Effective Internal Variable" (EIV) in each phase is made, are shown as dotted lines and labelled "SEC(EIV+FEM). Finally, the effective energy of the thermoelastic composite can be estimated using the lower Hashin-Shtrikman bound for isotropic composites. More specifically use is made of the relation (68) in appendix B, where the effective tensor L is given by the appropriate Hashin-Shtrikman (HS) bound. The corresponding prediction is shown in long-dashed lines and labelled as "SEC(EIV+HS)".

The following observations can be made :

(1) All curves can be schematically decomposed into three different regimes.

The first regime corresponds to the initial response of the composite which is purely elastic. In this regime the only potential which matters is the free-energy and the effective response of the composite can be obtained by homogenization of a single quadratic potential. A good match in this first regime indicates that the elastic properties of the composite are accurately predicted. In the third part of the curve, corresponding to (relatively) large deformations, the response of the composite is purely dissipative and therefore governed by the homogenized dissipation potential of the composite, when elastic effects are neglected. In between those two regimes there is an intermediate, or transient, regime corresponding to the zone where both elastic and dissipative effects are coupled. This regime has precisely motivated the development of the present model. Obviously, the agreement between the model and the exact results cannot be satisfactory in this regime if it is not satisfactory in the other two regimes. As will be seen, the error in this regime is indeed smaller than in the asymptotic regime.

(2) The error due to the EIV approximation, measured as the difference between SEC(FEM) and SEC(EIV+FEM), is seen to be small at all ratesensitivity exponents. Therefore the replacement of the full field α n by an effective internal variable α (r) n , after due account of the field fluctuations, is legitimate. The great advantage of this substitution is that it leads to a variational problem for a N-phase thermoelastic composite, which is, often, amenable to analytical resolution.

(3) The error introduced by the linearization step, reflected for instance in the difference between the "exact" curve Exact(FEM) and the "secant" predictions SEC(FEM), leads to an overestimation of the response of the composite. This is due to the inequality in ( 42). (4) Use of the HS lower bound, which leads to an underestimation of the effective energy of the linear thermoelastic composite, tends to compensate the error due to the linearization scheme. A similar compensation of errors has already been observed for nonlinear composites governed by a single potential (see for instance [START_REF] Moulinec | Intraphase strain heterogeneity in nonlinear composites: a computational approach[END_REF]. ( 5) The exact response of the composite shows that, roughly speaking, the size of the transient regime (where both elastic and dissipative effects are coupled) depends strongly on the rate-sensitivity exponent m and increases with increasing m.

Influence of the time-step

To evaluate the influence of the time-step on the accuracy of the predictions, the model SEC(EIV+HS) has been used with different time steps ranging from 0.01s to 0.2s (the macroscopic strain-rate was 10 -2 s -1 ). The corresponding predictions of the present model in the transient regime are shown in figure 2d. As can be seen the time-step has very little influence on the results and no influence on the asymptotic response.

Two-dimensional composites

Microstructures

The class of microstructures considered in this section approaches the composite cylinder assemblage of Hashin for which one of the Hashin-Shtrikman bound is an accurate approximation (it gives an exact result for the effective bulk modulus and a sharp estimate for the shear modulus). It has previously been used by [START_REF] Moulinec | Intraphase strain heterogeneity in nonlinear composites: a computational approach[END_REF], [START_REF] Moulinec | Homogenization for nonlinear composites in the light of numerical simulation[END_REF] and [START_REF] Idiart | Macroscopic behavior and field fluctuations in viscoplastic composites: second-order estimates versus full-field simulations[END_REF] in studies assessing the accuracy of different homogenization methods for nonlinear composites governed by a single potential. It consists of the assemblage of self-similar cylinders generated by a method described in [START_REF] Moulinec | Intraphase strain heterogeneity in nonlinear composites: a computational approach[END_REF]. Ideally an infinite number of sizes would be required to achieve exactly a composite cylinder assemblage so as to fill-in the whole space. To keep the computational complexity of the problem accessible, only three different sizes of cylinders have been considered. Two different volume fractions c (1) = 0.21 and c (1) = 0.41 have been investigated. Using results in [START_REF] Moulinec | Intraphase strain heterogeneity in nonlinear composites: a computational approach[END_REF] or [START_REF] Idiart | Macroscopic behavior and field fluctuations in viscoplastic composites: second-order estimates versus full-field simulations[END_REF] only one configuration has been considered for each volume fraction, the configuration which appeared to be the more representative of the average of the results obtained in those previous studies. The corresponding microstructures are shown in figure 3.

Material data

The individual constituents are nonlinear viscoelastic with a power-law dissipation potential in the form ( 22) having the same exponent in both phases. The elasticity of each phase is isotropic and characterized by a Young's modulus E and a Poisson ratio ν. The material data are as follows :

-Fibers : E = 100 GPa, ν = 0.45 σ 0 = 5 (hard fibers) or σ 0 = 0.2 GPa (weak fibers), ε0 = 10 -2 s -1 , -Matrix : E = 100 GPa, ν = 0.45, σ 0 = 1 GPa, ε0 = 1 s -1 . The rate-sensitivity exponent m varies from 1 to 0.1 and is identical in both phases.

Loading conditions

The loading applied to the representative volume element is an imposed macroscopic in-plane shear strain :

E(t) = E 11 (t) (e 1 ⊗ e 1 -e 2 ⊗ e 2 ) , (55)
with a constant strain-rate Ė11 = 5 10 -1 s -1 .

Computational methods

Exact results: The method used in the full-field numerical simulation is based on Fast Fourier Transforms (see [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear properties of composites[END_REF] and [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] for a presentation of the method). More details about the application of this method to nonlinear viscoelastic composites can be found in [START_REF] Idiart | Macroscopic behavior and field fluctuations in viscoplastic composites: second-order estimates versus full-field simulations[END_REF] section 4.

Proposed approximate model: Use has been made of the algorithm described in section 4.5. The effective properties of the linear comparison thermoelastic composite are obtained using the relevant Hashin-Shtrikman estimate. As shown in [START_REF] Moulinec | Intraphase strain heterogeneity in nonlinear composites: a computational approach[END_REF] this bound is an accurate, but not exact, prediction of the linear elastic properties of the microstructures shown in figure 3.

A major consequence of the use of the Hashin-Shtrikman prediction, is that the strain field, as estimated by this model, does not fluctuate in the inclusion phase. Within this model, the field of internal variables α is therefore uniform in the inclusion phase at all times :

α =< α > 1 , α n =< α n > 1 , in phase 1.
According to its expression (44), the secant viscosity of the inclusion phase reduces to

η (1) = η (1) sct (< α > 1 -< α n > 1 ) eq ∆t .
The equations ( 34) and ( 35) reduce to :

θ (1) = 1, α (1) n =< α n > 1 .
Consequently, the only true unknowns in the problem are θ (2) , α (2) n and < α > 1 .

In all examples investigated in the present study, the iterative algorithm of section 4.5 has always converged to a solution. We believe that this solution is unique (at least in rate-dependent problems) but we have no proof of this assertion. Convergence is reached when the relative errors between two iterates for θ (2 ), α (2) n and < α > 1 is less than a given threshold (typically 10 -8 ). With this criterion and for the examples of section 5.3, m = 0.2, c (1) = 0.41, convergence was attained in 20 to 50 iterations (typically), the slowest convergence being observed in the transient regime where elastic and viscous effects are of the same order of magnitude.

Discussion of the results

The responses of the two microstructures of figure 3 under monotone shear deformation are shown in figures 4 and 5 where the overall stress is plotted as a function of the overall strain. The full-field simulations are shown as solid lines whereas the predictions of the present model are shown as dashed lines.

As in the three-dimensional case, the main trends of the stress-strain curve for these composites under monotone deformation at constant strain-rate are typically an initial linear-elastic response for incipient strains, then a transient part and finally a plateau corresponding to the purely viscous response of the composite. This plateau is characterized by a stress Σ ∞ 11 .

Figure 4 corresponds to the case of "hard" fibers (σ 1 0 /σ 2 0 = 5), whereas figure 5 corresponds to "weak" fibers (σ 1 0 /σ 2 0 = 0.2). In both figures two rate-sensitivity exponents are considered, m = 1 (corresponding to linear viscoelasticity) and m = 0.2. The following comments can be made :

(1) In the linear case (m = 1) the error introduced by the model is essentially due to the HS bound. There is an additional error introduced by the EIV approximation, when the heterogeneous field α n is replaced by a piecewise uniform effective internal variable, but this additional error is small as already seen in section 5.1. (2) In all cases the largest error is found in the asymptotic regime. It is therefore possible to better quantify the error introduced by the model by comparing the asymptotic stresses Σ ∞ 11 . This is done in figure 6 where the asymptotic stress is plotted as a function of the nonlinearity exponent. The accuracy of the model in the asymptotic regime is essentially dictated by the choice of the linearization scheme for the dissipation potential ϕ and is not much influenced by the other approximations made by the model (time-discretization, stationarity conditions, effective internal variable).

(3) Regarding the linearization technique, it is now well-recognized that the "modified secant method" used here (which is a special case of the variational procedure of Ponte Castañeda) has the advantage of simplicity but can lead to significant errors for large contrast between the phases or high volume fraction of the constituents. This is especially true for fiber-weakened composites as can be seen in figure 6c andd. The error remains acceptable for fiber-reinforced composites, but is mainly due to a "compensation of errors" between the underestimation brought by the HS lower bound and the overestimation due to the linearization scheme itself (contained in inequality ( 42)). Improving on the "secant linearization" is the objective of part II of this study (companion paper).

Conclusion

This paper is devoted to the effective response of composites, or more generally of heterogeneous materials, whose individual constituents exhibit both reversible and irreversible effects in their mechanical behavior. The main findings are:

(1) Upon time-discretization of the constitutive differential equations, an incremental variational principle has been derived, by means of which the problem can be reduced to the minimization of a single nonquadratic incremental potential. (2) The spatial dependence of the incremental potential stems not only from the variation of the potential from phase to phase, as is naturally the case in heterogeneous systems, but also from the dependence of this incremental potential on the field of internal variables at the previous time-step which is nonuniform even within each phase. (3) A strategy, inspired by the variational procedure of Ponte Castañeda, has been used to linearize the nonquadratic condensed potential and also to define an "effective internal variable" per phase at each time-step. (4) Comparisons with full-field simulations show that the present model is good as long as the variational procedure is accurate in the purely dissipative setting, when elastic deformations are neglected. If this is the case, the present model accounts in a very satisfactory manner for the coupling between reversible and irreversible effects and is therefore an accurate model for treating nonlinear viscoelastic and elasto-viscoplastic materials.

(5) In certain situations, the variational procedure is not accurate in the purely dissipative limit and examples of such situations are given. They motivate the second part of this study in which a more refined scheme, still based on the condensed variational potential derived in this first part, but with a different linearization strategy based on an anisotropic linear viscoelastic composite, is proposed.

within each phase :

α(x, t) = N r=1 α (r) (t)χ (r) (x). ( 56)

In particular :

α n (x) = N r=1 α (r) n χ (r) (x), (57)
with a similar expression for α n+1 (for simplicity the index n+1 will be omitted in the sequel). The restriction brought by ( 56) on the choice of the trial fields α leads to an upper bound for w ∆ and therefore to an upper bound for the effective incremental potential J at each time step. The resulting prediction will be stiffer than the actual response of the composite.

The stationarity conditions with respect to α (r) in the variational problem (16) read as :

∂w (r) ∂α (ε(x), α (r) ) + ∂ϕ (r) ∂ α α (r) -α (r) n ∆t χ (r) (x) = 0, ( 58
)
while the conditions expressing the stationarity with respect to ε read as :

σ(x) = ∂w (r) ∂ε (ε(x), α (r) ) in phase r. (59)
It follows from (58) and from the fact that α (r) and α (r) n are independent of x, that :

A (r) = ∂ϕ (r) ∂ α α (r) -α (r) n ∆t = - ∂w (r) ∂α (ε(x), α (r) ) r .
As a consequence of assumption (H2), σ is a linear function of ε and α with uniform moduli in phase r. Therefore, the average stress σ (r) over phase r is related to the average strain ε (r) and to α (r) by the same relation. The final set of evolution equations for the model obtained under the specific assumption (56) reads as : r) , α (r) ), A (r) = -∂w (r) ∂α (ε (r) , α (r) ). ( 60)

σ (r) = ∂w (r) ∂ε (ε (
A (r) = ∂ϕ (r) ∂ α ( α(r) ), or equivalently α(r) = ∂ϕ (r) * ∂A (A (r) ), (61) This is exactly the set of equations used in the TFA.

Fig. 1 .

 1 Fig. 1. Three-dimensional microstructure. (a): Periodic unit-cell (right) approximated by a cylindrical circular unit-cell (left). (b) Two-dimensional mesh used in the computations.

Fig. 2 .

 2 Fig. 2. Axisymmetric unit cell. Elastic particles in an nonlinear viscoelastic matrix. Particle volume fraction c (1) = 0.25. Effective response of the composite for (a): m = 1, (b): m = 0.2, (c): m = 0.1.(d): Sensitivity of the predicted response to the time-step ∆t for m = 0.1.

Fig. 3 .

 3 Fig. 3. Composite cylinder assemblages. (a): c (1) = 0.21, (b) c (1) = 0.41.

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4. Response of fiber-reinforced (hard fibers) nonlinear viscoelastic composites under monotone loading for different volume fractions and different nonlinearity exponents m.(a): c (1) = 0.21, m = 1.(b): c (1) = 0.21, m = 0.2. (c): c (1) = 0.41, m = 1.(d): c (1) = 0.41, m = 0.2.

Note that, in order to avoid confusions, the present notations are slightly different from those used in[START_REF] Lahellec | Effective behavior of linear viscoelastic composites: a time-integration approach[END_REF]. The free energy of the thermoelastic composite, denoted here by w (r) 0 , was denoted by w (r) ∆ in(Lahellec and Suquet,

2006).

Acknowledgements

This study is part of program supported by the Centre National de la Recherche Scientifique under grant CNRS/NSF n • 14555. The authors thank R. Masson for refreshing their mind about the work of Mialon (1986).

Appendix A: Link with the Transformation Field Analysis (TFA)

In an attempt to approximate w ∆ by a piecewise uniform potential, the trial field α involved in the variational problem ( 16) can be chosen to be uniform Appendix B: Effective thermoelastic energy for N-phase and twophase composites Consider a N-phase thermoelastic composite whose phase r is characterized by the energy r) .

Its effective energy reads [START_REF] Willis | Variational and related methods for the overall properties of composites[END_REF] :

where :

These expressions simplify further for two-phase composites (N = 2). According to Levin's relations [START_REF] Levin | Thermal expansion coefficients of heterogeneous materials[END_REF], the localization operators can be explicitly expressed in terms of the effective stiffness L of the composite:

It follows from these relations that the effective thermoelastic energy of the composite can be expressed in terms of the effective stiffness of the composite only :

where

Appendix C: First and second moments of the strain field in a thermoelastic composite

The first moments of the strain field ε in each phase of a linear thermoelastic composite can be expressed as (see for instance [START_REF] Willis | Variational and related methods for the overall properties of composites[END_REF] : r) . ( 69) r) is the strain localization tensor giving the average strain in phase r due to the macroscopic strain E in the absence of thermal stresses (ρ (r) = 0). Similarly, a (r) is the average strain in phase r due to the thermal stresses (E = 0). Both tensors depend on the elastic moduli of the phases and on the microstructure of the composite. Explicit expressions corresponding to specific microstructures can be found in the literature. These tensors can also be evaluated numerically.

The second moment of the strain field is obtained by derivation of the effective energy of the composite with respect to the elastic moduli L (r) of the individual phases (see for instance Ponte Castañeda and Suquet, 1998[START_REF] Buryachenko | Multiparticle effective field and related methods in micromechanics of composite materials[END_REF][START_REF] Ponte Castañeda | Second-order homogenization estimates for nonlinear composites incorporating field fluctuations. I -Theory[END_REF] :

Detailed expressions for the effective energy in two-phase composites are given in appendix B.