Jean-François Barthélémy
email: jf.barthelemy@cerema.fr

Albert Giraud
email: albert.giraud@univ-lorraine.fr

Francis Lavergne
email: francis.lavergne@cerema.fr

Julien Sanahuja
email: julien.sanahuja@edf.fr

François Barthélémy

The Eshelby inclusion problem in ageing linear viscoelasticity

Keywords: Micromechanics, Ageing linear viscoelasticity, Viscoelastic Eshelby tensor, Viscoelastic Hill tensor, Viscoelastic inclusion problem

The present paper focuses on the Eshelby inclusion problem which is revisited in the framework of ageing linear viscoelasticity. All known results established in linear viscoelasticity thanks to the correspondence principle are recovered as particular cases of a general solution extended to ageing. A closed form solution is presented for the Hill and Eshelby tensors related to an ellipsoidal inclusion embedded in an infinite anisotropic medium. The case of the isotropic medium is investigated in detail and related solutions are presented for spherical and ellipsoidal inclusions. A numerical procedure which operates in time domain and based on the trapezoidal rule to determine Volterra integral operators allows to efficiently calculate Hill and Eshelby tensors for a wide range of behaviours including in particular time-dependent Poisson ratios. Validation and verification of the new developed solution are presented for the ageing spherical inclusion on the basis of recent published results and it is completed by comparisons with finite element simulations for the general and novel case of an oblate spheroid embedded in an ageing linear viscoelastic matrix.

Introduction

This paper revisits the Eshelby inclusion problem [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] in the framework of ageing linear viscoelasticity, on the basis of the ageing correspondence principle proposed by [START_REF] Read | Stress Analysis for Compressible Viscoelastic Materials[END_REF], [START_REF] Lee | Stress analysis in linear viscoelastic materials[END_REF], [START_REF] Lee | Stress analysis for linear viscoelastic materials[END_REF], [START_REF] Mandel | Sur les corps viscoélastiques linéaires dont les propriétés dépendent de l'âge[END_REF] and [START_REF] Mandel | Cours de Mécanique des Milieux Continus[END_REF]. It must be emphasized that it includes as a particular case the non-ageing linear viscoelasticity and all results developed in this context may be recovered. A novel solution is presented for the Hill polarization tensor, and the Eshelby tensor related to an ellipsoidal inclusion embedded in an ageing linear viscoelastic material and it extends results presented in recent papers ( [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF], [START_REF] Lavergne | Homogenization schemes for aging linear viscoelastic matrix-inclusion composite materials with elongated inclusions[END_REF]).

In the context of non-ageing linear viscoelasticity, and by using the correspondence principle, [START_REF] Hashin | Viscoelastic Behavior of Heterogeneous Media[END_REF], [START_REF] Hashin | Complex moduli of viscoelastic composites -I. General theory and application to particulate composites[END_REF], [START_REF] Hashin | Complex moduli of viscoelastic composites -II. Fiber reinforced materials[END_REF], [START_REF] Christensen | Viscoelastic properties of heterogeneous media[END_REF] and [START_REF] Christensen | Theory of Viscoelasticity : An Introduction, 2nd Edition[END_REF] have presented generalization of linear elastic problem to the linear viscoelastic one, using Stieltjes integrals combined with Laplace-Carson transform. Such an approach has been widely used to extend homogenization schemes based on Eshelby tensors in elasticity to linear viscoelasticity (see reviews of [START_REF] Pichler | Upscaling of viscoelastic properties of highly-filled composites: Investigation of matrix-inclusion-type morphologies with power-law viscoelastic material response[END_REF], [START_REF] Lavergne | Estimation of creep strain and creep failure of a glass reinforced plastic by semi-analytical methods and 3D numerical simulations[END_REF], [START_REF] Berbenni | A new internal variables homogenization scheme for linear viscoelastic materials based on an exact Eshelby interaction law[END_REF], and [START_REF] Sevostianov | Effective properties of linear viscoelastic microcracked materials: Application of Maxwell homogenization scheme[END_REF]). In many cases, exact solutions may be obtained in Laplace Transform Domain (TD, see [START_REF] Wang | The Influence of Inclusion Shape on the Overall Viscoelastic Behavior of Composites[END_REF] and [START_REF] Li | Anisotropic stress-strain relations and complex moduli of a viscoelastic composite with aligned spheroidal inclusions[END_REF]) but, due to the complexity of relations in TD, analytical inversions of Laplace transforms are only restricted to simple cases. Hashin has presented in [START_REF] Hashin | The inelastic inclusion problem[END_REF] an analytical transient solution in the case of a spherical isotropic elastic inclusion embedded in a non compressible isotropic viscoelastic matrix. One may find in [START_REF] Laws | Self-Consistent Estimates for the Viscoelastic Creep Compliances of Composite Materials[END_REF] the complete linear viscoelastic solution for ellipsoidal inclusions in anisotropic materials, and numerical results obtained thanks to a numerical inversion of Laplace-Carson solutions performed by using collocation method ( [START_REF] Schapery | Approximate methods of transform inversion for viscoelastic stress analysis[END_REF] and [START_REF] Rekik | Optimization of the collocation inversion method for the linear viscoelastic homogenization[END_REF] for a recent extension). Exact inversion of Laplace transform has been presented in [START_REF] Rougier | Représentation spectrale en viscoélasticité linéaire des matériaux hétérogènes[END_REF] for the problem of a spherical inclusion, considering isotropic non compressible viscoelastic Maxwell materials for inclusion and matrix and self consistent homogenization scheme. It may be noticed that, still in the context of non-ageing linear viscoelasticity, fraction-exponential operators which are extensively used in rheology and linear viscoelasticity (see [START_REF] Mainardi | Fractional Calculus and Waves in Linear Viscoelasticity. An Introduction to Mathematical Models[END_REF]) have been applied to Eshelby inclusion problem and related homogenization approach, to obtain solutions in more general cases of ellipsoidal or penny-shaped cracks in a viscoelastic matrix (see [START_REF] Levin | Effective properties of viscoelastic media with crack-like inclusions[END_REF], [START_REF] Sevostianov | Effective properties of linear viscoelastic microcracked materials: Application of Maxwell homogenization scheme[END_REF] and [START_REF] Sevostianov | Creep and relaxation contribution tensors for spheroidal pores in hereditary solids: fraction-exponential operators approach[END_REF]). Such operators have the advantage to accurately fit a wide range of behaviours and to conduce to analytical expressions for inverse Laplace transform.

Although theoretically very useful to solve non-ageing linear viscoelastic problems thanks to the correspondence principle, the Laplace-Carson transform suffers from important drawbacks mainly due to the fact that it remains difficult to ensure both accuracy and stability of the numerical algorithms available for inversion. Indeed the method based on the correspondence principle first requires to calculate transforms of some creep or relaxation functions of different phases but closed-form expressions may not exist due to the complexity of the initial functions or because the latter come from experimental points. The last step of the method consists in calculating the inverse of Laplace-Carson transforms in order to come back to the time domain. This step almost systematically needs to resort to numerical algorithms as those described in [START_REF] Abate | A Unified Framework for Numerically Inversion Laplace[END_REF] which generally have to be specifically adapted to the function to inverse. It proves in particular very difficult to reach a certain level of accuracy of the inverse while reducing the sensitivity to the errors in the Laplace domain. To avoid the shortcomings of the Laplace-Carson transform, homogenization methods staying in the time domain have recently been designed ( [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF], [START_REF] Masson | Incremental homogenization approach for ageing viscoelastic polycrystals[END_REF], [START_REF] Berbenni | A new internal variables homogenization scheme for linear viscoelastic materials based on an exact Eshelby interaction law[END_REF]). Among them, one may cite the time-incremental internal variable homogenization method based on Green function techniques and related integral representation presented in [START_REF] Coulibaly | Micromechanical modeling of linear viscoelastic behavior of heterogeneous materials[END_REF], [START_REF] Dinzart | Homogenization of the viscoelastic heterogeneous materials with multi-coated reinforcements: an internal variables formulation[END_REF] and [START_REF] Berbenni | A new internal variables homogenization scheme for linear viscoelastic materials based on an exact Eshelby interaction law[END_REF]. In this method, the effective behaviour as well as the evolution laws of the averaged stresses per phase are solved incrementally in the time domain through a time-differential equation, without need to analytical or numerical inverse of Laplace-Carson transforms.

Regarding ageing viscoelastic composite materials, some results have been obtained in the case of periodic homogenization applied to multilayered media in [START_REF] Maghous | Periodic homogenization in thermoviscoelasticity: case of multilayered media with ageing[END_REF] but mean-field schemes based on Eshelby's solution were previously restricted to particular ageing behaviours for the viscoelastic phases, such as the time-shift behaviour introduced to describe plastic materials and glasses ( [START_REF] Struik | Physical aging in amorphous polymers and other materials[END_REF], [START_REF] Sullivan | Creep and physical aging of composites[END_REF], [START_REF] Odegard | Physical aging of epoxy polymers and their composites[END_REF]). In this particular case, the correspondence principle can be applied after the change of variable defining the equivalent time ( [START_REF] Zheng | A new constitutive equation for the long-term creep of polymers based on physical aging[END_REF], [START_REF] Lavergne | Estimation of creep strain and creep failure of a glass reinforced plastic by semi-analytical methods and 3D numerical simulations[END_REF]). More recently, comprehensive homogenization methods have been uncovered. On the one hand, the incremental homogenization approach [START_REF] Masson | Incremental homogenization approach for ageing viscoelastic polycrystals[END_REF] consists in approximating the ageing viscoelastic behaviour to a non-ageing Dirichlet series on each time step. The correspondence principle is then applied on each time step to solve the homogenization problem. The scheme designed in [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF] is also able to estimate the effective response of viscoelastic composites: it relies on the variational technique proposed in [START_REF] Castañeda | New variational principles in plasticity and their application to composite materials[END_REF] to rigorously approximate non-uniform eigenstrains due to internal variables by piecewise uniform eigenstrains. On the other hand, [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] has proposed a closed-form solution to handle the case of spherical inhomogeneities, the ageing viscoelastic behaviours of all phases being isotropic without any particular restriction on the bulk and shear relaxation functions. The Volterra operator and its inverse appear in this solution. Using the properties of the Volterra operator, the viscoelastic Eshelby inclusion problem has been solved for ellipsoidal inclusions [START_REF] Lavergne | Homogenization schemes for aging linear viscoelastic matrix-inclusion composite materials with elongated inclusions[END_REF], providing that the ageing viscoelastic material features a time-independent Poisson ratio in a sense defined in [START_REF] Hilton | The significance of (an)isotropic viscoelastic Poisson ratio stress and time dependencies[END_REF]. Nevertheless, the viscoelastic Eshelby inclusion problem for ellipsoidal inclusions and an ageing viscoelastic material featuring a time-dependent Poisson ratio could only be solved by mean of approximations ( [START_REF] Masson | Incremental homogenization approach for ageing viscoelastic polycrystals[END_REF], [START_REF] Lavergne | Predicting the aging viscoelastic behavior of fiber-reinforced concrete by a multi-scale homogenization scheme[END_REF]) and an universal closed-form solution to the viscoelastic Eshelby inclusion problem was yet to be found. All the references listed in the present paragraph share a common feature: the strains in the Eshelby inclusion, or in the inhomogeneity, remain uniform at all times, as for the solutions to the elastic and the non-ageing viscoelastic Eshelby inclusion problems.

In the present article, after presenting the hypotheses of the ageing linear viscoelastic inclusion problem, a comprehensive closed-form solution is disclosed. Then, its expression in the case of an isotropic infinite medium is derived and successfully compared to existing results in particular situations. In the last section, the efficient numerical procedure designed in [START_REF] Bazant | Numerical determination of long-range stress history from strain history in concrete[END_REF], [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] and [START_REF] Lavergne | Homogenization schemes for aging linear viscoelastic matrix-inclusion composite materials with elongated inclusions[END_REF] is applied for further verifications.

The problem of an ellipsoidal inclusion of arbitrary shape in an ageing linear viscoelastic medium

The ageing linear viscoelastic behaviour

As in [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF], the studied domain is the whole three-dimensional space and the hypothesis of small perturbation is adopted. The difference with [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] relies on the fact that the constitutive law of the material is of the ageing linear viscoelastic type. This means that the strain and stress tensor histories, respectively ε(t) and σ(t), are related by means of a Stieltjes integral [START_REF] Salençon | Viscoélasticité pour le calcul des structures[END_REF]

ε(t) = t t ′ =-∞ L(t, t ′ ) : dσ(t ′ ) = t t ′ =-∞ L(t, t ′ ) : σ(t ′ ) dt ′ (1)
where L is the creep compliance tensor of the fourth order. In the general ageing framework, it depends on two independant time variables and is such that L(t, t ′ ) = 0 if t < t ′ due to the causality principle. The non-ageing case corresponds to a dependence of t and t ′ only through their difference tt ′ . The resolution of the Eshelby problem in the latter case is already available thanks to the correspondence principle and the resolution of the associated elastic problem in the Laplace-Carson domain (see references in introduction). The present paper deals with the more general ageing case.

The linear relationship (1), more simply denoted by ε = L

• : σ ([7]
, [START_REF] Lavergne | Homogenization schemes for aging linear viscoelastic matrix-inclusion composite materials with elongated inclusions[END_REF]), generalizes to the tensorial framework the so-called Volterra operator "•" [START_REF] Volterra | Theory of Functionals and of Integral and Integro-differential equations[END_REF] between scalar functions

Y = F • X ⇔ Y(t) = t t ′ =-∞ F(t, t ′ ) dX(t ′ ) (2)
When explicit references to the time variables are needed for clarity, (1) can also be written

Y(t) = F(t, •) • X(•)
where "•" represents the dummy variable involved in the Volterra integration. When a second dummy variable is needed in any expression in the following, the symbol "⋆" will be used.

Adopting the convention of summation of repeated indices, (1) writes in components

ε = L • : σ ⇔ ε i j = L i jkl • σ kl (3)
Note that L obviously verifies the minor symmetries L jikl = L i jlk = L i jkl .

Similarly a Volterra operator can be built between a second-order kernel K(t, t ′ ) and a time-dependent vector u(t)

T = K • • u ⇔ T(t) = t t ′ =-∞ K(t, t ′ ) • du(t ′ ) (4)
It is worth recalling that time discontinuities are allowed in all the functions involved in the previous relationships and time derivatives must be considered in the framework of the theory of distributions as introduced by [START_REF] Schwartz | Théorie des distributions[END_REF] (see also [START_REF] Gel'fand | Generalized functions -Vol I Properties and Operations[END_REF]).

The scalar Volterra kernel (t, t ′ ) → H(tt ′ ) where H is the Heaviside function 1 operates as an identity in [START_REF] Read | Stress Analysis for Compressible Viscoelastic Materials[END_REF]. For the sake of commodity, H shall also denote this function of two time variables in the following, i.e. H • X = X. As a consequence, the identities of (1) and ( 4) are respectively HI and H1, where I is the fourth-order identity tensor acting over symmetric second-order tensors (I i jkl = (δ ik δ jl + δ il δ jk )/2) and 1 is the second-order identity tensor (1 i j = δ i j ) with δ i j denoting the Kronecker symbol2 .

The operations between a function of two time variables and a function of a single time variable can easily be extended to associative, distributive with respect to the addition, but non-commutative operations between two functions of two time variables, either in scalar or tensorial case. As stated in [START_REF] Maghous | Periodic homogenization in thermoviscoelasticity: case of multilayered media with ageing[END_REF], the commutativity between two scalar functions of two time variables holds only for non-ageing kernels. A fourth-order tensorial example writes (A

• : B)(t, t ′ ) = A(t, •) • : B(•, t ′ ) = t τ=-∞ A(t, τ) : ∂B ∂τ (τ, t ′ ) dτ (5)
The definition of the inverse of a (scalar or tensorial) kernel in the sense of Volterra follows immediately. To avoid any confusion with the classical inversion, the notation T -1 • for the inverse of the (scalar or tensorial) kernel T is adopted. For instance, the inverse of the creep compliance L, called the relaxation tensor and denoted by

C = L -1 • , satisfies C • : L = L • : C = HI.

The Eshelby inclusion problem

Let E denote an ellipsoid in R3 defined by three radii ρ 1 , ρ 2 and ρ 3 (alternatively denoted by a, b, and c when useful for readability enhancement of formulas in Appendix C) and an associated orthonormal frame (e 1 , e 2 , e 3 ) such that the equation of E writes

A -1 • x ≤ 1 with A = ρ 1 e 1 ⊗ e 1 + ρ 2 e 2 ⊗ e 2 + ρ 3 e 3 ⊗ e 3 (6)
A time-dependent polarization field p(x, t) = p(t)χ E (x) is then introduced, where χ E is the characteristic function of E 3 . In other words, p(x, t) is uniform within E and 0 outside E. The whole space is occupied by a uniform ageing linear viscoelastic material of relaxation tensor C so that the constitutive law is expressed as

σ(x, t) = C(t, •) • : ε(x, •) + p(t)χ E (x) (7)
In absence of any body force, the momentum balance applied to [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF], i.e. div σ = 0 involves the term div ( pχ E ). Recalling that p does not depend on the position x and that grad χ E = -nδ ∂E where δ ∂E is the surface Dirac distribution ( [START_REF] Schwartz | Théorie des distributions[END_REF], [START_REF] Gel'fand | Generalized functions -Vol I Properties and Operations[END_REF]) supported by the boundary of the ellipsoid E and n is its unit outward normal vector, it comes that div

( pχ E ) = p • grad χ E = -p • nδ ∂E .
The displacement field u(x, t) solution to the ageing linear viscoelastic inclusion problem must finally satisfy the following equations

           div C(t, •) • : grad u(x, •) -p(t) • n(x) δ ∂E (x) = 0 lim x →∞ u(x, t) = 0 (8a) (8b)
which is rigorously the same as in the elastic case provided that the double contracted product ":" be replaced by the tensorial Volterra operator "

• :".

Resolution and construction of ageing viscoelastic Hill and Eshelby tensors

Elementary problems and Green tensor

This section aims at deriving the complete solution of the problem (8a)-(8b) thanks to a methodology adapted from the elastic case presented in [START_REF] Willis | Bounds and self-consistent estimates for the overall properties of anisotropic composites[END_REF] or [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]. For this purpose, it is useful to take advantage of the linearity of the problem and introduce the second-order ageing viscoelastic Green tensor G(x, t, t ′ ) associated to these equations: G(x, t, t ′ ) • P is the displacement field at time t satisfying (8b) and resulting from the application of a point force step P concentrated at the origin and occuring at t ′ , which means

div C(t, •) • : grad G(x, •, t ′ ) • P + H(t -t ′ )δ(x)P = 0 (9)
where δ stands here for the three-dimensional Dirac distribution associated with the origin ( [START_REF] Schwartz | Théorie des distributions[END_REF], [START_REF] Gel'fand | Generalized functions -Vol I Properties and Operations[END_REF]). It is worth precising again that t ′ is arbitrarily chosen as the starting time of the loading and the equation ( 9) finally depends on x and t as well as t ′ . Moreover, due to the fact that the problem is posed in an infinite domain, any translation of the loading, i.e. H(tt ′ )δ(xx ′ )P, results in the same translation in the first argument of the solution G, i.e. G(xx ′ , t, t ′ ) • P, which justifies the location of the elementary loading at the origin. In contrast, if the domain was bounded, it would have been necessary to consider unit loadings at any point x ′ and introduce this additional variable to the Green function.

After eliminating P, [START_REF] Hashin | Viscoelastic Behavior of Heterogeneous Media[END_REF] writes in components

C i jkl (t, •) • G kn, jl (x, •, t ′ ) + H(t -t ′ )δ in δ(x) = 0 ( 10
)
where indices following the comma ( j and l) represent derivatives with respect to the corresponding components of x (x j and x l ).

Two equivalent strategies can then be adoped to solve [START_REF] Hashin | Complex moduli of viscoelastic composites -I. General theory and application to particulate composites[END_REF]: one based on Fourier transforms as in [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] and another one, inspired by [START_REF] Willis | Bounds and self-consistent estimates for the overall properties of anisotropic composites[END_REF] and followed in the sequel, based on the plane-wave expansion (see Appendix A for a brief proof)

δ(x) = - 1 8π 2 ξ =1 δ ′′ (ξ • x) dS ξ (11)
Note that δ in the left hand side is the three-dimensional Dirac distribution wheras δ ′′ in the right hand side corresponds to the second derivative of the scalar Dirac distribution ( [START_REF] Schwartz | Théorie des distributions[END_REF], [START_REF] Gel'fand | Generalized functions -Vol I Properties and Operations[END_REF]).

Considering [START_REF] Hashin | Complex moduli of viscoelastic composites -II. Fiber reinforced materials[END_REF] as a linear decomposition, it becomes interesting to find the solutions G ξ (x, t, t ′ ) of the following equations parametrized by the unit vector ξ

C i jkl (t, •) • G ξ kn, jl (x, •, t ′ ) + H(t -t ′ )δ in δ ′′ (ξ • x) = 0 (12)
such that the relationship between G and G ξ stems from the superposition principle

G(x, t, t ′ ) = - 1 8π 2 ξ =1 G ξ (x, t, t ′ ) dS ξ (13)
The solution to [START_REF] Christensen | Viscoelastic properties of heterogeneous media[END_REF] can easily be built by analogy with the elastic case

G ξ (x, t, t ′ ) = -δ(ξ • x) ξ • C • ξ -1 • (t, t ′ ) = -δ(ξ • x) N ξ (t, t ′ ) ( 14
)
where the viscoelastic counterpart of the elastic acoustic tensor is identified in [START_REF] Pichler | Upscaling of viscoelastic properties of highly-filled composites: Investigation of matrix-inclusion-type morphologies with power-law viscoelastic material response[END_REF] by

K ξ (t, t ′ ) = ξ • C(t, t ′ ) • ξ and its inverse (in the sense of Volterra) by N ξ (t, t ′ ) = K ξ -1 • (t, t ′ ) = ξ • C • ξ -1 • (t, t ′ ).
Both are second-order Volterra kernels. In order to verify that G ξ given by ( 14) satisfies [START_REF] Christensen | Viscoelastic properties of heterogeneous media[END_REF], the derivatives of its components are calculated

G ξ kn, jl (x, t, t ′ ) = -ξ j ξ l δ ′′ (ξ • x) N ξ kn (t, t ′ ) ( 15
)
which successively implies

C i jkl (t, •) • G ξ kn, jl (x, •, t ′ ) = -ξ j C i jkl (t, •) ξ l • N ξ kn (•, t ′ ) δ ′′ (ξ • x) (16) = -K ξ ik (t, •) • N ξ kn (•, t ′ ) δ ′′ (ξ • x) (17) = -H(t -t ′ ) δ in δ ′′ (ξ • x) ( 18
)
The definition of the inverse in the sense of Volterra N ξ = K ξ -1 • in (17) leads to [START_REF] Wang | The Influence of Inclusion Shape on the Overall Viscoelastic Behavior of Composites[END_REF] which is indeed identical to [START_REF] Christensen | Viscoelastic properties of heterogeneous media[END_REF].

The determination of G is finally completed by introducing ( 14) in ( 13)

G(x, t, t ′ ) = 1 8π 2 ξ =1 δ(ξ • x) N ξ (t, t ′ ) dS ξ (19)
It should also be noticed that [START_REF] Li | Anisotropic stress-strain relations and complex moduli of a viscoelastic composite with aligned spheroidal inclusions[END_REF] allows to verify that lim

x →∞ G(x, t, t ′ ) = 0 since δ(ξ • x) = δ(ξ • e)/
x where e = x/ x .

Solution to the Eshelby inclusion problem

The solution to (8a)-( 8b) is obtained by remarking that the loading -p(t) • n(x)δ ∂E (x) can be viewed as a superposition of translated elementary Dirac contributions. Consequently the displacement field writes (omitting time dependences)

u(x, t) = - x ′ ∈∂E G(x -x ′ , t, •) • • p(•) • n(x ′ ) dS x ′ (20)
Note that (8b) is satisfied thanks to the behaviour at infinity of the Green kernel.

Using [START_REF] Li | Anisotropic stress-strain relations and complex moduli of a viscoelastic composite with aligned spheroidal inclusions[END_REF], [START_REF] Hashin | The inelastic inclusion problem[END_REF] becomes

u(x, t) = - 1 8π 2 x ′ ∈∂E ξ =1 δ ξ • (x -x ′ ) N ξ (t, •) • • p(•) • n(x ′ ) dS ξ dS x ′ (21)
Equivalently, [START_REF] Laws | Self-Consistent Estimates for the Viscoelastic Creep Compliances of Composite Materials[END_REF] yields in components

u i (x, t) = - 1 8π 2 x ′ ∈∂E ξ =1 δ ξ • (x -x ′ ) N ξ ik (t, •)n l (x ′ ) dS ξ dS x ′ • p kl (•) (22)
Applying Fubini's theorem and Stokes' theorem to the integral over ∂E leads to

u i (x, t) = 1 8π 2 ξ =1 x ′ ∈E δ ′ ξ • (x -x ′ ) N ξ ik (t, •)ξ l dΩ x ′ dS ξ • p kl (•) (23)
Finally the strain tensor is given by a relationship analogous to the elastic case ε(x, t) = -P(x, t, •)

• : p(•) [START_REF] Rougier | Représentation spectrale en viscoélasticité linéaire des matériaux hétérogènes[END_REF] in which the fourth-order Volterra kernel P(x, t, t ′ ) writes

P(x, t, t ′ ) = - 1 8π 2 ξ =1 x ′ ∈E δ ′′ ξ • (x -x ′ ) Γ(ξ, t, t ′ ) dΩ x ′ dS ξ (25) with Γ(ξ, t, t ′ ) = ξ s ⊗ N ξ (t, t ′ ) s ⊗ ξ = ξ s ⊗ ξ • C • ξ -1 • (t, t ′ ) s ⊗ ξ (26)
where the notation s ⊗ denotes a tensor product with symmetrization over the two adjacent indices (here, on the one hand the two first and on the other hand the two last indices of Γ(ξ, t, t ′ )).

The end of the reasoning relies on the observation that the dependence of P on x in [START_REF] Mainardi | Fractional Calculus and Waves in Linear Viscoelasticity. An Introduction to Mathematical Models[END_REF] does not differ from that of the elastic case. In particular, the integrations are simplified thanks to the following change of variables

x ′ = A • y ′ and ξ = A -1 • ζ A -1 • ζ (27)
implying (see Appendix B for a proof of the second equation)

dΩ x ′ = det A dΩ y ′ and dS ξ = det A -1 A -1 • ζ 3 dS ζ (28)
and consequently

dΩ x ′ dS ξ = dΩ y ′ dS ζ A -1 • ζ 3 ( 29
)
Introducing ( 27) and ( 29) in [START_REF] Mainardi | Fractional Calculus and Waves in Linear Viscoelasticity. An Introduction to Mathematical Models[END_REF], P(x, t, t ′ ) becomes

P(x, t, t ′ ) = - 1 8π 2 ζ =1 Γ        A -1 • ζ A -1 • ζ , t, t ′        y ′ ≤1 δ ′′        ζ • ( A -1 • x -y ′ ) A -1 • ζ        dΩ y ′ dS ζ A -1 • ζ 3 ( 30
)
Observing that Γ (26) is positively homogeneous of degree 0 and that δ ′′ (z/λ) = λ 3 δ ′′ (z) for λ > 0 (obtained by differentiating twice δ(z/λ) = λδ(z)), P (30) is finally simplified

P(x, t, t ′ ) = - 1 8π 2 ζ =1 Γ( A -1 • ζ, t, t ′ ) y ′ ≤1 δ ′′ ζ • ( A -1 • x -y ′ ) dΩ y ′ dS ζ (31)
The integration over the unit ball with respect to the variable y ′ is then performed by means of cylindrical coordinates (r, θ, z) around the axis directed by ζ (i.e. z = ζ • y ′ ), which means dΩ y ′ = r dr dθ dz with θ varying from 0 to 2π, z from -1 to 1 and r from 0 to √ 1z 2 for each z. Hence [START_REF] Coulibaly | Micromechanical modeling of linear viscoelastic behavior of heterogeneous materials[END_REF] becomes

P(x, t, t ′ ) = - 1 8π ζ =1 Γ( A -1 • ζ, t, t ′ ) 1 z=-1 (1 -z 2 )δ ′′ (ζ • A -1 • x -z) dz dS ζ (32)
which can also write by application of the derivative rules of distributions ( [START_REF] Schwartz | Théorie des distributions[END_REF], [START_REF] Gel'fand | Generalized functions -Vol I Properties and Operations[END_REF])

P(x, t, t ′ ) = 1 4π ζ =1 Γ( A -1 • ζ, t, t ′ ) 1 z=-1 δ(ζ • A -1 • x -z) dz dS ζ (33)
This expression still depends on x. Nevertheless, similarly to the elastic case, if x belongs to the ellipsoid E of equation ( 6), the scalar ζ • A -1 • x belongs to the interval [-1, 1] and the integral on z is 1 whatever the position x within E. An important consequence is that the well known result of uniformity of ε and σ within E in the elastic case still holds in the ageing linear viscoelastic framework

∀ x ∈ E ε(x, t) = -P(t, •) • : p(•) , σ(x, t) = -C(t, •) • : P(•, ⋆) • : p(⋆) (34)
The strain tensor ε is uniform within E but depends on t and the polarization tensor P extending to ageing viscoelasticity the elastic Hill polarization tensor (but acting here as a Volterra kernel or in other words a function of two time variables) writes

P(t, t ′ ) = 1 4π ζ =1 Γ( A -1 • ζ, t, t ′ ) dS ζ (35)
or alternatively thanks to ( 26), ( 27) and (B.4)

P(t, t ′ ) = det A 4π ξ =1 ξ s ⊗ ξ • C • ξ -1 • (t, t ′ ) s ⊗ ξ A • ξ 3 dS ξ (36)
This general result recalls that P depends only on the shape and orientation of the ellipsoid and on the relaxation C of the reference medium. It is valid for any ellipsoid (even with three different radii) and possibly for an anisotropic tensor C. Nevertheless, as in elasticity, only few cases can lead to analytical expressions of P. In section 4, that of an isotropic relaxation tensor C is considered and a complete analytical expression of P is derived.

By analogy with the elastic case, the ageing viscoelastic Eshelby tensor is defined by S = P • : C. It is useful to relate the uniform strain tensor history within the ellipsoid subjected to a uniform eigenstrain history ε * (t)

∀ x ∈ E ε(x, t) = S(t, •) • : ε * (•) = P(t, •) • : C(•, ⋆) • : ε * (⋆) (37)
obtained thanks to the previous reasoning by introducing the polarization tensor p = -C

• : ε * in (34).

Extension to the inhomogeneity problem

The resolution of the inhomogeneity problem may also be inspired by that of the elastic case combined with the previous result as in [START_REF] Lavergne | Homogenization schemes for aging linear viscoelastic matrix-inclusion composite materials with elongated inclusions[END_REF]. Indeed, instead of (8b) let us consider now that the infinite domain is subjected to a remote boundary condition of the type

u(x, t) ∼ x →∞ E(t) • x ( 38
)
where E(t) is a given history of the macroscopic strain tensor. The relaxation tensor outside the ellipsoid E is still C but the material within E is now characterized by a different relaxation tensor denoted by C E . No polarization is considered here but the strain history solution ε(t) within E of the present problem can be built by considering the superposition in an infinite uniform medium of the uniform tensor E(t) and the response to a fictitious polarization of the form (C E -C)

• : ε within E assuming by anticipation that ε(t) remains uniform within E. Therefore, omitting the time variables for the sake of conciseness, it comes that

∀ x ∈ E ε = E -P • : (C E -C) • : ε (39)
and eventually

∀ x ∈ E ε = A • : E with A = HI + P • : (C E -C) -1 • (40)
where A, which depends on two time variables, is the ageing viscoelastic counterpart of the concentration tensor prevailing in the elastic framework [START_REF] Zaoui | Continuum Micromechanics: Survey[END_REF]. The uniformity of ε within E is then also ensured in the ageing viscoelastic framework: ε within E only depends on time.

It is worth emphasizing again the fact that the result ( 40) is valid whatever the shape and orientation of the ellipsoid and whatever the anisotropy of the matrix and inhomogeneity behaviours. Nevertheless the case of an anisotropic matrix behaviour may prevent from determining P (36) in a simple form. On the contrary, when the matrix is isotropic, the determination of P is easy to obtain as shown in the next section, allowing then the calculation of the concentration tensor (40) even if C E is anisotropic (see section 6.1 for practical calculations).

Analytical expression of the viscoelastic Hill tensor associated to an isotropic infinite medium

The isotropy of the reference medium means that the tensor kernel C(t, t ′ ) can be decomposed by means of a bulk relaxation kernel k(t, t ′ ) and a shear one µ(t, t ′ )

C(t, t ′ ) = 3 k(t, t ′ ) J + 2 µ(t, t ′ ) K (41)
where J = (1 ⊗ 1)/3 and K = I -J denote the orthogonal projectors onto the spaces of respectively spherical and deviatoric symmetric second-order tensors. Introducing the conjugated orthogonal projectors ξ ⊗ ξ and 1ξ ⊗ ξ for any unit vector ξ, the following relationships

ξ • J • ξ = 1 3 ξ ⊗ ξ and ξ • K • ξ = 2 3 ξ ⊗ ξ + 1 2 (1 -ξ ⊗ ξ) (42)
allow to write the isotropic acoustic tensor and its inverse in the form

K ξ = k + 4 3 µ ξ ⊗ ξ + µ(1 -ξ ⊗ ξ) (43
)

N ξ = K ξ -1 • = 3 (3 k + 4 µ) -1 • ξ ⊗ ξ + µ -1 • (1 -ξ ⊗ ξ) (44)
where K ξ , N ξ , k and µ implicitly depend on two time variables not recalled here for conciseness.

Introducing [START_REF] Schwartz | Théorie des distributions[END_REF] in [START_REF] Odegard | Physical aging of epoxy polymers and their composites[END_REF] yields

P = 3 (3 k + 4 µ) -1 • Q + µ -1 • (R -Q) ( 45
)
where P depends on two time variables through the kernels k and µ and Q and R are constant tensors defined by

Q = det A 4π ξ =1 ξ ⊗ ξ ⊗ ξ ⊗ ξ A • ξ 3 dS ξ (46) = 1 4π ζ =1 ξ ⊗ ξ ⊗ ξ ⊗ ξ dS ζ with ξ = A -1 • ζ A -1 • ζ (47)
and

R = det A 4π ξ =1 ξ s ⊗ 1 s ⊗ ξ A • ξ 3 dS ξ (48) = 1 4π ζ =1 ξ s ⊗ 1 s ⊗ ξ dS ζ with ξ = A -1 • ζ A -1 • ζ (49)
Remarkably the tensors Q and R depend only on the shape and orientation of the ellipsoid through A and not on the material properties. They are consequently rigorously identical as those that can be used to express the Hill polarization tensor in the elastic framework. In the latter case, replacing the relaxation functions k and µ by bulk and shear moduli and the inverse in the sense of Volterra by simple scalar inverse, (45) also provides the expression of P in the elastic case. For practical implementation, the components of Q and R are explicited in Appendix C.

In the case of a spherical inclusion, the tensors Q and R are simply derived by projection of the integrands of ( 47) and (49) onto J and K [START_REF] Bornert | Homogénéisation en mécanique des matériaux 1. Matériaux aléatoires élastiques milieux périodiques[END_REF]

Q = 1 3 J + 2 15 K and R = 1 3 I = 1 3 (J + K) (50)
It follows from ( 45) and ( 41) that the viscoelastic Hill and Eshelby tensors associated to a spherical inclusion respectively write

P = (3 k + 4 µ) -1 • • HJ + 3 5 (k + 2 µ) • µ -1 • K (51)
and

S = P • : C = 3 (3 k + 4 µ) -1 • • kJ + 2 5 (k + 2 µ)K (52)

Consistency with known results

The elastic and non-ageing viscoelastic cases

All the previous results obtained in the very general framework of ageing linear viscoelastic behaviour and ellipsoidal inclusion are consistent with their well established counterparts in elasticity or non-ageing viscoelasticity.

The consistency with the elastic behaviour is easily verified. Indeed, this behaviour corresponds to a relaxation tensor of the form C(t, t ′ ) = H(t, t ′ )C 0 where C 0 is the elastic stiffness tensor. As a consequence, [START_REF] Odegard | Physical aging of epoxy polymers and their composites[END_REF] becomes

P(t, t ′ ) = H(t -t ′ ) P 0 with P 0 = det A 4π ξ =1 ξ s ⊗ (ξ • C 0 • ξ) -1 s ⊗ ξ A • ξ 3 dS ξ ( 53
)
where P 0 is the well known elastic Hill tensor [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]. The other relationships [START_REF] Zheng | A new constitutive equation for the long-term creep of polymers based on physical aging[END_REF], [START_REF] Lavergne | Predicting the aging viscoelastic behavior of fiber-reinforced concrete by a multi-scale homogenization scheme[END_REF] or [START_REF] Gel'fand | Generalized functions -Vol I Properties and Operations[END_REF] give the corresponding elastic expressions by replacing the inverses and the tensor contractions in the sense of Volterra by classical inverses and tensor contractions and by multiplying the resulting fourth-order tensors by H.

The non-ageing viscoelastic case is almost as simple as the elastic one. As already mentioned in introduction, the correspondence principle can be applied to the set of equations characterizing the inclusion problem and the resolution in Laplace Transform Domain relies on the results already available in the elastic domain. But the objective here is to verify that the solution previously obtained by considering a general ageing behaviour is consistent with the particular non-ageing one. This verification is in fact straightforward by taking advantage of the properties of the Laplace-Carson transform (LCT) [START_REF] Salençon | Viscoélasticité pour le calcul des structures[END_REF]. Indeed the LCT of a Volterra product with a non-ageing kernel, i.e. depending on the pair (t, t ′ ) only through the difference tt ′ , becomes the corresponding product between the LCT of each terms. In addition the LCT of identities and inverses in the sense of Volterra are the corresponding (tensor or scalar) identities and inverses in the usual sense applied in the Laplace Transform Domain. Hence, as expected, all the classical elastic results, in particular the expressions of the Hill and Eshelby tensors and those of the concentration rules, are retrieved in the Laplace domain by applying the LCT to all the previous equations of this article.

The ageing linear viscoelastic spherical inhomogeneity with isotropic phases

This particular case has been solved in [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] using a direct displacement approach which eventually provides an explicit expression of the concentration tensor relating the strain tensor history within the sphere to the remote strain history. It is shown here that this tensor is exactly retrieved by [START_REF] Lavergne | Predicting the aging viscoelastic behavior of fiber-reinforced concrete by a multi-scale homogenization scheme[END_REF] thanks to (51) as a function of the relaxation bulk and shear moduli of the matrix (k and µ) and those of the spherical inhomogeneity (k E and µ E )

A(t, t ′ ) = A J (t, t ′ ) J + A K (t, t ′ ) K (54)
with

A J = (3 k E + 4 µ) -1 • • (3 k + 4 µ) (55)
A K = 5 9 k + 8 µ + 6(k + 2 µ) • µ -1 • • µ E -1 • • (3 k + 4 µ) (56)
On the one hand the expression of A J obtained in [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] is strictly identical to (55) but on the other hand the expression of A K provided in [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] requires some algebraic manipulations developed in Appendix D to lead to the more condensed relationship (56).

Note that the well known elastic expressions recalled in [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] are immediately retrieved by replacing the relaxation functions in (55) and (56) by the elastic bulk and shear moduli and the Volterra products and inverses by usual products and inverses within the set of real numbers.

It is worth emphasizing here that the scope of ( 40) is much wider in terms of inhomogeneity shape and behaviour since it can easily be applied to an anisotropic ellipsoidal inhomogeneity.

The case of a time-independent Poisson ratio

In [START_REF] Lavergne | Homogenization schemes for aging linear viscoelastic matrix-inclusion composite materials with elongated inclusions[END_REF], the ageing viscoelastic Eshelby inclusion problem is solved for materials featuring a relaxation tensor C of the form described by [START_REF] Hilton | The significance of (an)isotropic viscoelastic Poisson ratio stress and time dependencies[END_REF] C(t, t ′ ) = f (t, t ′ )C 0 (57)

where f (t, t ′ ) is a scalar function such that f (0, 0) = 1 and C 0 is the instantaneous stiffness tensor at time t = 0. When C 0 is isotropic, the relationship (57) implies that µ(t, t ′ )/k(t, t ′ ) does not depend on (t, t ′ ), which means a constant viscoelastic Poisson ratio. Nevertheless, the following conclusion remains valid if C 0 is anisotropic.

On the one hand, for such a material, it is shown in [START_REF] Lavergne | Homogenization schemes for aging linear viscoelastic matrix-inclusion composite materials with elongated inclusions[END_REF] that the ageing viscoelastic Eshelby tensor S is

S(t, t ′ ) = H(t -t ′ )S 0 (58)
where S 0 is the elastic Eshelby tensor associated to C 0 and an ellipsoid of the same shape.

On the other hand, introducing (57) in the comprehensive expression (36) yields a polarization tensor of the form f -1 • P 0 where the constant tensor P 0 writes as in [START_REF] Stehfest | Algorithm 368: Numerical inversion of {L}aplace transforms [D5[END_REF]. Applying then the formula [START_REF] Zheng | A new constitutive equation for the long-term creep of polymers based on physical aging[END_REF] providing S with (57) finally leads to the same expression as (58).

Practical implementation and validation

Numerical procedure to evaluate Volterra operators

The expressions of the Hill, Eshelby and strain concentration tensors established in this paper involve the Volterra operator and its inverse. In most practical cases, analytical evaluation is not tractable. Nevertheless, taking advantage of the algorithm proposed in [START_REF] Bazant | Numerical determination of long-range stress history from strain history in concrete[END_REF], it is possible to numerically evaluate the Volterra operator and its inverse. This approach, successfully used in [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] on the concentration problem on a spherical inhomogeneity, is briefly recalled here for scalar functions.

Considering that the value of all the involved functions of a single time variable is zero before t 0 and that the time interval of interest is [t 0 , t n ], a time sampling t 0 ≤ t 1 ≤ . . . ≤ t n is introduced and the time integral appearing in the Volterra operator ( 2) is discretized using the trapezoidal rule

Y i = i j=0 F(t i , t j ) + F(t i , t j-1 ) 2 (X j -X j-1 ) (∀ 0 ≤ i ≤ n) = F(t i , t i ) + F(t i , t i-1 ) 2 X i + i-1 j=0 F(t i , t j-1 ) -F(t i , t j+1 ) 2 X j (59)
in which the conventions t -1 = t 0 , X i = X(t i ), Y i = Y(t i ) for 0 ≤ i ≤ n and X -1 = 0 have been adopted.

Thus, it follows that the function of two different times F(t, t ′ ) playing the role of a Volterra kernel in ( 2) is now transformed into a lower triangular matrix [F] of rank n + 1 relating the vector

[Y] = [Y 0 , . . . , Y n ] to [X] = [X 0 , . . . , X n ] [Y] = [F][X] with F i j =                      0 if j > i F(t i , t i ) + F(t i , t i-1 ) /2 if j = i F(t i , t j-1 ) -F(t i , t j+1 ) /2 if j < i (60)
still with the convention t -1 = t 0 .

As already presented in [START_REF] Lavergne | Homogenization schemes for aging linear viscoelastic matrix-inclusion composite materials with elongated inclusions[END_REF], the previous methodology can be extended to tensor functions in order to build the numerical counterparts to expressions such as (3), (4) or even [START_REF] Lavergne | Predicting the aging viscoelastic behavior of fiber-reinforced concrete by a multi-scale homogenization scheme[END_REF] especially in the case of geometrical or physical anisotropy of the inhomogeneity. Indeed the expressions (59) are also valid for tensor variables X, Y and F provided that the scalar multiplications be replaced by adequate tensor contractions. In addition, tensor contractions can be seen as matrix-vector multiplications involving components in a particular vector frame. It is worth recalling here that the Voigt or the Kelvin-Mandel notations allow to represent strain or stress tensors as vectors of 6 components depending on time and to represent the relaxation, creep or even concentration tensor kernels as 6 × 6 matrices of functions of two time variables. As a consequence and taking into account the matrix or vector representations of tensors in a given frame, the matrix-vector product (60) remains also valid but should be considered in the sense of a block matrix operation since the i th component of [X] (resp. [Y]) contains itself a vector of X (resp. Y) at time t i and each component of the matrix [F] represents a block matrix computed with the relationships (60) applied on the tensor kernel F. For example, (3) becomes a matrix-vector product in a space of dimension 6(n + 1) in which the matrix is lower triangular and contains (n + 1) × (n + 1) blocks of 6 × 6 submatrices. Alternatively the block matrix operations can also be built by means of a 6 × 6 matrix of (n + 1) × (n + 1) block matrices if the components corresponding to space and time are reordered. The inverse of a tensor kernel in the sense of Volterra simply becomes a matrix inverse in this numerical procedure.

Thus, practical implementation of the results derived in this paper is straightforward and only requires

• to time-discretize relaxation functions into matrices using (60),

• to implement expressions of the Hill and concentration tensors, replacing the Volterra operator and its inverse by matrix multiplication and inversion, using any software package implementing numerical matrix algebra.

Validation of the strain concentration tensor with respect to FEM computations

Apart form the well known elastic or non ageing cases (section 5.1), and the spherical inhomogeneity with ageing case (section 5.2), reference results in more general situations can be obtained resorting to numerical computations.

The concentration problem is solved here in ageing linear viscoelasticity, for an oblate spheroid of axis e 3 and aspect ratio 0.2, using the Code_Aster [START_REF] Edf | Code_Aster: analyse des structures et thermomécanique pour des études et recherches[END_REF] finite element code. A spheroid of radii a = b = 5, c = 1 is embedded into a box Ω which half lengthes are χ Ω times the spheroid radii. Uniform strain type boundary conditions following a macroscopic strain step are imposed, and the evolution of the average strain in the inhomogeneity is sought. The six independent components of the concentration tensor are obtained by considering six elementary orientations of the macroscopic strain. Taking advantage of the axisymmetric nature of the geometry, some computations can be performed in 2D axisymmetric conditions. Besides, taking advantage of the three symmetry planes, only one eight of the domain needs to be meshed (one half in the 2D axisymmetric case). Thus, 2D and 3D quadratic meshes have been prepared, using the Salome [START_REF] Cea | Salome: the open source integration platform for numerical simulation[END_REF] platform, ensuring similar element sizes in critical zones of the computation domain as detailed hereafter. The numerical parameters affecting accuracy are the relative size χ Ω of the computing domain (which is not infinite, in contrary to the Eshelby problem), and mesh refinement.

A preliminary convergence analysis is performed, considering (see details in table 1) • three relative domain sizes: χ Ω = 10, 20 and 40,

• three sets of meshes: coarse, medium and fine. The meshes are refined around the inhomogeneity and especially around the highest curvature zone of the interface (figure 1).

This analysis is restricted to elastic behaviours (see elastic characteristics of inhomogeneity and matrix in table 2), to use the classical Eshelby result as a reference. The relative errors on components 1133 and 3333 of the strain concentration tensor are plotted on figure 2. The relative domain size χ Ω is found to have a much larger influence on accuracy than mesh refinement. Similar trends are found on the other components of the strain concentration tensor. Thus, the ageing linear viscoelastic numerical computations, subsequently used as references, are performed on the largest domain (χ Ω = 40) and with the medium mesh, as a compromise to keep computing times reasonable.

As it is the only ageing linear viscoelastic behaviour available in Code_Aster, the Granger model [START_REF] Granger | Comportement différé du béton dans les enceintes de centrales nucléaires[END_REF] is adopted for both matrix and inhomogeneity. A constant Poisson ratio is assumed:

L(t, t ′ ) = L E (t, t ′ ) [(1 -2ν)J + (1 + ν)K] (61)
χ Ω = 10 The uniaxial compliance is here defined with two Kelvin chains, and an arbitrary ageing function f a (t ′ ):

χ Ω = 20 χ Ω =
L E (t, t ′ ) = 1 E + f a (t ′ ) 2 i=1 s i 1 -e -t-t ′ τ i with f a (t ′ ) = e -t ′ τa 2 (62)
Arbitrary material parameters are chosen with the inhomogeneity stiffer than the matrix (see table 2). Bulk and shear compliance functions of both matrix and inhomogeneity phases are plotted in figure 3. The six independent components of the strain concentration tensor are plotted in figures 4 and 5. The latter are functions of (t, t ′ ) and correspond to the uniform strain arising in the inhomogeneity at t for a unit strain step occurring at t ′ at the domain boundary. The semi-analytical result put in evidence in this paper agrees very well with reference FEM computations for all the six components and is able to capture all the details of non-monotonic evolutions.

ν E s 1 s 2 τ 1 τ 2 τ a matrix 0.
Regarding convergence analysis of FEM results, trends similar to the elastic case are found: Figure 2: Relative error on selected components of the strain concentration tensor, with elastic behaviours, between FEM computations and the analytical reference result from Eshelby, using 3 different mesh refinements and 3 different domain sizes.

• the relative domain size χ Ω has a large impact on accuracy (see details in the right part of figures 4 and 5), especially on components 1122, 1133 and 3311,

• the mesh refinement has a much lesser impact on accuracy (figure 6), suggesting that ageing linear viscoelastic computations on the fine mesh would not provide significant improvements.

Performance of the present numerical method based on quadrature in time domain compared to the correspondence principle in the non-ageing case

In order to assess the performance of the numerical procedure based on the trapezoidal rule ( [START_REF] Bazant | Numerical determination of long-range stress history from strain history in concrete[END_REF], [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF]) presented in this paper to calculate the viscoelastic concentration tensor of the Eshelby problem, a comparison with an alternative technique based on the correspondence principle (Laplace-Carson) is proposed for the non-ageing case. To this aim, the hypotheses of the concentration problem presented in the previous section devoted to FEM computations are considered again except that the ageing functions in (62) are replaced by 1.

On the one hand, the Laplace-Carson transforms of the elementary behaviours modelled as Kelvin chains have very simple analytical forms. The concentration problem becomes elastic in the Laplace domain and the application of the Gaver-Stehfest algorithm (see [START_REF] Stehfest | Algorithm 368: Numerical inversion of {L}aplace transforms [D5[END_REF]) allows to come back to the time domain and to express the six independent components of the concentration tensor at any time. If N denotes the number of times at which the latter is calculated, it comes that the complexity of this technique is O(N) since the number of calls of the elastic problem with the Gaver-Stehfest algorithm at each time is limited.

On the other hand, the procedure developped in the present paper allows to compute the whole history (with respect to the chosen discrete set of times) of the concentration tensor in only one calculation. As shown in the paper the latter requires to inverse lower triangular block matrices of rank 6N which is the most time consuming step of the procedure. It follows that the expected complexity is at best O(N 2 ).

Both techniques have been implemented in Python 2.7 codes and performed on only one core of an Intel i7-3840QM CPU @ 2.80GHz with 8Gb of RAM. The Gaver-Stehfest algorithm has been used with parameter M = 8 (summation on 2M = 16, see [START_REF] Stehfest | Algorithm 368: Numerical inversion of {L}aplace transforms [D5[END_REF]). It must be precised that both methods perfectly coincide, even for the lowest number of times (N = 150 between t = 0 and t = 10). The evolutions of the concentration tensor are not presented because they are identical to the ones appearing in right figures 4 and 5. The computing time is plotted against the size of time sampling in figure 7. The results are consistent with the expected complexities evoked in the previous paragraphs. In particular, it appears that the method based on the use of the correspondence principle is more efficient. Nevertheless, this method is designed only for non-ageing behaviours whereas the approach proposed in this paper is more general since it can address the ageing case. Moreover, as recalled in introduction, it should be kept in mind that the algorithms allowing to calculate the inverse of the Laplace transform suffer from several shortcomings and should be carefully adapted to the function to inverse.

Concluding remarks

In this paper, the complete solution to the Eshelby problem of an ellipsoidal inclusion submitted to a uniform polarization history and embedded in an infinite matrix has been derived in the general framework of ageing linear viscoelasticity. In particular it has been proven that the strain state and subsequently the stress state are uniform within the ellipsoid and depend only on time. The structure of the solution is similar to that of the elastic case provided that tensor contractions be replaced by Volterra products: the strain tensor is related to the polarization tensor thanks to a polarization fourth-order Volterra kernel generalizing the elastic Hill polarization tensor. Furthermore, this polarization kernel writes as an integral over the unit sphere of a function depending on the shape of the ellipsoid and of the viscoelastic counterpart of the elastic acoustic tensor of the matrix. In addition, generalizations to ageing linear viscoelasticity of the Eshelby tensor and of the concentration tensor in the case of an ellipsoidal inhomogeneity are provided. Complete analytical expressions of the Hill polarization kernel are determined in the case of an isotropic matrix. The theoretical calculations of the paper are then compared to several published results in order to verify their consistency in particular cases (elasticity, non-ageing viscoelasticity, spherical inclusion). Finally, after recalling an efficient numerical procedure to evaluate Volterra operators, solutions in the case of an ageing spheroidal inhomogeneity in an ageing matrix are very satisfactorily compared to finite element simulations.

The generalization of Eshelby's results to ageing linear viscoelasticity as presented in this paper combined with the numerical procedure to evaluate Volterra operators opens a new way to tackle upscaling problems of random media made up with general ageing linear viscoelastic constituents of anisotropic shape. Indeed it becomes rather straightforward to implement classical schemes such as the Mori-Tanaka (for an isotropic matrix) or the self-consistent ones (by means of an iterative fixed-point procedure for a macroscopically isotropic material so as to keep an isotropic fictitious matrix at each iteration) in order to estimate the macroscopic properties of heterogeneous materials with ageing linear viscoelastic random microstructure. As in the elastic framework, the particular case of cracks deserves a special attention: it is currently under investigation and is the topic of a paper in preparation. which ends the proof by inserting (B.3) in (B.2). It is worth precising that the following relationships are similarly obtained

ζ = A • ξ A • ξ and dS ζ = det A A • ξ 3 dS ξ (B.4)
It may be noticed that this change of variable has been considered in [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] and the relationship between surface elements has been obtained by a reasoning on components (see (16.10), (17.09) and (17.10) in [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]).

Appendix C. Tensors Q and R and Newtonian potential

The dependence of Q and R on A implies that those tensors are orthotropic in the frame (e 1 , e 2 , e 3 ) defining the ellipsoid [START_REF] Mandel | Cours de Mécanique des Milieux Continus[END_REF]. The non-zero components of Q and R in this frame are given by (the summation over repeated indices is not applied in what follows)

Q iiii = 3(I i -ρ 2 i I ii ) 2 ∀ i ∈ {1, 2, 3} (C.1) Q ii j j = Q i ji j = Q i j ji = I j -ρ 2 i I i j 2 = I i -ρ 2 j I i j 2 ∀ i j ∈ {1, 2, 3} (C.2) and R iiii = I i ∀ i ∈ {1, 2, 3} (C.3) R i ji j = R i j ji = I i + I j 4 ∀ i j ∈ {1, 2, 3} (C.4)
The coefficients I i and I i j used in (C.1), (C.2), (C.3) and (C.4) are adapted from those provided in [START_REF] Kellogg | Potential theory[END_REF] and [1] (i.e. differ by a factor of 4π/3 for I i j with i j and by 4π for the others). Assuming that ρ 1 ≥ ρ 2 ≥ ρ 3 without loss of generality and renaming ρ 1 = a, ρ 2 = b and ρ 3 = c when useful for a better readability of the following formulas, different cases are considered The expression of A K provided in [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] writes with the notations of the present paper

• if a > b > c I 1 = a b c (a 2 -b 2 ) √ a 2 -c 2 (F -E) (C.5) I 3 = a b c (b 2 -c 2 ) √ a 2 -c 2        b √ a 2 -c 2 a c -E        (C.6) I 2 = 1 -I 1 -I 3 (C.7) I i j = I j -I i ρ 2 i -ρ 2 j ∀ i j ∈ {1, 2, 3} (C.8) I ii = 1 3         1 ρ 2 i - j i I i j         ∀ i ∈ {1, 2, 3} (
A K = H + 2(2H + 3D) • 2µ E • (2H + 3D) + µ • (6H -D) -1 • • (µ -µ E ) (D.1)
with

D = 2 3 (k + µ) -1 • • µ (D.2) Replacing 2H + 3D = 2(k + µ) -1 • • (k + 2µ) (D.
A K = H + 6 6µ E + M -1 • • (µ -µ E ) (D.5) with M = µ • (k + µ) -1 • • (9k + 8µ) • (k + 2µ) -1 • • (k + µ) (D.6)
The terms appearing in M in (D.6) do not commute but M advantageously rewrites with respect to ρ = k

• µ -1 • M = (ρ + H) -1 • • (9ρ + 8H) • (ρ + 2H) -1 • • (ρ + H) • µ (D.7)
Considering now that ρ, H, linear combinations of the latter as well as their inverses commute, (D.7) becomes

M = (ρ + 2H) -1 • • (9ρ + 8H) • µ = µ • (k + 2µ) -1 • • (9k + 8µ) (D.8)
The last expression of M in (D.8) is then used in (D.5) to end the proof

A K = H + 6 6µ E + µ • (k + 2µ) -1 • • (9k + 8µ) -1 •
• (µµ E ) (D.9)

= H + 6 9k + 8µ + 6(k + 2µ) • µ -1 • • µ E -1 • • (k + 2µ) • µ -1 • • (µ -µ E ) (D.10) = 5 9 k + 8 µ + 6(k + 2 µ) • µ -1 • • µ E -

Figure 1 :

 1 Figure 1: 2D axisymmetric (left) and 3D (right) medium meshes for χ Ω = 10: global view (top) and detailed view on inhomogeneity (bottom).

Figure 3 :

 3 Figure 3: Bulk (plain lines) and shear (dashed lines) compliance functions of matrix and inhomogeneity phases, for t ′ = 0, 10, 20, 30, 40, 50.

C. 9 )I 2 = I 3 = a a √ a 2 -2 a 2 -c 2 I 1 = 1 - 2 I 3 (C. 12 ) 2 -ρ 2 i••I 1 = I 2 = I 3 2

 92222112312221232 where F = F (θ, κ) and E = E(θ, κ) are respectively the elliptic integrals of the first and second kinds of amplitude and parameterθ = arcsin 1 -c 2 a 2 ; κ = a 2b 2 a 2c 2 (C.10) • if a > b = c (prolate spheroid) c 2c 2 arccosh (a/c) I 1i = I i1 = I i -I 1 a if a = b > c (oblatespheroid)I 1 = I 2 = c a 2 arccos (c/a)c √ a 2c 2 2 a 2c 2 3/2 (C.16)I 3 = 1 -2 I 1 (C.17)I 3i = I i3 = I 3 -I i ρ 2 ic 2 ∀ i ∈ {1, 2} (C.18) if a = b = c (sphere) ∀ i, j ∈ {1, 2, 3} (C.22)Appendix D. Proof of equality between (56) and the expression of A K provided in[START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF]

Figure 5 :Figure 6 :

 56 Figure5: Components of strain concentration tensor: Eshelby and FEM (medium mesh) results, for t ′ = 0, 10, 20, 30, 40, 50. Right: zoom on small t, for t ′ = 0; FEM results for the medium mesh and the 3 investigated domain sizes are plotted with dotted/dashed lines.

Figure 7 :

 7 Figure 7: Comparison of computation time required to compute the viscoelastic concentration tensor with the approach proposed in this paper (quadrature in time domain), or with the correspondence principle (Laplace-Carson) resorting to numerical Laplace transform inversion using the Gaver-Stehfest algorithm [53].

Table 1 :

 1 Number of nodes and elements in the considered quadratic meshes.

	40

Table 2 :

 2 Material properties of matrix and inhomogeneity phases (arbitrary units).

  Components of strain concentration tensor: Eshelby and FEM (medium mesh) results, for t ′ = 0, 10, 20, 30, 40, 50. Right: zoom on small t, for t ′ = 0; FEM results for the medium mesh and the 3 investigated domain sizes are plotted with dotted/dashed lines.

						A 1111 (t, t ′ ) A 3311 (t, t ′ )					A 1111 (t, t ′ ) zoom A 3311 (t, t ′ ) zoom
		0.8 0.04			Eshelby χ Ω = 40 Eshelby χ Ω = 40				0.04 0.8	Eshelby Eshelby χ Ω = 40 χ Ω = 40
											χ Ω = 20 χ Ω = 20
	A 1111 A 3311	0.4 0.6 0.02						A 1111 A 3311	0.7 0.02	χ Ω = 10 χ Ω = 10
		0.2	0							0.6 0.00
	0 0.01 0 0 0.02 0.03 0 0.01 0.02 0.03 0.04 0 0.1 0.2 0.3 0.4 0.5 0.6 0 0 0.1 0.2 0.3 0 Figure 4: -0.02 0 A 1122 0 A 1133 A 3333 A 2323	0	20 20 20 20 20 20	40 40 A 1122 (t, t ′ ) 60 t 60 t 40 60 t A 1133 (t, t ′ ) 40 60 t 40 60 t A 3333 (t, t ′ ) 40 60 t A 2323 (t, t ′ )	80 80 80 80 80 80	100 100 100 100 100 100	A 1122 A 1133 A 3333 A 2323	• 0.5 0.015 • (3 k + 4 µ) 0 0.020 0.025 0 0.030 0.035 0.040 0 0 0.2 0.3 0.4 0.5 0.6 0 0.10 0.15 0.20 0.25 0	1 1 1 A 1122 (t, t ′ ) zoom 2 t 2 t 1 2 t A 1133 (t, t ′ ) zoom 2 t 1 2 t A 3333 (t, t ′ ) zoom 1 2 t A 2323 (t, t ′ ) zoom	(D.11) 3 3 3 3 3 3

H(t) = 1 if t ≥ 0 and H(t) = 0 if t < 0

δ i j = 1 if i = j and δ i j = 0 if i j

χ E (x) = 1 if x ∈ E and 0 otherwise

Acknowledgements

The first author (J.-F. Barthélémy) expresses his gratitude to Dr. S. Cariou for fruitful discussions and her encouragement during the course of this work.

The second author (A. Giraud) gratefully acknowledges the support from TAMER European project (Trans-Atlantic Micromechanics Evolving Research Materials containing inhomogeneities of diverse physical properties, shapes and orientations), FP7 Project TAMER IRSES-GA-2013-610547.

Appendix A. Proof of the plane-wave expansion (11)

The starting point of the proof is the relationship

obtained by introducing spherical coordinates (θ, ϕ) of ξ such that z = ξ • x = x cos θ. The change of variable (θ, ϕ) → (z, ϕ) allows to write the right hand side of (A.1) under the form

which leads to the relationship (A.1) thanks to the property of the Dirac distribution δ ( [START_REF] Schwartz | Théorie des distributions[END_REF], [START_REF] Gel'fand | Generalized functions -Vol I Properties and Operations[END_REF]).

Finally taking the Laplacian with respect to the variable x to each side of (A.1) and recalling that ξ = 1 and ∆ (1/ x ) = -4πδ (three-dimensional Poisson formula [START_REF] Schwartz | Théorie des distributions[END_REF]), the relationship (11) is retrieved.

Appendix B. Change of variable over the unit sphere

This appendix aims at studying the change of variable from the unit sphere to itself ζ → ξ = A -1 •ζ/ A -1 • ζ with A a symmetric second-order tensor and more specifically at proving the relationship (28) between the infinitesimal surface elements dS ξ and dS ζ .

By differentation, the infinitesimal vectors are related by

This corresponds to the application of a transformation gradient A -1 / A -1 • ζ followed by an orthogonal projection onto the plane of normal ξ. The infinitesimal surface elements are then related by

The definition of the unit vectors ξ and ζ implies