Renald Brenner

Renaud Masson

Olivier Castelnau

André Zaoui

A "quasi-elastic" affine formulation for the homogenised behaviour of nonlinear viscoelastic polycrystals and composites

Keywords: Nonlinear viscoelasticity, Polycrystal, Homogenisation

Introduction

During the last decades, several attempts have been made to predict the overall response of composites or polycrystals from that of their constituent phases when the latter obeys a nonlinear viscoelastic (or rate-dependent elastoplastic) behaviour. After the first approximate treatments proposed in the seventies by [START_REF] Brown | A self-consistent polycrystalline model for creep under combined stress state[END_REF] or [START_REF] Zaoui | Effets de la désorientation des grains sur la viscoplasticité de polycristaux C[END_REF], several developments for small strains [START_REF] Weng | Self-consistent determination of time dependent behavior of metals[END_REF] or finite strains [START_REF] Nemat-Nasser | Rate-dependent finite elstoplastic deformation of polycrystals[END_REF][START_REF] Harren | The finite deformation of rate-dependent polycrystals-II. A comparison of the self-consistent and Taylor methods[END_REF] were based on an extension of Kröner's conception of the self-consistent approach for elastoplasticity [START_REF] Kröner | Zur plastischen Verformung des Vielkristalls[END_REF], according to which the viscoplastic strain is considered as a stress-free strain in the sense of Eshelby. It has been proved [START_REF] Zaoui | Investigation of inter-and intragranular plastic inhomogeneity[END_REF]) that, such as for elastoplasticity, this approach still deals with elastic interactions between the phases and leads then to too stiff overall responses.

An important specificity with respect to elastoplasticity lies in the fact that the viscoelastic coupling is responsible for a complex "long range memory effect" [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF] which, for instance, shows itself in the non-Maxwellian overall behaviour 1 of a two-phase material with Maxwellian constituents. Various propositions have been made to deal with this characteristic feature. Generally, they are based either on ad hoc concentration equations [START_REF] Weng | A self-consistent relation for time-dependent creep of polycrystals[END_REF][START_REF] Turner | Self-consistent modelling of nonlinear visco-elastic polycrystals: an approximate scheme[END_REF][START_REF] Molinari | On the self-consistent modeling of elastic-plastic behavior of polycrystals[END_REF] or on an internal variables procedure [START_REF] Paquin | Integral formulation and self-consistent modelling of elastoviscoplastic behavior of heterogeneous materials[END_REF]. To extend to the nonlinear context the pioneering work of [START_REF] Laws | Self-consistent estimates for the viscoelastic creep compliances of composite materials[END_REF] proposed for linear viscoelasticity, an alternative approach based on the use of the correspondence principle [START_REF] Mandel | Mécanique des Milieux Continus[END_REF] and the Laplace transform technique has been proposed by [START_REF] Masson | Self-consistent estimates for the rate-dependent elastoplastic behavior of polycrystalline materials[END_REF]. This theory, which extends to rate-dependent elastoplasticity the tangent linearisation treatment of the local constitutive equations introduced by [START_REF] Molinari | A self-consistent approach of the large deformation polycrystal viscoplasticity[END_REF] for viscoplasticity, makes use of a homogenisation procedure proposed by [START_REF] Rougier | Self-consistent modelling of elastic-viscoplastic polycrystals[END_REF]. The resulting so-called "affine" model retains the viscoelastic nature of the interactions between the constituents and the inherent long range memory effect. For linear behaviours, it is consistent with the previous results of [START_REF] Laws | Self-consistent estimates for the viscoelastic creep compliances of composite materials[END_REF]. For nonlinear behaviours, it tends in the steady-state regime to the affine procedure for viscoplasticity [START_REF] Masson | An affine formulation for the prediction of the effective properties of nonlinear composites and polycrystals[END_REF]. Nevertheless, although tractable, it has to make use of a complex and computer time-consuming numerical resolution of the inverse Laplace transformation through a collocation method. This method has already been promoted by [START_REF] Schapery | Approximate methods of transform inversion for viscoelastic stress analysis[END_REF] (and also frequently applied in linear viscoelasticity, e.g., [START_REF] Turner | Self-consistent modeling of visco-elastic polycrystals: application to irradiation creep and growth[END_REF]) who proposed at the same time a powerful approximation referred to as the "direct method".

The main practical aim of this paper is a systematic application of this last method to the affine modelling of the overall behaviour of nonlinear viscoelastic polycrystals and composites. We first (Section 2) recall the essential points of the affine procedure for this class of materials. In Section 3, we analyse different ways to make use of the "direct method" for homogeneous and heterogeneous linear viscoelasticity and different applied strain and stress paths. We then propose (Section 4) a "quasi-elastic" approximation of the affine formulation for nonlinear heterogeneous viscoelasticity and we show in Section 5 that it provides a very satisfying tool for the prediction of the creep and relaxation behaviour of polycrystals. Further progress are briefly discussed in the conclusion.

The affine procedure for nonlinear viscoelastic materials

The linearisation procedure

Homogenisation of the nonlinear behaviour of heterogeneous materials generally resorts to the linearisation of the local constitutive equations before applying some scale change scheme. This is also the case for nonlinear viscoelasticity. For such a behaviour, variational procedure (see Ponte [START_REF] Ponte Castañeda | The effective mechanical properties of nonlinear isotropic composites[END_REF][START_REF] Ponte Castañeda | Exact second-order estimates for the effective mechanical properties of nonlinear composite materials[END_REF][START_REF] Suquet | Overall potentials and extremal surfaces of power law or ideally plastic composites[END_REF] aiming at the construction of bounds or estimates cannot be invoked since the constitutive equations do not derive from a single potential, neither at the local nor at the global scale. In this context, for which only estimates can be investigated, an alternative approach allowing to describe the coupling between conservative (elastic) and dissipative (viscous) regimes has to be adopted. The affine procedure [START_REF] Masson | Self-consistent estimates for the rate-dependent elastoplastic behavior of polycrystalline materials[END_REF], developed in that goal, has been chosen for the present study.

Let us consider a multiphase nonlinear viscoelastic medium, occupying a volume Ω, whose phases obey constitutive equations which read, for infinitesimal strains, at any time t

ε(x, t) = s(x) : σ (x, t) + g x, σ (x, t) , ∀x ∈ Ω, (1)
with ε and σ the strain and stress rates, s the elastic compliance and g the nonlinear viscous (or viscoplastic) strain rate deriving from the potential u. For the sake of simplicity (but this is by no means a limitation), a sole dependence of g on the local stress state σ (x, t) is considered:

g x, σ (x, t) = ∂u ∂σ x, σ (x, t) . ( 2
)
The material is supposed to be subjected to a known macroscopic mechanical loading for t ∈ [0; +∞[. For brevity, the position vector x is omitted afterwards.

The main difficulty of the linearisation step is related to the differential nature of relation (1). Indeed, unlike elastoplasticity or viscoplasticity, it does not allow us to define directly instantaneous moduli or compliances since the strain rate at each time t depends on the whole stress history from t = 0. Consequently, associated to the linearisation of the strain rate ε at a prescribed time t = ζ , it is necessary to know its expression from t = 0 to t = ζ . This linearisation, proposed by [START_REF] Rougier | Self-consistent modelling of elastic-viscoplastic polycrystals[END_REF], reads

ε(t) = s : σ (t) + m(ζ ) : σ (t) + ε0 (t, ζ ), m(ζ ) = ∂g ∂σ (ζ ), (3) with ε0 (t, ζ ) = g(σ (t)) -m(ζ ) : σ (t) ∀t ζ, ε0 (t, ζ ) = g(σ (ζ )) -m(ζ ) : σ (ζ ) ∀t ζ (4)
which is an approximation of a more general procedure proposed by [START_REF] Pouya | Linéarisation et homogénéisation en viscoélasticité[END_REF]. It can be noted that the stress-free strain rate ε0 (t, ζ ) does not depend on the mechanical loading after t = ζ . With use of the notion of the "linear composite comparison medium" (Ponte Castañeda, 1991), we can say that according to (3) and ( 4), this procedure refers to a linear thermoviscoelastic heterogeneous medium. Relation (3) can be put in the more convenient form of a Stieljes convolution product:

                       ε(t) = t 0 q ζ (t -u) : dσ (u) + ε 0 (t, ζ ), q ζ (t) = s + m(ζ )t, ε 0 (t, ζ ) = t 0 ε0 (θ, ζ ) dθ.
(5)

The stress-free strain ε 0 is a function of the whole stress history thus accounting for the viscoelastic coupling phenomena within the material. At last, we recall that the choice of a sole dependence of the potential u on the stress σ in ( 2) is not a restriction of the procedure. The introduction of internal variables which are necessary to take into account microstructural evolutions (hardening, texture evolution, etc.) within this scheme is fully manageable. For a complete derivation, the reader is referred to [START_REF] Masson | Self-consistent estimates for the rate-dependent elastoplastic behavior of polycrystalline materials[END_REF].

The mean field approximation

Hereafter, the material is assumed to be composed of N phases (r), occupying volumes Ω r , with a given microstructure described by the characteristic functions χ r (x) equal to 1 if x ∈ Ω r and 0 otherwise. In the following, the microstructure is only considered in a statistical way through the first-order (e.g., volume fractions) and second-order (e.g., covariances) correlation functions.

To derive the homogenised properties, the linearisation procedure defined above would have to be applied at each point x so as to determine the local stress fields σ (x). Obviously, such a problem cannot be addressed without some simplification, since it would require a complete knowledge of the characteristic functions in order to obtain the stress concentration tensors at each point x. Consequently, the usual mean field approximation has to be used. It consists, using the statistical description of the microstructure, in the determination of the phase stress and strain averages, σ r and ε r , e.g., for Hashin-Shtrikmantype estimates, through the use of Eshelby's solution [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] of the inclusion problem. For this step, the affine method classically defines the linearised behaviour (compliance and stress-free strain) with respect to the average fields over each phase. This choice, originally made by [START_REF] Hill | Continuum micro-mechanics of elastoplastic polycrystals[END_REF] for elastoplasticity, implies that the stress and strain phase averages are related by the constitutive equations (1). At any linearisation time ζ , the affine behaviour (5) reduces to

ε r (t) = [q ζ * σ r ](t) + ε 0 r (t, ζ ), q ζ (t) = s + m(ζ )t, m(ζ ) = ∂g ∂σ (σ r (ζ )), (6)
where * denotes the Stieltjes convolution product and introducing the notation . r = • r . It has to be mentioned that refined definitions of the linearised behaviour taking into account, for instance, the intraphase second moment of the fields as proposed by [START_REF] Brenner | A modified affine theory for the overall properties of nonlinear composites[END_REF] could be used in (6). This possibility will be explored in a future work.

To obtain the homogenised properties at time t = ζ , the correspondence principle can be used. It helps to transform the linear thermo-viscoelastic behaviour into a symbolic thermoelastic one for which the classical homogenisation techniques apply. For that goal, the use of the Laplace-Carson transform, namely

f (p) = p +∞ 0 f (t) e -pt dt, ( 7
)
where p is the complex variable, is especially convenient. Applied to (6), it leads to

εr (p) = qζ (p) : σ r (p) + ε0 r (p), qζ (p) = s + 1 p m(ζ ). ( 8
)
The symbolic overall properties are then classically given by

Q ζ (p) = qζ (p) : B r (p, ζ ) , E 0 (p, ζ ) = t B r (p, ζ ) : ε0 r (p, ζ ) , ( 9
)
where • denotes the volume average over Ω, B r is the symbolic stress concentration tensor over phase (r), Q ζ is the symbolic effective compliance and E 0 is the symbolic effective stress-free strain. The B r tensors depend on the homogenisation scheme adopted, according to the structural morphology of the considered material.

Once the overall symbolic properties are obtained, it is necessary to proceed to their Laplace-Carson inversion in order to get the homogenised behaviour in the direct time space. Except for very particular cases (see Section 3), analytical expressions of the different symbolic quantities ( Q(p), B r (p), etc.) are not available. As a consequence, their inverse transforms must be performed numerically. To do that, Masson and Zaoui chose the collocation method [START_REF] Schapery | Approximate methods of transform inversion for viscoelastic stress analysis[END_REF] which approximates the transient part of any viscoelastic function by a finite Dirichlet series. From the values of the symbolic function to invert at M collocation points p in the Laplace-Carson space, the inverse transform is computed using a least square minimisation with respect to the M coefficients of the series. The main advantage of this method lies in the control of the accuracy of the inversion which can be improved with an increasing number of collocation points. For the considered case, it requires in practice about 20 relaxation times equispaced over six decades on a logarithmic time scale [START_REF] Turner | Self-consistent modeling of visco-elastic polycrystals: application to irradiation creep and growth[END_REF]. The counterpart of the method is that it requires, at each time ζ , as many resolution operations of the homogenisation equations (9) as collocation points. Though tractable, this implies CPU time-consuming computations. These observations motivate the study of an alternative approximate treatment in view of faster computations. This treatment is now first defined for linear viscoelasticity (Section 3) and then applied to the affine modelling of nonlinear heterogeneous viscoelasticity (Section 4).

An approximate treatment for linear viscoelasticity

The direct inversion method

The "direct inversion method" is an approximate method of Laplace-Carson (LC) transform inversion initially proposed by [START_REF] Schapery | Approximate methods of transform inversion for viscoelastic stress analysis[END_REF]. It is based on a mathematical property of the LC transform (7) which, using the change of variables w = log 10 (pt), reads alternatively

         f (p) = +∞ -∞ f 10 w p g(w) dw,
g(w) = 10 w e -10 w log e 10.

(10)

The key-point of the direct inversion method lies in the form of the function g(w) which is quasi-null everywhere except for -2 w 1 (Fig. 1) and satisfies It can thus be advantageously replaced, as a first approximation, by a Dirac function δ(ww 0 ). This leads to the following relation

f (p) f (t) t =10 w 0 /p or, inversely, f (t) f (p) p=10 w 0 /t . (12)
Thus, knowing the value of f for a single value of p = 10 w 0 /t is sufficient for an approximate estimation of f (t). To determine the optimal value of w 0 , it is necessary to consider the shape of f . For functions presenting a quasi-linear variation versus log(t), Schapery obtained the value w 0 = -C/ log e 10 with C the Euler constant. The method becomes exact if the considered function is rigorously proportional to log(t). Schapery further proposed to use an approximate value for C leading to 10 w 0 = 0.5. Therefore, the initially proposed direct method provides a reasonable approximation if the function f (t) presents a quasi-linear behaviour over a range of 3 decades on a logarithmic time scale, due to the form of g(w). For polycrystals, according to [START_REF] Turner | Self-consistent modelling of nonlinear visco-elastic polycrystals: an approximate scheme[END_REF], the stress and strain responses as well as creep compliances are smooth functions of time and can therefore be transformed in accordance with the so-defined direct method. The authors further concluded that this approximate treatment is limited to responses that do not exhibit rapid variations in time, such as creep deformation tests.

In this section, restriction is made to linear viscoelasticity in order to reconsider, on simpler situations, the most appropriate value for w 0 for the coming applications which are mainly focussed on the creep response of nonlinear viscoelastic heterogeneous materials. The way other deformation paths (i.e., more complex than creep) may be treated will also be indicated.

Creep response

For sake of simplicity, we consider in the following an isotropic material whose constituents obey an isotropic incompressible Maxwellian behaviour defined by

ε(t) = m e i σ (t) + m v i σ (t). ( 13
)
The scalar constants m e i and m v i are, respectively, the elastic and viscous shear compliances of phase i, ε is isochoric and σ is the stress deviator. Thus, the relaxation spectrum reduces to a single relaxation time τ i = m e i /m v i . The overall behaviour of this heterogeneous material is known to be no more Maxwellian, due to the long range memory effect. Relation ( 13) can be given in the following integral form

ε(t) = m ve i * σ (t) (14) with m ve i (t) = m e i + tm v i .
Using the LC transform technique, relation ( 14) can be transformed into

ε(p) = m ve i (p) σ (p) ( 15
)
with m ve i (p) = m e i + (1/p)m v i or equivalently:

σ (p) = 2 μi (p)ε(p), μi (p) = p 2m e i (p + 1/τ i ) . ( 16
)
The viscoelastic behaviour ( 14) is transformed into a symbolic elastic behaviour (15) in the LC space in which classical linear homogenisation schemes apply. The response is obtained in the LC space and it has then to be inverse transformed. For this last step, we analyse the consequences of the approximate method evoked above. Let us consider a creep loading under homogeneous boundary conditions at t = 0 of the form Σ = S 0 H (t), with H (t) the Heaviside function and Σ the overall stress deviator. We start with the most simple case of a homogeneous material, i.e., consisting in a single phase, so that σ (t) = Σ . The strain response reads, for t > 0

ε(t) = m e + tm v S 0 LC -→ ε(p) = m e + 1 p m v S 0 , ( 17
)
where LC denotes the LC transform. So, ε is linear with respect to t and not to log(t) as considered by Schapery and Turner et al. It follows from (17) that the direct method leads in that case to the exact response only for 10 w 0 = 1 (i.e., p = 1/t in ( 12)) instead of 10 w 0 = 0.5 as proposed initially, which would give a twice too large creep rate. It must also be remarked that the choice 10 w 0 = 1 is the one made implicitly by [START_REF] Turner | Self-consistent modelling of nonlinear visco-elastic polycrystals: an approximate scheme[END_REF] whereas they were supposed to use 10 w 0 = 0.5 according to their reference to [START_REF] Schapery | Approximate methods of transform inversion for viscoelastic stress analysis[END_REF].

The effect of the direct inversion method applied on a heterogeneous material may be more difficult to analyse since, in that case, not only the shape of the local strain response has to be considered, but also the evolution with time of all other quantities such as local stress and localisation tensors. For this reason, we still consider the relatively simple case of a two-phase material. The local strain reads, according to the direct method:

ε(t) = m e i + t 10 w 0 m v i σ (t) with ε(t) = ε(p) p=10 w 0 /t , σ (t) = σ (p) p=10 w 0 /t . ( 18
)
Important is to recognise that the application of the direct method leads to a linear (time-dependent) problem that can be treated directly in the time space with classical homogenisation schemes, and it is therefore very efficient numerically. Note however that some long range memory effect is still considered in this so-called "quasi-elastic" approximation, as shown by the fact that the overall behaviour is not Maxwellian. It can also be remarked that it is equivalent to apply the direct method to the quantities homogenised in the LC space with local behaviour given by ( 15) or to carry out homogenisation directly in the time space with local behaviour defined by ( 18). The error introduced by the direct inversion method can be estimated by comparing the time derivative of ( 18)

ε(t) = m e i + t 10 w 0 m v i σ (t) + 1 10 w 0 m v i σ (t) (19)
with Eq. ( 13). The direct method introduces an extra time-dependent elastic compliance tm v i /10 w 0 and therefore the Maxwellian character of the local behaviour is lost. Since the effect of this extra term vanishes in the steady-state creep regime ( σ = 0), the steady-state strain rate predicted by the direct method is the exact one for 10 w 0 = 1. In the transient regime, the effect of this term can be investigated by calculating the error (F -F DM )/F , where F (t) refers to the creep function (defined as E(t) = F (t)Σ 0 H (t)) obtained by the rigourous LC inversion and F DM the same quantity when applying the direct method.

This error can be expressed easily if one considers the particular microstructure (say, a "parallel" microstructure) which would correspond to a Voigt-type estimation ε = E, with E the overall strain. In that case, the exact 1 and the approximate overall creep functions read respectively:

         F (t) = m e 1 m e 2 c 2 m e 1 + c 1 m e 2 1 + t τ + (τ 2 -τ )( τ -τ 1 ) τ 2 1 -e -t/ τ , F DM (t) = m e 1 m e 2 c 2 m e 1 + c 1 m e 2 1 + 10 -w 0 t τ + (τ 2 -τ )( τ -τ 1 ) τ 2 1 - 1 1 + 10 -w 0 t/ τ (20) with τ = (c 1 m e 2 τ 1 + c 2 m e 1 τ 2 )/(c 2 m e 1 + c 1 m e 2
) and τ = τ 1 τ 2 / τ . Corresponding results are plotted in Fig. 2 for 10 w 0 = 0.5 and 10 w 0 = 1, different ratios τ 2 /τ 1 of the relaxation times of the two phases, identical volume fractions c i of both phases, and k = m e 1 /m e 2 = 1. As discussed above for the homogeneous material, taking 10 w 0 = 1 greatly improves the approximation as compared to the original proposition of Schapery, and this for both transient and steady-state regimes. As shown in Fig. 2, the error increases, for 10 w 0 = 1, with increasing values of τ 2 /τ 1 . Indeed, it is worth noting that when the two phases have identical relaxation times, the ratio m ve 1 / m ve 2 of the symbolic compliances does not depend on p and equals m e 1 /m e 2 . It can then easily Fig. 2. Normalised difference between the exact creep function and the one predicted by the quasi-elastic approximation for 10 w 0 = 0.5, 10 w 0 = 1, and several ratios τ 2 /τ 1 , in the case of a two-phase material with "parallel" Maxwellian phases (c 2 = 0.5, k = 1).

1 In the sequel, the term "exact" supposes that the microstructure realises the adopted linear scheme.

be proved that in that case, the local stresses are constant during the whole creep test, and consequently the direct method gives the exact LC inversion as for a single phase material. The accuracy of the direct method is therefore especially dependent on the heterogeneity of the relaxation times of the different phases.

For the coming applications to polycrystals, a more realistic microstructure has to be considered. As an illustration, the overall behaviour and local mechanical fields are estimated through the classical self-consistent scheme, for which the quasielastic approximation can be compared to the exact solution. According to [START_REF] Rougier | Représentation spectrale en viscoélasticité linéaire des matériaux hétérogènes[END_REF] or [START_REF] Beurthey | Structural morphology and relaxation spectra of viscoelastic heterogeneous materials[END_REF], the overall relaxation shear modulus μ(t) can be derived from its Laplace-Carson transform μ(p) which obeys the second-order equation

X 2 + (2 -5c)k(p + 1/τ 1 ) 3(p + 1/τ 2 ) + 5c -3 3 X - 2k(p + 1/τ 1 ) 3(p + 1/τ 2 ) = 0 (21) with X = μ/ μ1 and c = c 2 . The solution μ(p) reads          12 m e 1 μ(p) = (5c -3)p p + 1/τ 1 + (5c -2)p p + 1/τ 2 + K p [P (p)] 1/2 (p + 1/τ 1 )(p + 1/τ 2 ) , K = (2 -5c) 2 k 2 + 2k(6 -5c)(5c + 1) + (5c -3) 2 1/2 , P (p) = (p + 1/θ 1 )(p + 1/θ 2 ), (22)
where -1/θ 1 and -1/θ 2 are the real negative solutions of the second-order equation

K 2 p 2 + 2 (2 -5c) 2 k 2 /τ 1 + k(6 -5c)(5c + 1)(1/τ 1 + 1/τ 2 ) + (5c -3) 2 /τ 2 p + (2 -5c) 2 (k/τ 1 ) 2 + 2k(6 -5c)(5c + 1)/τ 1 τ 2 + (5c -3) 2 /τ 2 2 = 0. ( 23
)
From ( 22) and ( 23), the overall response to any stress Σ (t) or strain E(t) history can be derived, after inversion of the Laplace-Carson transformation, through the relation Σ (p) = 2 μ(p) E(p) whereas the per phase average variables εi (p) and σ i (p) are obtained from the macroscopic quantities through the average strain Âi (p) or stress B i (p) concentration tensors, which reduce here, due to isotropy, to the scalar functions Âi (p) and B i (p), e.g.,

εi (p) = Âi (p) E(p), σ i (p) = 2 μi (p)ε i (p), Âi = 5 μ 3 μ + 2 μi . ( 24
)
The exact inverse transforms of the symbolic creep functions F (p) and fi (p), defined as

     Σ (t) = S 0 H (t), E(p) = F (p)S 0 , εi (p) = fi (p)S 0 , F (p) = 1 2 μ(p) , fi (p) = Âi (p) 2 μ(p) = 5 2[3 μ(p) + 2 μi (p)] ( 25
)
are given in Appendix A. An illustrative example of the overall and average phase responses for creep is given in Fig. 3 for τ 2 = 10τ 1 , by comparing the exact solutions to their approximations according to the direct inversion method. It is clear that the initial choice of Schapery (10 w 0 = 0.5) does not give satisfactory results, a much better one being obtained when taking 10 w 0 = 1. As shown in Fig. 3(b), the direct method is not only able to predict the overall behaviour with a very good accuracy, but also phase average responses are very well estimated. For the presented case, the strain difference between exact and approximate solutions never exceeds 5%, neither in the transient nor in the steady-state regime. The origin of this very good agreement lies in the almost constant local stresses for t > 0.

More generally, according to the known property of the LC transform for the initial and final values, i.e.,

     lim |p|→+∞ f (p) = lim t →0 + f (t), lim p→0 f (p) = lim t →+∞ f (t), ( 26
)
the direct method yields exact limits for f (t) when t → 0 + and t → +∞ whatever the choice of ω 0 . But for sufficiently smooth functions f (t), the same property stands near t = 0 for the time derivative ḟ (t) only for 10 w 0 = 1, since the asymptotic behaviour of such functions obeys the relation: where O(x) means "of the order of x". A similar conclusion stands for f (t)/t for large values of t . Nevertheless, as commented below, the practical importance of these properties depends on the values of the involved relaxation times with respect to the observation time interval.

   if f (t) = f (0 + ) + ḟ (0 + )t + O t 2 , t > 0, then f (p) = f (0 + ) + ḟ (0 + ) p + O 1 p 2 , Re(p) > 0, (27)

Relaxation

Let us now consider the response of the same material when subjected to a relaxation strain path at t = 0, i.e., to an overall applied strain E(t) = E 0 H (t).

If the material is homogeneous, the stress response reads, for t > 0

σ (t) = E 0 m e e -t/τ LC -→ σ (p) = E 0 m e τp 1 + τp . ( 28
)
It is clear that, in that case, the direct method as defined above does not lead to the exact result whatever the choice of w 0 .

A different value for w 0 may be proposed to minimise the gap between the exact and approximate stress responses, but this value must depend on the relaxation time τ of the material. Furthermore, the direct method provides at large t a much slower stress decrease as compared to the exact solution, whatever w 0 . Note however that, according to (27), only the value 10 w 0 = 1 provides the exact stress rate at t = 0. As for creep, the error introduced by the direct method can be calculated easily for a two-phase material with a "parallel" microstructure. It reads now (R -R DM )/R(0), with R the relaxation function defined by Σ (t) = R(t)E 0 H (t). The exact and approximate overall relaxation functions are given by

       R(t) = c 1 m e 1 e -t/τ 1 + c 2 m e 2 e -t/τ 2 , R DM (t) = c 1 m e 1 1 1 + 10 -w 0 t/τ 1 + c 2 m e 2 1 1 + 10 -w 0 t/τ 2 . ( 29
)
Results for τ 2 = 2τ 1 and τ 2 = 10τ 1 are shown in Fig. 4. As for the single phase material, the agreement is not satisfactory for both 10 w 0 = 0.5 and 10 w 0 = 1, especially at large reduced time t/τ 1 where the stress decrease is really too slow. According to Fig. 5, the same conclusion can be drawn for a disordered microstructure corresponding to the self-consistent scheme for which the exact inverse transforms of local and overall symbolic relaxation functions ri (p) and R(p) are derived in Appendix A. In that case, taking p = 0.5/t (i.e., 10 w 0 = 0.5) in Eq. ( 12) gives better predictions than p = 1/t at large t/τ 1 , but the slope at t/τ 1 = 0 becomes simultaneously too steep. The direct method with a relation between p and t of the form p = α/t , with α a constant scalar, seems therefore to be mainly limited to creep loading. As for creep, a better way to treat relaxation loading with the quasi-elastic approximation is to choose a function p = p(t) for which the direct method gives the exact response for a homogeneous material. Relation ( 12) is thus transformed into

     E(t) = E 0 H (t), Σ (p) = R(p)E 0 , σ i (p) = ri (p)E 0 , R(p) = 2 μ(p), ri (p) = 2 μi (p) Âi (p) = 10 μ(p) μi (p) 3 μ(p) + 2 μi (p) (30)
f (t) f (p) p=h(t )/t , (31)
where the function h(t) is adjusted with respect to the relaxation function of the single phase, Eq. ( 28):

h(t) = t e -t/τ τ (1 -e -t/τ ) . ( 32
)
It must be remarked that h(t) now depends on the relaxation time τ of the considered phase. Thus, for the application to heterogeneous materials with phases exhibiting different τ i , τ in Eq. ( 32) could be replaced by an average quantity, say τ . Several definitions for τ may be adopted. Among a few attempts, the best results were obtained here for τ defined as τ = τ . Results are presented in Figs. 4 and5. Clearly, the proposed treatment greatly improves the predictions at both local and global scales, as compared to the original proposition p = 10 w 0 /t . The relaxation functions quickly decrease in good agreement with the exact solutions, and the overall response is particularly close to the exact one. However, local response are not predicted as precisely as for creep. Indeed, as shown in Fig. 5, the relaxation of phase 1 is too slow whereas that of phase 2 is too fast. This may originate from the choice of a unique τ for the two phases. Finally, it is worth noting that, as for creep, the proposed treatment also leads to the exact solution if τ 1 = τ 2 , whatever the value of k. The error introduced by the quasi-elastic approximation has clearly to be related to the local relaxation time ratio τ 1 /τ 2 .

Other deformation paths

The good agreement obtained for a relaxation loading shows that the direct inversion method is not limited to creep. For example, if the material is homogeneous and subjected to a constant strain rate loading of the form ε(t) = Ė0 H (t), the stress response reads, for t > 0

σ (t) = Ė0 m v 1 -e -t/τ LC -→ σ (p) = Ė0 m v 1 1 + τp . ( 33
)
Thus, the same function h(t) as the one given in (32) can be used. As for relaxation, an average τ must be defined for heterogeneous materials and here again, the exact solution is attained if all phases exhibit the same relaxation time. A similar procedure could also be used for more complex loading paths (cycling loading, for instance).

Conclusion

The application of the direct inversion method to simple isotropic materials with (linear) Maxwellian local behaviour has allowed us to analyse its domain of validity. First of all, it has been remarked that the initial Maxwellian character of the local behaviour is lost when the direct method is applied. In the case of creep loading, this results in an additional time-dependent local elastic compliance. Due to the shape of the single phase creep function, an improved treatment with 10 w 0 = 1 has been proposed, which leads to the exact asymptotic regime for heterogeneous materials at t = 0 and t → ∞. Similarly, we have shown how other deformation paths can be treated by using a relation between p and t of the form p = h(t)/t , where the function h(t) has to be adjusted on the response of a single phase material. In the case of relaxation, an average relaxation time τ must however be defined and for this several choices are possible. It has been shown that the direct method provides the exact local and global responses if the relaxation times of the different phases are identical, and therefore its accuracy is expected to become worst with increasing spread of local relaxation times. A very good agreement at both local and global scales has been found in the presented numerical applications, suggesting that the direct method may be also very efficient for nonlinear problems. Since the inversion of the LC transform becomes immediate, the direct method is very powerful numerically. The counterpart of this is that, contrary to the collocation method, the accuracy of the method cannot be improved by means of improved numerical convergence. It has been noticed that the long range memory effects are still considered, that is the overall response is nowhere assumed to be Maxwellian, in accordance with [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF].

The "quasi-elastic" approximation of the affine formulation

Consequences of the direct method

We now go back to the nonlinear frame and study the application of the direct method in that context. The local constitutive equations have the form (1), i.e., they define a nonlinear Maxwell-type behaviour. According to the results obtained in Section 3, the overall creep responses are aproximated with w 0 = 0, whereas other choices are more appropriate for other loading paths. Whatever this choice, the symbolic homogenisation relations (9) now read

Q ζ (t) = q ζ (t) : B r (t, ζ ) , E 0 (t, ζ ) = t B r (t, ζ ) : ε 0 r (t, ζ ) . ( 34
)
The unknown tensors B r (t, ζ ) are obtained by solving the homogenisation problem in the LC space (first equation of ( 9)) and applying the variables change p = h(t)/t , with h(t) = 1 for a creep loading. Equivalently, as shown above, homogenisation can be carried out directly in the time space. As for the prestrain field ε 0 r (t, ζ ), it does not need this approximate inversion since it is explicitly defined in the time space according to (4) and ( 5). Nevertheless, it is emphasised that this field does not identify with the one that would be obtained with the "classical" affine formulation of [START_REF] Masson | Self-consistent estimates for the rate-dependent elastoplastic behavior of polycrystalline materials[END_REF] since the average stresses per phase are generally different. This observation is closely linked to the notion of linear comparison medium. Indeed, the new definition of the homogenised compliance and prestrain (34) corresponds to an effective behaviour that reads

E(t) = Q ζ (t) : Σ(t) + E 0 (t, ζ ). ( 35
)
It thus refers to a thermoelastic (with time-dependent characteristics), instead of thermoviscoelastic, linear comparison medium. That is why this approach is called "quasi-elastic", the convolution products of the initial problem being replaced by algebraic ones. Nevertheless, as explained above, it retains the essential features of the nonlinear viscoelastic behaviour. Particularly, the elastic and viscous behaviours are not treated separately during the transition scale, e.g., the linearised overall behaviour is not supposed to be of the Maxwell-type. Besides, the definition of the eigenstrains ε 0 r is left unchanged, which keeps its hereditary nature.

The "quasi-elastic" affine procedure

Formulation

Using the above results, the "quasi-elastic" affine procedure can be derived as follows. Let us consider, e.g., a macroscopic stress loading path Σ(t) on the time interval [0, ζ ]. We aim at the determination of the effective response E(ζ ). Due to the hereditary character of the behavior, we must assume that the problem was solved until t = ζδt . The phase stress averages are thus known on [0, ζδt] and the unknowns are σ r (ζ ). For a given linearisation time ζ , the linearised behavior reads, at

any time t ζ ,                    ε r (t) = q ζ (t) : σ r (t) + ε 0 r (t, ζ ), q ζ (t) = s + m(ζ )t, ε 0 r (t, ζ ) = t 0 ε0 r (θ, ζ ) dθ, ε0 r (θ, ζ ) = g σ r (θ) -m(ζ ) : σ r (θ). ( 36
)
Due to the quasi-elastic approximation, it is only required to apply relation (36) for t = ζ . The Hashin-Shtrikman-type homogenisation procedure leads to the following interaction equation

   ε r (ζ ) -E(ζ ) = -Q ζ (ζ ) : σ r (ζ ) -Σ(ζ ) , Q ζ (ζ ) = S E -1 -I -1 : Q * ζ (ζ ) (37)
with Q ζ the interaction tensor, S E the Eshelby tensor and Q * ζ the compliance of the reference medium which depends on the chosen homogenisation scheme. The quasi-elastic affine formulation is entirely defined by ( 34), ( 36) and (37). Thus, its implementation only requires the resolution of a set of nonlinear equations which reads F (σ r (ζ )) = 0. Roughly, it reduces the CPU time by a factor 10 with respect to the resolution with use of the collocation method. The problem is solved with a classical fixed-point iterative algorithm. Note that an additional inner iterative resolution is required to solve (34) for the particular case of the self-consistent scheme.

Comparison with former works

A previous application of the quasi-elastic approximation (hereafter denoted TTW) was made by [START_REF] Turner | Self-consistent modelling of nonlinear visco-elastic polycrystals: an approximate scheme[END_REF] to derive the overall properties of nonlinear viscoelastic polycrystals. The differences between their formulation and our proposition are significant enough to be emphasised. Starting from the good agreement obtained with the quasi-elastic method in the linear case, for the particular case of a creep loading, as already mentioned by [START_REF] Laws | Self-consistent estimates for the viscoelastic creep compliances of composite materials[END_REF]McLaughlin (1978), Yurner et al. (1994) proposed an empirical extension to the nonlinear case. Their proposed interaction equation, namely

εr (t) -Ė(t) = -Q TG (t) : σ r (t) -Σ(t) -n ˙ Q TG (t) : σ r (t) -Σ(t) , ( 38
)
where ˙ Q TG identifies in the viscoplastic regime with the accommodation tensor of the tangent model of [START_REF] Molinari | A self-consistent approach of the large deformation polycrystal viscoplasticity[END_REF], uses the assumption of a power-law viscous flow with a unique stress sensitivity coefficient n. It can be checked that such a law leads to a tangent description of the steady-state regime. This property is the sole justification of the proposition (38) which is, roughly speaking, an empirical interpolation between a purely elastic and a tangent viscoplastic regime. On the contrary, according to our approach, the local compliance in the steady-state limit is expressed as q(t) = s + mt with m independent on t since lim t →+∞ σ r (t) = 0. The behaviour becomes purely viscoplastic (or nonlinear viscous) and the formulation obtained, still referring to a thermoelastic linear medium, identifies with the affine procedure for viscoplastic material which differs significantly from the tangent model [START_REF] Masson | An affine formulation for the prediction of the effective properties of nonlinear composites and polycrystals[END_REF]. Especially, for highly nonlinear behaviours, the tangent model tends to the Reuss bound [START_REF] Lebensohn | A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: application to zirconium alloys[END_REF], which is not the case of the affine procedure. However, according to [START_REF] Molinari | Averaging models for heterogeneous viscoplastic and elastic viscoplastic materials[END_REF], the tangent model may be seen as one particular procedure within a more general approach.

Apart from this remark which only concerns the purely nonlinear viscous regime, a specific difference exists in the transient regime where the viscoelastic coupling takes place. Indeed, the TTW approach does not take into account the influence of the whole stress history, since the behaviour is linearised as follows

ε r (ζ ) = s + m sec (ζ )ζ : σ r (ζ ) (39)
with m sec the viscous linearised secant compliance. For a power-law viscous flow with stress sensitivity n, Eq. ( 39) gives [START_REF] Hutchinson | Bounds and self-consistent estimates for creep of polycrystalline materials[END_REF])

ε r (ζ ) = s + m(ζ )ζ : σ r (ζ ) + ε 0 r (ζ ), ε 0 r (ζ ) = (1 -n)m sec (ζ )ζ : σ r (ζ ) (40)
with the tangent compliance m = nm sec . It is to be stressed that the TTW formulation (40) differs from the affine one presented in this paper, since it relies on an instantaneous definition of the stress-free strain ε 0 r (ζ ) whereas the quasi-elastic affine approach accounts for the hereditary nature of the behaviour. Consequently, the TTW approach neglects the cumulated influence of the whole stress history on the current linearised behaviour.

Application to polycrystals

We now consider the elastoviscoplastic behaviour of various polycrystals. We assume that the material is made of randomly distributed grains with an equiaxed shape. Consequently, all the grains with the same crystallographic orientation define a phase with a given volumic fraction and an isotropic distribution. In the following, we compare the elastoviscoplastic response of this material as predicted by the self-consistent "classical" affine formulation [START_REF] Masson | Self-consistent estimates for the rate-dependent elastoplastic behavior of polycrystalline materials[END_REF] and the self-consistent quasi-elastic affine treatment proposed here.

Constitutive behaviour

We assume that the viscoplastic deformation is only due to dislocation glide on crystallographic slip systems. According to the kinematics of plasticity, it follows:

g = s γs R s . ( 41
)
The Schmid tensor R s , characterizing the orientation of the slip system (s), is given by

R s = 1 2 (n s ⊗ m s + m s ⊗ n s ), (42)
with n s and m s the unit vectors normal to the slip plane and along the slip direction. In the sequel, we consider a classical power-law type behaviour of the form

γs = γ 0 s τ s τ 0 s n , ( 43
)
with γ 0 s the reference slip rate and τ 0 s the reference shear stress of system (s) and n the material stress sensitivity. The resolved shear stress is defined as τ s = R s : σ . For simplicity, the reference shear stresses and Schmid tensors will be kept constant here. In other words, we shall not consider hardening effects and texture evolution.

We consider polycrystals with two types of crystalline structures: face-centered cubic (FCC) and hexagonal close-packed (HCP). For FCC polycrystals, glide occurs on the twelve octahedral slip systems {111} 110 with the same reference slip rate γ 0 s and shear stresses τ 0 s . Such polycrystals consequently obey a weakly anisotropic plastic behaviour. For HCP polycrystals, we assume three different slip modes to be active: prismatic {1 100} 11 20 , basal {0001} 11 20 and first-order pyramidal c + a {1 101} 1 21 3 . All systems present the same reference slip rate but we consider a ratio of 6 on the reference shear stresses between the prismatic system and the other ones. The prismatic system is taken as the easy-glide system. Contrarily to the FCC case, these polycrystals present a highly anisotropic plastic behaviour. They can be thought as zirconium-type polycrystals. For both kinds of polycrystals, we assume a uniform power-law exponent n = 5, which is representative of hightemperature deformation mechanisms. Besides, the elasticity is isotropic with a Poisson ratio of 1/3. For the FCC case, we have τ 0 /E y = 10 -3 and for the HCP case, we have 2τ 0

pr /E y = 10 -3 (where the index pr stands for "prismatic") with E y the Young modulus. Finally, the polycrystals under consideration here are untextured.

Creep responses

We performed computations for tensile creep tests with instantaneous loading at t = 0. Fig. 6 shows the evolution of the overall equivalent creep strain E eq , normalised by the equivalent elastic strain E e eq , for FCC and HCP polycrystals. It can be observed that the affine effective response predicted with the use of the direct method with p = 1/t lies very close to the one obtained with the collocation method in the transient and steady-state regimes. Thus, our approach is shown to be relevant even in the case of a highly anisotropic viscoplastic behaviour (i.e., with a high contrast on the relaxation times). To help understanding the accuracy of the "quasi-elastic" approximation, the results of a decoupled model identifying separately the elastic strains and the viscous strains are also reported in Fig. 6. This approach assumes that the overall behaviour is of the Maxwell-type. Consequently, in the case of a creep test with an instantaneous elastic loading at t = 0, it is unable to describe the transient regime and leads directly for t > 0 to the steady-state viscous response. Comparing the FCC and HCP cases, it is observed that the discrepancy between the "full" affine scheme and this crude approximation increases with the viscoplastic anisotropy whereas the "quasi-elastic" formulation remains valid in both cases.

For the HCP polycrystal, it is particularly interesting to study the predictions of the phase field averages, since the interphase fields heterogeneity is very pronounced. Due to the transverse isotropy of the phase responses around the loading direction, computations can be restricted to one spherical standard triangle. The stress and strain responses are reported in Fig. 7 for four orientations in the standard triangle. The evolution of the stresses shows that the transient response of the "soft" and "hard" orientations is correctly captured by our approach and that the relaxation times are very similar in the two approaches. It can be stressed that the decoupled model would predict, as for the effective behaviour, a constant strain rate (versus time) in each phase for t > 0 and therefore also a constant stress. Apart from this remark, it is seen that the large interphase strain heterogeneity, which is physically relevant with respect to the grain-to-grain incompatibility, is satisfactorily restituted with the present proposition. It is recalled that, in the viscoplastic regime, this formulation coincides with the affine formulation. From that point of view, our approximate treatment is expected to be more accurate than former approximate approaches which are known to lead in the steady-state regime either to the Taylor bound [START_REF] Weng | Self-consistent determination of time dependent behavior of metals[END_REF][START_REF] Weng | A self-consistent relation for time-dependent creep of polycrystals[END_REF], as shown by [START_REF] Masson | Self-consistent estimates for the rate-dependent elastoplastic behavior of polycrystalline materials[END_REF], or the Reuss bound [START_REF] Turner | Self-consistent modelling of nonlinear visco-elastic polycrystals: an approximate scheme[END_REF] for a highly nonlinear viscoplastic potential. The good agreement of our treatment with the affine formulation has important practical implications. Indeed, since the computation time is divided by a factor of about 10, it allows an identification of the model on a wide experimental database and makes it a powerful tool for the study of polycrystal plasticity (see [START_REF] Brenner | Influence de la microstructure sur le comportement en fluage thermique d'alliages de zirconium : analyse experimentale et mise en oeuvre de méthodes d'homogénéisation[END_REF][START_REF] Letouzé | Residual strain distribution in Zircaloy-4 measured by neutron diffraction and estimated by homogenization techniques[END_REF][START_REF] Brenner | Thermal creep of Zr-Nb1%-O alloys: experimental analysis and micromechanical modelling[END_REF].

Relaxation responses

As shown previously, the quasi-elastic approach has not to be limited to creep. As a first illustration, the response of the FCC polycrystal to a uniaxial relaxation loading has been studied. According to the results obtained in Section 3, we apply the direct method with p = h(t)/t and h(t) defined by (32). For that goal, it is necessary to determine a reference relaxation time for the polycrystal. This step is not straightforward due to the local anisotropy of the material. Among other possibilities, we defined at each linearisation time ζ the "relaxation time" of each phase τ r and the reference relaxation time τ by

τ r = s r ij ij m r ij ij , τ = τ r . ( 44
)
Fig. 8 shows the overall stress responses according to the direct method with p = h(t)/t and p = t and to the classical affine formulation. It can be seen, as expected, that the choice of p = h(t)/t leads to a better agreement in the transient regime, especially at short reduced times, than the choice p = 1/t . Though our definition of the reference relaxation time is very simple and certainly not optimal, this result shows a simple but efficient way to manage with another type of loading than creep with an appropriate definition of the function h(t). Extension to constant strain rate behaviour, for instance, is straightforward (see Section 3.4). To the best of our knowledge, the proposed approach is the first to extend the quasi-elastic treatment to other loadings than creep.

In view of these results, it can be appreciated that our approach reconciles an affine formulation which is known to be one of the most accurate at the time for rate-dependent plasticity [START_REF] Masson | Self-consistent estimates for the rate-dependent elastoplastic behavior of polycrystalline materials[END_REF] with an efficient numerical resolution which renders the computation as simple as for any other nonlinear homogenisation scheme since it uses a thermoelastic linear comparison medium. Additionally, it is recalled that this formulation does not imply any restriction on the chosen linear homogenisation scheme (and so on the form of the interaction equation) or on the form of the nonlinear Maxwell-type behaviour of the constituents: it can be used with more refined forms of the local viscoplastic potential than the one used in this section.

Conclusion

Several major difficulties have to be overcome for the application of homogenisation techniques to nonlinear viscoelastic heterogeneous materials. First of all, a linearisation procedure must be defined. In the present study, the proposition of [START_REF] Rougier | Self-consistent modelling of elastic-viscoplastic polycrystals[END_REF] has been adopted. Next, per phase uniform linear viscoelastic compliances have to be calculated, and for this a reference stress must be defined. The classical (and most simple) choice defining the linearised compliances with respect to the stress phase average has been chosen here. Finally, the mechanical fields obtained in the Laplace-Carson space at each linearisation time have to be transformed back into the physical time space. This last point is more "technical", but it has a great practical importance because it renders the homogeneisation scheme quite CPU time-consuming.

Although an accurate inversion can be carried out numerically, the application to anisotropic polycrystals, with a necessarily associated step of coefficients identification, is hardly workable. For this, we have reconsidered the approximate direct inversion method initially proposed by [START_REF] Schapery | Approximate methods of transform inversion for viscoelastic stress analysis[END_REF] which leads to a reduction of CPU time by one order of magnitude. Its accuracy has been first investigated on a linear heterogeneous material with maxwellian phases. It has been shown that the method is exact if the phases exhibit identical relaxation times, and the improved treatment proposed here renders the direct method applicable to any loading path. Applied to nonlinear polycrystals in the framework of the "classical affine" linearisation procedure, the direct method is shown to describe the exact asymptotic behaviour at t = 0 and t → ∞. The initial thermoviscoelastic linear comparison material is transformed into a thermoelastic one, so that homogenisation can be completely carried out in the physical time space. Moreover, the long range memory effects associated to the viscoelastic coupling are still taken into account. Illustrative examples for FCC and HCP nonlinear polycrystals show that the direct method predicts almost the same mechanical responses in creep and relaxation as the collocation treatment at both local and global scales, and therefore reconciles the accuracy of the affine formulation with a very efficient numerical resolution.

It is emphasized that the "direct method" can also be used in the context of the generalised affine linearisation proposed by [START_REF] Pouya | Linéarisation et homogénéisation en viscoélasticité[END_REF] and/or when the reference stress accounts for intraphase heterogeneity. These aspects are currently under development.

Appendix. Analytical derivation of the creep and relaxation behaviour of the two-phase material of Section 3

The derivation of the original time functions of the Laplace-Carson tranformed functions F (p), fi (p), R(p) and ri (p) defined in ( 25) and ( 30), where μi (p) and μ(p) are given by ( 16) and ( 22) respectively (with the convention τ 1 < τ 2 ), can make use of the Mellin-Bromwich inversion formula

f (t) = 1 2iπ f (p) p
e pt dp, (A.1)

where the straight line , parallel to the y-axis, must leave all the poles and critical points of f (p)/p at left. This line may be complemented on its left by a half-circle centered on it at y = 0 with a radius growing to infinity; the contribution of this half-circle to the integral in (A.1) vanishes according to Jordan's lemma. Due to the term [P (p)] 1/2 in ( 22), which needs to be uniformly defined, a cut must be defined on the real negative axis between -1/θ 1 and -1/θ 2 , with -1/τ 1 < -1/θ 1 < -1/θ 2 < -1/τ 2 < 0. The computation then reduces to that of some residues on the x-axis outside the cut [-1/θ 1 , -1/θ 2 ] and of an integral on this interval. The searched functions read as follows (when k = 1 and k = τ 1 /τ 2 ):

• Creep f 1 (t) = m e 1 τ 1 τ 2 4(1 -c)(kτ 2 -τ 1 ) k(2 -5c) τ 1 + 5c -7 τ 2 - 5k(τ 2 -τ 1 ) (kτ 2 -τ 1 )τ 1 + K τ 1 (θ 1 θ 2 ) 1/2 θ 1 + θ 2 2 + τ 1 (τ 2 -τ 1 ) kτ 2 -τ 1 + k(2 -5c) τ 1 + 5c -7 τ 2 + K (θ 1 θ 2 ) 1/2 t τ 1 - K π J 1 τ 1 + (τ 2 -τ 1 )(J 2 -J 3 ) kτ 2 -τ 1 , b = kτ 2 -τ 1 (k -1)τ 1 τ 2 , J 1 = -1/θ 2 -1/θ 1 [-P (x)] 1/2
x 2 e tx dx,

J 2 = -1/θ 2 -1/θ 1 [-P (x)] 1/2
x e tx dx,

J 3 = -1/θ 2 -1/θ 1 [-P (x)] 1/2
x + b e tx dx, (A.2)

f 2 (t) = m e 1 τ 1 τ 2 4ck(τ 1 -kτ 2 ) - k(2 + 5c) τ 1 + 5c -3 τ 2 - 5k(τ 2 -τ 1 ) (kτ 2 -τ 1 )τ 2 + K τ 2 (θ 1 θ 2 ) 1/2 θ 1 + θ 2 2 + kτ 2 (τ 2 -τ 1 ) kτ 2 -τ 1 + - k(2 + 5c) τ 1 + 5c -3 τ 2 + K (θ 1 θ 2 ) 1/2 t τ 2 - K π J 1 τ 2 + k(τ 2 -τ 1 )(J 2 -J 3 ) kτ 2 -τ 1 , (A.3) F (t) = (1 -c)f 1 (t) + cf 2 (t) = m e 1 4k k(2 -5c) + 5c -3 + k(2 -5c) τ 1 + 5c -3 τ 2 t + K (θ 1 θ 2 ) 1/2 θ 1 + θ 2 2 + t - KJ 1 π . (A.4) • Relaxation r 1 (t) = 1 3m e 1 (1 -c)
3 -5c e -t/τ 1 + K(J 4 -J 3 )τ 1 τ 2 2π(τ 2τ 1 ) ,

x = |x| + x 2 , J 4 = -1/θ 2 -1/θ 1 [-P (x)] 1/2
x + 1/τ 1 e tx dx, (A.5)

r 2 (t) = 1 3m e 1 c
5c -2 k e -t/τ 2 -K(J 5 -J 3 )τ 1 τ 2 2π(τ 2τ 1 ) , J 5 = -1/θ 2 -1/θ 1 [-P (x)] 1/2 x + 1/τ 2 e tx dx, (A.6)

Fig. 1 .

 1 Fig. 1. Shape of the function g(w) defined by Eq. (10).

Fig. 3 .

 3 Fig. 3. Local (f ) and overall (F ) creep functions predicted by the quasi-elastic approximation (dashed lines) as compared to the exact results (continuous lines) for a two-phase material according to the self-consistent scheme, for τ 2 /τ 1 = 10, c 2 = 0.5, k = 1. (a) 10 w 0 = 0.5 and (b) 10 w 0 = 1.

Fig. 4 .

 4 Fig. 4. Normalised difference between the exact relaxation function and the one predicted by the quasi-elastic approximation for two ratios τ 2 /τ 1 in the different cases p = 0.5/t, p = 1/t, and p = h(t)/t. Two-phase material with "parallel" Maxwellian phases, c 2 = 0.5, and k = 1.

Fig. 5 .

 5 Fig. 5. Local (r) and overall (R) relaxation functions predicted by the quasi-elastic approximation (dashed lines) as compared to the exact results (continuous lines) for a two-phase material according to the self-consistent scheme and τ 2 /τ 1 = 10, c 2 = 0.5, k = 1, in the cases (a) p = 0.5/t, (b) p = 1/t, and (c) p = h(t)/t.

Fig. 6 .

 6 Fig. 6. Normalised overall equivalent strain for a tensile creep test for (a) FCC (Σ 33 = 4τ 0 ) and (b) HCP (Σ 33 = 2τ 0 pr ) polycrystals as predicted by the affine model (continuous line), its quasi-elastic approximation (dashed lines) with p = 1/t and a decoupled (Maxwell-type) model (dotted lines).

Fig. 7 .

 7 Fig. 7. Intraphase average stress and strain responses for four different orientations during a tensile creep test (Σ 33 = 2τ 0 pr ) for an HCP polycrystal as predicted by the affine model (continuous line) and its quasi-elastic approximation (dashed lines) with p = 1/t.

Fig. 8 .

 8 Fig.8. Normalised overall equivalent stress during a tensile relaxation test (E 33 = 4τ 0 /E y ) for a FCC polycrystal as predicted by the affine model (continuous line) and its quasi-elastic approximation (dashed lines) with p = 1/t et p = h(t)/t.

Acknowledgement

One of us (R.B.) gratefully acknowledges the Commissariat à l'Energie Atomique for the partial funding of this work.

Note that the integrals J 1 to J 5 may be expressed as functions of Φ 1 , degenerate hypergeometric series in two variables [START_REF] Gradshteyn | Table of Integrals, Series and Products[END_REF]. They have been evaluated numerically for the illustrations reported in Section 3.

When k = 1, Eqs. (A.2) to (A.7) can still be used with K = 5 and J 3 = 0.


Stéphan Beurthey

André Zaoui

Structural morphology and relaxation spectra of viscoelastic heterogeneous materials

We analytically derive the relaxation spectra of a two-phase isotropic material whose phases are isotropic Maxwell media, according to the classical and to the generalized self-consistent schemes. Whereas these spectra are continuous in both cases, they exhibit strong differences which can be associated with the different underlying morphology, either symmetrical (polycrystal-type) in the first case or asymmetrical (composite-type) in the second case. The treatment is extended to the (N + 1)-phase model which allows us to deal with coated inclusions or with an interphase between the matrix and the inclusions: the interphase is shown to strongly modify the resultant spectrum. More general cases are then considered for different kinds of constitutive behaviour as well as for coated fibre reinforced composites. As a whole, the spectral analysis method appears to be an efficient tool for the investigation of the connection between structural morphology and the overall behaviour of viscoelastic heterogeneous materials. spectral analysis / structural morphology / heterogeneous materials

Introduction

The influence of internal morphological characteristics on the overall behaviour of micro-inhomogeneous materials is a matter of increasing work and interest. Up to now, this topic has been mainly developed in the framework of linear elasticity (since [START_REF] Kröner | Bounds for effective elastic moduli of disordered materials[END_REF]. Basic questions such as the way to take into account the inclusion/matrix morphology of composite materials or the disordered constitution of polycrystals have been given satisfying practical answers in this case, especially through the classical self-consistent scheme (C.S.C.S., [START_REF] Hershey | The elasticity of an isotropic aggregate of anisotropic cubic crystals[END_REF] on the one hand, for intricate morphologies where no phase plays any specific morphological role, and through the composite spheres (or cylinders) assemblage [START_REF] Hashin | The elastic moduli of heterogeneous materials[END_REF] and the associated generalized self-consistent scheme (G.S.C.S. or three-phase model, [START_REF] Christensen | Solutions for effective shear properties in three-phase sphere and cylinder models[END_REF] on the other hand, when one phase is disposed as a well-connected matrix and the other ones are dispersed in it. The latter case has been generalized in different ways, either by taking multilayered inclusions into account through so-called (N + 1)-phase models in order to deal with coated particles or fibres [START_REF] Hervé | N-layered inclusion-based micromechanical modelling[END_REF][START_REF] Hervé | Elastic behaviour of multiply coated fibre-reinforced composites[END_REF] or by considering arbitrary composite patterns as basic units of heterogeneous materials through 'morphological pattern-based' bounding or estimating treatments [START_REF] Bornert | Morphologically representative pattern-based bounding in elasticity[END_REF].

Nevertheless the generalization of such approaches to nonlinear behaviour has to overcome strong difficulties which originate in the fact that the linearized (either secant or tangent) moduli of the phases cannot be treated as uniform per phase [START_REF] Zaoui | Structural morphology and constitutive behaviour of microheterogeneous materials[END_REF]. Whereas many current works are devoted to rigorous or approximate treatments of this problem, this paper proposes to explore another route within the same general field in order to enlarge the knowledge of the connection between micro-morphology and the overall behaviour of heterogeneous materials. This approach is based on the extension of known treatments for the case of linear elasticity to linear (non ageing) viscoelasticity. The difficult coupling of elasticity and viscosity, which makes straightforward extensions of Green's techniques inefficient due to the simultaneous occurrence of derivatives of different orders of the mechanical variables in the considered constitutive equations, can then be easily treated with the help of the Laplace transformation technique and the associated 'correspondence principle' (Mandel, 1966, Laws and[START_REF] Laws | Self-consistent estimates for the viscoelastic creep compliance of composite materials[END_REF]: this principle allows one to transform an inhomogeneous viscoelastic problem into a symbolical elastic one which may have a known solution and to go back to real variables by inversion of the Laplace transformation. Such a method makes it possible to study actual viscoelastic heterogeneous materials, such as reinforced polymer matrices or polymer blends, to investigate the mechanical influence of interphases between (possibly coated) inclusions and the matrix or to serve as a basis for further extensions to nonlinear viscoelasticy or rate-dependent elastoplasticity. An additional interest can be found in the use of the spectral analysis of the overall behaviour through the derivation of the relaxation spectra of such materials which is expected to give access to significant information on the micro-morphology. This point can be illustrated as follows: we consider a multiphase material consisting of Maxwellian constituents with local constitutive equations of the form

ε = a : σ + b : σ (1)
with ε and σ the strain and strain tensors and a and b material constant fourth order tensors. Whereas the overall behaviour is known to be linear viscoelastic, it is no more Maxwellian due to the coupling between elasticity and viscosity and the associated 'long range memory effect'. The overall constitutive equations are expected in the form [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF]:

Ė = A eff : + B eff : ˙ + t 0 J(t -τ ) : ˙ (τ ) dτ, (2)
where E and are the macroscopic strain and stress tensors, A eff and B eff are the macroscopic (or effective) analogues of a and b and J(tτ ) is an interaction kernel which reflects the complexity of the delayed mechanical interactions between the phases.

How can this kernel be represented by relaxation spectra? How does it depend on the considered morphology? What do the associated relaxation spectra look like according to this morphology? How can they be attached to a polycrystal-type or to a composite-type morphology? How does an interphase affect the shape of the relaxation spectra? What happens when coated particles or fibres are considered? These are the main questions which are addressed in the following, through the use of the C.S.C.S., of the G.S.C.S. and of the (N + 1)-phase model for simple linear viscoelastic constituents. This is performed through the use of the Laplace(-Carson) transformation technique in order to convert the analysis of the coupling of elasticity and viscosity into the resolution of a symbolical linear elastic problem.

In the first part (Section 2), we determine the shape of the shear relaxation spectra when each phase has only one relaxation time, as a function of the morphology, by comparing the predictions of the 'classical' and the 'generalized' self-consistent models. We then consider (Section 3) the '(N + 1)-phase model' which is based on the configuration of a multilayered spherical inclusion surrounded by the equivalent homogeneous medium. Each phase may have several relaxation times. Thus the existence of an interphase located between the inclusions and a continuous matrix can be taken into account. The comparison of the relaxation spectra indicates an impressive 'signature' of the interface. Other applications are suggested. Finally, we report similar results for cylindrical geometry (Section 4).

Shear relaxation spectra of two-phase isotropic Maxwellian materials

Assuming local isotropy and incompressibility for the sake of simplicity, we write the constitutive equations for each phase in the form:

ė = a i s + b i ṡ, i = 1, 2, (3)
where a i and b i are scalar constants and s and e are the stress and strain deviators. After Laplace-Carson transformation, defined by

f * (p) = p ∞ 0 f (t) e -pt dt, (4)
we get:

s * (p) = 2µ * i (p) e * (p), µ * i (p) = p 2b i (p + 1/t i ) , ( 5
)
with p the complex variable, t i = b i /a i the (single) relaxation time of phase (i) with the intensity 1/2b i . The two-phase material is supposed to be isotropic too. Due to incompressibility, it is completely defined by the overall shear relaxation modulus µ eff (t). According to the theory of spectral analysis, we are looking for this modulus in the form

µ eff (t) = ∞ 0 g eff (τ ) e -t/τ dτ, (6)
where g eff (t) gives access to the effective shear relaxation spectrum -which reduces for each phase (i) to one single line at t = t i with the intensity 1/2b i . Two self-consistent models will be considered concurrently, the classical and the generalized ones. Owing to the correspondence principle, µ eff * (p) is obtained from µ * 1 (p) and µ * 2 (p) through the same equations as those which relate µ eff , µ 1 and µ 2 in the elastic case. These classical quadratic equations have the form:

A.X 2 + B.X + C = 0, (7)
where X = µ eff /µ 1 . For the C.S.C.S. we have:

A = 1, B= 2 -5c 3
µ 2 µ 1 + 5c -3 3 , C= - 2 3 µ 2 µ 1 (8)
with c the volume fraction of phase (2). For the G.S.C.S., with phase (1) as the continuous matrix, A, B and C are second order polynoms in (µ 2 /µ 1 ) and tenth order polynoms in c 1/3 (see Appendix 1). For the sake of simplicity, we define A such that it is non-negative for any c. In both cases, equation ( 7) has only one real positive root.

For the present viscoelastic case, we define X * (p) as:

X * (p) = µ eff * (p) µ * 1 (p) . ( 9
)
Now X * (p) is a solution of an equation of the form [START_REF] Rougier | Représentation spectrale en viscoélasticité linéaire des matériaux hétérogènes[END_REF]:

A * X * 2 + B * X * + C * = 0, ( 10
)
where A * , B * and C * depend on µ * 2 (p)/µ * 1 (p) and c in the same way as A, B and C depend on µ 2 /µ 1 and c. As a result, µ eff * (p) can be expressed as the sum of two terms:

µ eff * (p) = f * 1 (p) ± f * 2 (p), f * 1 (p) = - pB * 4b 1 (p + 1/t 1 )A * , f * 2 (p) = p(B * 2 -4A * C * ) 1/2 4b 1 (p + 1/t 1 )A * . ( 11
)
The first term f * 1 (p) is a rational fraction in p, so that f 1 (t) is a sum of decreasing exponential functions: two for the C.S.C.S., contributing to the spectrum by two single lines at times t 1 and t 2 , but three for the G.S.C.S., contributing to the spectrum by three single lines, one at the matrix relaxation time t 1 and two other ones at times θ 1 and θ 2 which lie somewhere between t 1 and t 2 . Note for the moment that their respective intensities may be negative according to the value of c.

After some reduction using ( 9) and ( 11), f * 2 (p) can be written for the C.S.C.S. as:

f C * 2 (p) = k C p[P C (p)] 1/2 (p + 1/t 1 )(p + 1/t 2 ) , P C (p) = (p + 1/τ 1 )(p + 1/τ 2 ), k C = (2 -5c) 2 b 2 1 + 2(6 -5c)(5c + 1)b 1 b 2 + (5c -3) 2 b 2 2 12b 1 b 2 , ( 12
)
where τ 1 and τ 2 are positive times lying between t 1 and t 2 . For the G.S.C.S., we get:

f G * 2 (p) = k G p[P G (p)] 1/2 (p + 1/t 1 )(p + 1/θ 1 )(p + 1/θ 2 ) , P G (p) = p + 1/τ 1 p + 1/τ 2 p + 1/τ 3 p + 1/τ 4 , ( 13
)
where τ 1 , τ 2 , τ 3 and τ 4 are positive times such that

τ 1 , τ 2 ∈ [t 1 , θ 1 ] and τ 3 , τ 4 ∈ [θ 1 , θ 2 ] and k G is a constant given in Appendix 1.
We can easily derive f 2 (t) from ( 12) or ( 13) by use of the inverse Laplace-Carson transformation which is defined in the complex plane by:

f 2 (t) = 1 2iπ f * 2 (p) p e pt dp, (14)
where , parallel to the y-axis, must have all the poles and critical points of f * 2 (p) on the left; it can be taken here as the y-axis itself since all these poles are real negatives.

The theorem of residues can be used after complementing with a half-circle centered at the origin and located in the left half-plane Re(p) 0, with a radius R growing to infinity; in addition, in order to get a uniform representation of f * 2 (p), one has to define adequate cuts on the real negative axis on [-1/τ 1 , -1/τ 2 ] for the C.S.C.S. and on [-1/τ 1 ; -1/τ 2 ] and [-1/τ 3 ; -1/τ 4 ], for the G.S.C.S. (figure 1). Suffice it now to make use of Jordan's lemma in order to get rid of the integral on the half-circle when R → ∞ and we are left with f 2 (t) as a sum of some residues and integrals on the cuts: the residues are associated with decreasing exponentials, and then with single lines of the spectra, whereas the integrals on the cuts will be responsible for continuous parts of the spectra.

More precisely, we get [START_REF] Rougier | Représentation spectrale en viscoélasticité linéaire des matériaux hétérogènes[END_REF], with t 1 < t 2 : for the C.S.C.S.:

f C 2 (t) = ± k C π - -1/τ 2 -1/τ 1 -P C (x) (x + 1/t 1 )(x + 1/t 2 ) e tx dx -π P C (-1/t 1 ) 1/t 2 -1/t 1 e -t/t 1 + π P C (-1/t 2 ) 1/t 1 -1/t 2 e -t/t 2 , ( 15
)
for the G.S.C.S.:

f G 2 (t) = ± k G π -1/τ 2 -1/τ 1 -P G (x) (x + 1/t 1 )(x + 1/θ 1 )(x + 1/θ 2 ) e tx dx - -1/τ 4 -1/τ 3 -P G (x) (x + 1/t 1 )(x + 1/θ 1 )(x + 1/θ 2 ) e tx dx + π P G (-1/t 1 ) (1/θ 1 -1/t 1 )(1/θ 2 -1/t 1 ) e -t/t 1 -π P G (-1/θ 1 ) (1/t 1 -1/θ 1 )(1/θ 2 -1/θ 1 ) e -t/θ 1 + π P G (-1/θ 2 ) (1/t 1 -1/θ 2 )(1/θ 2 -1/θ 1 ) e -t/θ 2 . ( 16
)
With help of the change of variable τ = -1/x, (15) and ( 16) can be put in the standard form (6) and yield the expected spectra after addition to f 1 (t) and determination of the front sign in (15) and ( 16) according to the condition µ eff (0) = lim |p|→∞ (µ eff (p)) > 0. The resulting spectra consist of a continuous part and discrete lines whose intensity can be checked to be always positive or null. For the C.S.C.S. we get one continuous spectrum defined by the function

g C (τ ) = k C t 1 t 2 π √ τ 1 τ 2 √ (τ -τ 1 )(τ 2 -τ ) τ (τ -t 1 )(t 2 -τ ) , τ ∈ [τ 1 , τ 2 ], (17)
and one or two discrete lines located at τ = t 1 with the intensity (3 -5c)/6b 1 for c 3/5a n da t τ = t 2 with the intensity (5c -2)/6b 2 for c 2/5. For the G.S.C.S., the continuous spectrum is divided into two parts defined by

g G (τ ) = (τ )k G t 1 θ 1 θ 2 π τ 1 τ 2 τ 3 τ 4 (τ -τ 1 )(τ -τ 2 )(τ -τ 3 )(τ 4 -τ ) τ (τ -t 1 )(θ 1 -τ )(θ 2 -τ ) , τ ∈ [τ 1 , τ 2 ] ∪ [τ 3 , τ 4 ], (τ) = 1 if τ ∈ [τ 1 , τ 2 ], (τ) = -1 if τ ∈ [τ 3 , τ 4 ]. (18)
Additional single lines may exist or not at t 1 , θ 1 and θ 2 depending on the volume fraction c and on the parameters

t 1 , t 2 , b 1 and b 2 .
Illustrative examples are reported in figure 2a and 2b for the C.S.C.S. and the G.S.C.S. respectively, for identical values of the material constants c, t 1 , t 2 , b 1 and b 2 so as to stress the influence of the sole morphology.

Other examples would lead to similar conclusions, namely:

-for both models, the resulting spectra are definitely different from Maxwell-type spectra, which would reduce to discrete lines: the continuous part reflects mainly the 'long-range memory effect' and the complicated underlying viscoelastic interactions between the phase domains as expressed, according to the self-consistent scheme, by the (simple or composite) inclusion/matrix interactions; -in addition, the material structural morphology is shown to have a prominent influence on the shape of the continuous spectrum: whereas the symmetric morphology associated with the C.S.C.S. leads to an unbroken continuous spectrum, the composite-type morphology which is accounted for by the G.S.C.S. makes this continuous spectrum split into two separate parts and enhances the mechanical role of the matrix phase. So, these qualitative differences of the spectrum shapes, which are much more visible than the differences of the relaxation functions or of the Cole-Cole diagrams [START_REF] Rougier | Étude du comportement sous irradiation: modélisation micromécanique de l'élastoviscoplasticité[END_REF], can be suggested to be used as 'morphological discriminators' when analysing experimental data.

Generalized Maxwellian bodies: (N + 1)-phase self-consistent modelling

The foregoing analysis may be extended in two (possibly combined) directions: first, more general linear constitutive equations may be considered so as to deal with generalized Maxwellian constituents which are still isotropic but can be compressible and exhibit a finite number of relaxation times; second, the (N + 1)-phase self-consistent model may be used by replacing the two-phase composite sphere of the G.S.C.S. by a composite sphere with an arbitrary number of concentric shells with different moduli. Such a model has already been defined in the case of linear elasticity [START_REF] Hervé | N-layered inclusion-based micromechanical modelling[END_REF]) so that its extension to linear viscoelasticity is straightforward: it is proposed here in the context of spectral analysis in order to emphasize, as it has been done above, the influence of the material structural morphology on the overall response.

The corresponding morphology is related to the dispersion of multiply coated inclusions in a continuous matrix, including the important case of interphases located between the matrix and the inclusions, since the difficult study of the characteristics and role of such interphases is a matter of intensive research in the field of polymer science. In addition, the (N + 1)-phase self-consistent scheme can be taken as a simplified model when studying nonlinear composites by approximating the non uniformity of the phase secant or tangent moduli by a discretized radial variation [START_REF] Bornert | Self-consistent approaches and strain heterogeneities in two-phase elastoplastic materials[END_REF][START_REF] Beurthey | Modélisations du comportement d'alliages de polymères[END_REF].

The following treatment is restricted to linear viscoelasticity: it is concerned with the spectral analysis of the (N + 1)-phase self-consistent model consisting of an N -layered spherical inclusion embedded in an infinite matrix made of the searched homogenized material. Each phase has an isotropic (generalized)-Maxwellian behaviour with several relaxation times. We are looking for the overall relaxation spectra, both for the shear and the bulk moduli, with special attention paid to the case of a two-phase material with a third phase (the 'interphase') located between the inclusions and the matrix. The mathematical treatment is a straightforward extension of the preceding one, with some additional help of numerical computations [START_REF] Beurthey | Modélisations du comportement d'alliages de polymères[END_REF].

Shear relaxation spectra

For the sake of simplicity, we restrict ourselves to phases (i) with a constant Poisson ratio ν i (i = 1, N). Their Laplace-Carson-transformed shear modulus µ * i (p) reads:

µ * i (p) = p 2b 1 i (p + 1/t 1 i ) + • • • + p 2b n i i (p + 1/t n i i ) = pN * i (p) (p + 1/t 1 i ) • • • (p + 1/t n i i ) , ( 19
)
where N * i (p) is a polynom of degree (n i -1) in p. Through the use of the correspondence principle, the overall shear modulus µ eff(N) * (p) is calculated from the known quadratic equation of the elastic problem [START_REF] Hervé | N-layered inclusion-based micromechanical modelling[END_REF] after replacing real moduli by transformed complex ones (see Appendix 2):

A * N (p)X * 2 + B * N (p)X * + C * N (p) = 0, X * = µ * eff(N) (p)/µ * N (p), ( 20
)
where A * N (p), B * N (p) and C * N (p) depend on µ * i (p), ν i and on the external radii R i of the shells (i) located between R i-1 and R i (figure A1). Note that some relaxation times may be common to two adjacent shells (i) and (i + 1): let j i,i+1 be the number of such values. The polynoms A * N (p), B * N (p) and C * N (p) are then of degree d N in p, with d N given by:

d N = 2 N-1 i=1 (n i + n i+1 -j i,i+1 -1). ( 21
)
Except for these differences, the treatment is quite similar to the preceding one. Instead of (11), we have to set:

µ eff(N) * (p) = f (N) * 1 (p) ± f (N) * 2 (p), f (N) * 1 (p) = - B * N (p)µ * N (p) 2A * N (p) , f (N) * 2 (p) = [B * N (p) 2 -4A * N (p)C * N (p)] 1/2 µ * N (p) 2A * N (p) . ( 22
)
In all our numerical calculations, the polynom A * N (p) happened to have d N real negative roots (-1/θ k ) with θ k ∈ [t min , t max ], where t min and t max are the smallest and the largest relaxation times of the constituent phases, so that the rational fraction f (N) * 1 (p) can be decomposed into (d N + n N ) simple elements: thus the original f (N) 1 (t) is a sum of (d N + n N ) decreasing exponential functions with either positive or negative amplitudes, depending on the material parameters and phase volume fractions. As for f (N) * 2 (p), we also found numerically that the polynom [B * N (p) 2 -4A * N (p)C * N (p)] has 2d N real negative roots (-1/τ k ) with τ k ∈ [t min , t max ] and no pole inserted between two consecutive roots (τ 2j -1 , τ 2j ), so that we can write: Here again, we have to define d N cuts on the real negative axis (figure 3) in view of a uniform representation of the square root in ( 23). The same mathematical treatment as above leads to the following expression of f (N) 2 (t):

f (N) * 2 (p) = µ * N (p) K (N) P (N) * (p) (p + 1/θ 1 ) • • • (p + 1/θ d N ) , K (N) = (B * N (0) 2 -4A * N (0)C * N (0))τ 1 • • • τ 2d N 2A * N (0)θ 1 • • • θ d N , P (N) * (p) = (p + 1/τ 1 ) • • • (p + 1/τ 2d N ). ( 23
f (N) 2 (t) = ± K (N) π d N j =1 (-1) j +1 -1/τ 2j -1/τ 2j-1 -P (N) (x)µ * N (x) (x + 1/θ 1 ) • • • (x + 1/θ d N ) e tx dx + N i=1 π -P (N) (-1/t n N i ) e -t/t n N i 2b n N i (1/θ 1 -1/t n N i ) • • • (1/θ d N -1/t n N i ) + d N j =1 (-1) j π -P (N) (-1/τ j )µ * N (τ j ) e -t/τ j (1/τ 1 -1/τ j ) • • • (1/τ j -1 -1/τ j )(1/τ j +1 -1/τ j ) • • • (1/τ d N -1/τ j ) (24)
with the sign ± chosen according to the condition µ eff (0) = lim |p|→∞ (µ eff (p)) > 0. After adding f (N) 1 (t) and f (N) 2 (t) and operating in (24) the change of variable τ = -1/x, µ eff(N) (t) can be put in the standard form (6). We get (d N + n N ) discrete lines with a positive or null amplitude at times θ j (j

∈ [1, d N ]) and t n N i (i ∈ [1, n N ]
), the latter ones corresponding to the matrix phase (N), and a continuous spectrum split into d N separate parts on the intervals [τ 2j -1 , τ 2j ] (j ∈ [1, d N ]) with the amplitude:

g (N) (τ ) = K (N) θ 1 • • • θ d N π τ 1 • • • τ 2d -(τ -τ 1 ) • • • (τ -τ 2d N ) τ (τ -θ 1 ) • • • (τ -θ d N ) n N k=1 t k N 2b k N (τ -t k N ) , τ ∈ τ 1 , τ 2 ∪ • • • ∪ τ 2j -1 , τ 2j ∪ • • • ∪ τ 2d N -1 , τ 2d N . ( 25
)
Note that in (25) too the matrix phase (N) plays a prominent role with respect to the other ones through the term

( n N k=1 t k N 2b k N (τ -t k N )
). Now again, it is obvious that the overall behaviour is no more a (generalized) Maxwellian one since the long range memory effect leads to continuous relaxation spectra. When one relaxation time is added to some phase, the number of separate parts of the continuous spectrum may be modified through the modification of d N and the shape of the whole spectrum may be completely changed. When one Maxwellian layer is added, this gives rise to two additional parts of the continuous spectrum and to two new lines. Among many possible applications, we have chosen two simple illustrative examples:

-on figure 4, we use the three-phase self-consistent scheme (N = 2) for the case already treated on figure 2b but we give an additional relaxation time to each phase. We find six (instead of two) separate parts for the continuous spectrum and six (instead of three) discrete lines with a non zero amplitude; it turns out that, if such a spectrum was got from experiment, it could be used to derive the total number of relaxation times only if the number of phases and their inclusion/matrix-type morphology is known. -on figure 5, we still start from the same example treated on figure 2b but we transform a small part of the matrix phase into a different Maxwellian material inserted between the matrix and the inclusion as a thin interphase layer, with the volume fraction c 2 = 1.5% and material constants intermediate between those of the two main phases. According to our model (with N = 3), the relaxation spectrum is modified drastically, exhibiting a well-identified signature of the interphase, whereas at the same time the overall resulting creep or relaxation functions would be almost unchanged. Here again, we can suggest from these illustrations that the spectral analysis of available experimental data on viscoelastic materials could provide useful information on the constitutive micromorphology of these materials when it is derived in connection with micromechanical models such as the one proposed hereabove.

ν 1 = ν 2 = 0.5; c = 0.5; b 1 1 = 10; b 2 1 = 5; b 1 2 = 1; b 2 2 = 2; t 1 1 = 10; t 2 1 = 5; t 1 2 = 1; t 2 2 = 2.

Bulk relaxation spectra

The bulk relaxation spectra may be derived much more easily since, in the elastic case [START_REF] Hervé | N-layered inclusion-based micromechanical modelling[END_REF], the overall bulk modulus k eff(N) is given in closed form instead of as the root of a quadratic equation. From that, it can be inferred that the resulting spectra in the viscoelastic case reduce to discrete lines when the constituent phases are (generalized)-Maxwellian bodies. When Poisson's ratios are constant in time, we start from the relation:

k * i (p) = 2(1 + ν i ) 3(1 -2ν i ) µ * i (p), i ∈ [1, N]. ( 26
)
Through the use of the elastic solution, of the correspondence principle and of the expression (19) for the shear moduli, we easily find (see Appendix 2) the Laplace-Carson-transformed effective bulk modulus k eff(N) * (p) in the form of the rational fraction

k eff(N) * (p) = S (N) * (p) G (N) * (p) µ * N (p), ( 27
)
k eff(N) (t) = n N k=1 F k e -t/t k N + d N /2 =1 E e -t/ξ (28)
with F k , E > 0. Thus the bulk relaxation spectrum reduces to discrete lines, as illustrated in figure 6 for a 5-phase model with four Maxwellian phases. Note that, for the bulk modulus only, the overall behaviour is a generalized Maxwellian one, but with relaxation times which, except for the matrix phase, do not coincide with those of the constituents.

Transverse shear relaxation spectra according to the '(N + 1)-phase cylindrical self-consistent model'

A straightforward extension of the foregoing analysis may be performed similarly to what has been done already in elasticity [START_REF] Hervé | Elastic behaviour of multiply coated fibre-reinforced composites[END_REF] by dealing with circular cylindrical instead of spherical geometries through the '(N + 1)-phase cylindrical self-consistent scheme' so as to consider multiply coated fibre-(instead of particle-) reinforced composites with phases obeying a transversely isotropic (generalized) Maxwellian behaviour (which of course includes the case of linear elasticity for elastic fibres). Since the method is unchanged, we only give in Appendix 3 the new expression of the polynoms A * N (p), B * N (p) and C * N (p) for the cylindrical geometry when the fibres lie along x 1 from which, through Equations ( 20)-( 25), one can easily derive the shear relaxation spectra for the transverse shear modulus µ eff(N) 23 (t) when the Poisson ratios

ν i = C i 23 (t)/[C i 23 (t) + C i 22 (t)
] actually do not depend on time t (here C i (t) is the fourth order relaxation tensor of phase (i)). The mathematical structure of these equations indicates that split continuous spectra will be obtained for the transverse shear modulus, whereas discrete lines only are to be expected for the four other moduli, as already known for the case N = 2 [START_REF] Bourgeois | Incompressible fibre-reinforced composites: Effective elastoplastic or viscoelastic behaviour[END_REF].

Conclusion

The foregoing analysis gives an additional unambiguous illustration of the importance of the 'long range memory effect' in viscoelastic heterogeneous materials, resulting in the fact that the overall behaviour is generally much more complex than the one of the constituents; moreover it leads to suggest to experimentalists in the field of polymer blends and reinforced polymers to make efforts in order to extract the relaxation spectra from their experimental data so as to get useful informations on the micromorphology of their materials (including the question of the presence of an interphase). Since a direct investigation of such spectra presently seems to be out of reach, this could be done by numerical deconvolution treatments applied to experimental creep/relaxation or complex moduli data.

In addition, the hereabove defined linear viscoelastic (N + 1)-phase model can be used for a simplified treatment of the nonlinear behaviour of polymer blends [START_REF] Beurthey | Modélisations du comportement d'alliages de polymères[END_REF] or composites, in a way similar to what has already been done in the field of elastoplasticity [START_REF] Bornert | Self-consistent approaches and strain heterogeneities in two-phase elastoplastic materials[END_REF][START_REF] Zaoui | Structural morphology and constitutive behaviour of microheterogeneous materials[END_REF].

A * N (p) = 4R 10 N (1 -2ν N )(7 -10ν N )Z * 12 + 20R 7 N 7 -12ν N -8ν 2 N Z * 42 + 12R 5 N (1 -2ν N ) Z * 14 -7Z * 23 + 20R 3 N (1 -2ν N ) 2 Z * 13 + 16(4 -5ν N )(1 -2ν N )Z * 43 , B * N (p) = 3R 10 N (1 -2ν N )(15ν N -7)Z * 12 + 60R 7 N (ν N -3)ν N Z * 42 -24R 5 N (1 -2ν N ) Z * 14 -7Z * 23 -40R 3 N (1 -2ν N ) 2 Z * 13 -8(1 -5ν N )(1 -2ν N )Z * 43 , C * N (p) = -R 10 N (1 -2ν N )(7 + 5ν N )Z * 12 + 10R 7 N 7 -ν 2 N Z * 42 + 12R 5 N (1 -2ν N ) Z * 14 -7Z * 23 + 20R 3 N (1 -2ν N ) 2 Z * 13 -8(7 -5ν N )(1 -2ν N )Z * 43 with Z * αβ = L (N-1) * α1 L (N-1) * β2 -L (N-1) * β1 L (N-1) * α2 , (α, β) ∈ [1, 4], L (N-1) * = N-1 i=1 M (i) *
and

M (i) * =                c i 3 R 2 i (3b i -7c i ) 5(1 -2ν i ) -12α i R 5 i 4(f i -27α i ) 15(1 -2ν i )R 3 i 0 (1 -2ν i+1 )b i 7(1 -2ν i ) -20(1 -2ν i+1 )α i 7R 7 i -12α i (1 -2ν i+1 ) 7(1 -2ν i )R 5 i R 5 i α i 2 -R 7 i (2a i + 147α i ) 70(1 -2ν i ) d i 7 R 2 i (g i + 12α i (7 -10ν i+1 ) -7e i ) 35(1 -2ν i ) - 5 6 (1 -2ν i+1 )α i R 3 i 7(1 -2ν i+1 )α i R 5 i 2(1 -2ν i ) 0 e i (1 -2ν i+1 ) 3(1 -2ν i )                , a i = -(7 -10ν i )(7 + 5ν i+1 ) p + 1/t j +1 i • • • p + 1/t n i i N i+1 (p) +(7 + 5ν i )(7 -10ν i+1 ) p + 1/t j +1 i+1 • • • p + 1/t n i+1 i+1 N i (p), b i = 4(7 -10ν i ) p + 1/t j +1 i • • • p + 1/t n i i N i+1 (p) + (7 + 5ν i ) p + 1/t j +1 i+1 • • • p + 1/t n i+1 i+1 N i (p), c i = (7 -5ν i+1 ) p + 1/t j +1 i • • • p + 1/t n i i N i+1 (p) + 2(4 -5ν i+1 ) p + 1/t j +1 i+1 • • • p + 1/t n i+1 i+1 N i (p), d i = (7 + 5ν i+1 ) p + 1/t j +1 i • • • p + 1/t n i i N i+1 (p) + 4(7 -10ν i+1 ) p + 1/t j +1 i+1 • • • p + 1/t n i+1 i+1 N i (p), e i = 2(4 -5ν i ) p + 1/t j +1 i • • • p + 1/t n i i N i+1 (p) + (7 -5ν i ) p + 1/t j +1 i+1 • • • p + 1/t n i+1 i+1 N i (p), f i = (4 -5ν i ) p + 1/t j +1 i • • • p + 1/t n i i N i+1 (p) -(4 -5ν i+1 )(7 -5ν k ) p + 1/t j +1 i+1 • • • p + 1/t n i+1 i+1 N i (p), g * i = 105(1 -ν i+1 ) p + 1/t j +1 i • • • p + 1/t n i i N i+1 (p), α i = p + 1/t j +1 i+1 • • • p + 1/t n i+1 i+1 N i (p) -p + 1/t j +1 i • • • p + 1/t n i i N i+1 (

p). • Bulk relation spectra

Polynoms S (N) * (p) and G (N) * (p) in ( 27) are given by: and

Q (N-1) * = N-1 i=1 T * i (R i ), T (i) * (R i ) =       2 1 + ν i 1 -2ν i 1 * i + 4.h * i 4 R 3 i h * i -1 * i 2R 3 i 1 + ν i+1 1 -2ν i+1 h * i - 1 + ν i 1 -2ν i 1 * i 2 1 + ν i+1 1 -2ν i+1 h * i + 4.l * i       , h * i = p + 1/t j +1 i • • • p + 1/t n i i T i+1 (p), S * (p) = 2 (1 + ν N ) (1 -2ν N ) R 3 N Q N-1 * 11 -4Q N-1 * 21 l * i = p + 1/t j +1 i+1 • • • p + 1/t n i+1 i+1 T i (p), G * (p) = 3 R 3 N Q N-1 * 11 + Q N-1 *
M (i) * =          a i q 2 i 1 q 2 i ρ i 1 q 2 k ρ i 0 2q 4 i b i q 4 i c i 2q 4 i d i -q 4 i ρ i -3q 2 i ρ i 0 q 2 i a k q 2 i ρ i -6d i 3ρ i 2b i c i          , a i = p + 1/t j +1 i+1 • • • p + 1/t n i+1 i+1 N i (p) + (3 -4ν i ) p + 1/t j +1 i • • • p + 1/t n i i N i+1 (p), b i = (2ν k+1 -3) p + 1/t j +1 i+1 • • • p + 1/t n i+1 i+1 N i (p) + (3 -2ν k ) p + 1/t j +1 i • • • p + 1/t n i i N i+1 (p), c i = (3 -4ν k+1 ) p + 1/t j +1 i+1 • • • p + 1/t n i+1 i+1 N i (p) + p + 1/t j +1 i • • • p + 1/t n i i N i+1 (p), d i = (1 -2ν k+1 ) p + 1/t j +1 i+1 • • • p + 1/t n i+1 i+1 N i (p) + (2ν k -1) p + 1/t j +1 i • • • p + 1/t n i i N i+1 (p), ρ i = p + 1/t j +1 i • • • p + 1/t n i i N i+1 (p) -p + 1/t j +1 i+1 • • • p + 1/t n i+1
i+1 N i (p),

q i = R i R i+1 .

Figure 1 .

 1 Figure 1. Integration paths used for the inversion of the Laplace transformation (t 1 < t 2 ); (a) C.S.C.S; (b) G.S.C.S.

Figure 2 .

 2 Figure 2. Comparison of the relaxation spectra as derived from the two considered models, in the case c = 0.5; b 1 = 1; b 2 = 10; t 1 = 1; t 2 = 10; (a) classical self-consistent model; (b) generalized self-consistent model.

Figure 3 .

 3 Figure 3. (a) Integration path used for the Laplace-Carson inversion; (b) roots location for (1/t i < 1/t i+1 ).

Figure 4 .

 4 Figure 4. Shear relaxation spectrum derived from the 3-phase model for phases with two relaxation times: ν 1 = ν 2 = 0.5; c = 0.5; b 1 1 = 10; b 2 1 = 5; b1 2 = 1; b 2 2 = 2; t 1 1 = 10; t 2 1 = 5; t 1 2 = 1; t 2 2 = 2.

Figure 5 .

 5 Figure 5. Shear relaxation spectrum derived from the 4-phase self consistent model, in the case of an interphase: ν i = 0.5; c 1 = 0.5; c 2 = 0.015; c 3 = 0.485; b 1 = t 1 = 10; b 2 = t 2 = 5; b 3 = t 3 = 1.

Figure 6 .

 6 Figure 6. Bulk relaxation spectrum of bulk modulus from the 5-phase self-consistent model (c 1 = 0.5; c 2 = 0.165; c 3 = 0.2; c 4 = 0.135; t 1 = b 1 = 10; t 2 = 2; t 3 = 1.5; t 4 = b 4 = 1; b 2 = b 3 = 1).

Figure A1 .

 A1 Figure A1. The N -layered spherical inclusion embedded in an infinite matrix.

Acknowledgements

Yves Rougier and Claude Stolz are gratefully acknowledged for their initial contribution to this study.

Appendix 1

• The incompressible G.S.C.S.

Coefficients A, B abd C in (7) for the elastic incompressible G.S.C.S. [START_REF] Hervé | Modelling the effective behaviour of nonlinear matrix inclusion composites[END_REF] read:

• Shear relaxation spectrum (G.S.C.S.)

Coefficient k G in (13) reads:

Appendix 2

• Shear relaxation spectra

In (19), we suppose the j first relaxation times of phases (i) and (i + 1) to be equal (j = j i,i+1 20) are then given by: