F Lavergne

K Sab

J Sanahuja

M Bornert

C Toulemonde

Homogenization schemes for aging linear viscoelastic matrix-inclusion composite materials with elongated inclusions

Keywords: aging, creep, homogenization, ellipsoid, concrete

An extension of the Mori-Tanaka and Ponte Castañeda-Willis homogenization schemes for linear elastic matrix-inclusion composites with ellipsoidal inclusions to aging linear viscoelastic composites is proposed. To do so, the method of Sanahuja [81] dedicated to spherical inclusions is generalized to ellipsoidal inclusions under the assumption of time-independent Poisson's ratio. The obtained time-dependent strains are successfully compared to those predicted by an existing method dedicated to time-shift aging linear viscoelasticity showing the consistency of the proposed approach. Moreover, full 3D numerical simulations on complex matrix-inclusion microstructures show that the proposed scheme accurately estimates their overall time-dependent strains. Finally, it is shown that an aspect ratio of aggregates in the range 0.3 to 3 has no significant influence on the time-dependent strains of composites with per-phase constitutive relations representative of a real concrete.

Introduction

Eshelby's solution [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] of an ellipsoidal inclusion in an elastic material has been used in various ways to upscale the behavior of composite materials. The Mori-Tanaka [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF][START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF] scheme and the Ponte Castañeda-Willis scheme [START_REF] Ponte Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF] have been designed to retrieve the elastic behavior of composites featuring spherical [START_REF] Weng | Some elastic properties of reinforced solids, with special reference to isotropic ones containing spherical inclusions[END_REF] or elongated inclusions [START_REF] Tandon | The effect of aspect ratio of inclusions on the elastic properties of unidirectionally aligned composites[END_REF]. Such mean-field homogenization schemes based on Eshelby's solution have been coupled to the correspondence principle [START_REF] Lee | Stress analysis for linear viscoelastic materials[END_REF][START_REF] Mandel | Cours de mécanique des milieux continus[END_REF] to estimate the time-dependent strains of nonaging viscoelastic materials [START_REF] Wang | Influence of inclusion shape on the overall viscoelastic behavior of composites[END_REF][START_REF] Brinson | Comparison of micromechanics methods for effective properties of multiphase viscoelastic composites[END_REF][START_REF] Lévesque | Numerical inversion of the Laplace-Carson transform applied to homogenization of randomly reinforced linear viscoelastic media[END_REF]: the Laplace-Carson transform turns the nonaging problem into a set of formal elastic problems in complex space. For composite materials made of elastic inclusions and a matrix modeled by a time shift method, such as many 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], the Laplace-Carson transform may still be applied in the equivalent time space [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]. Yet, inverting the Laplace-Carson transform is still a compromise between accuracy and stability since this operation is ill-conditioned. This is one of the reasons why modern homogenization methods operate in the time domain [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A timeintegration approach[END_REF][START_REF] Tran | A simple computational homogenization method for structures made of linear heterogeneous viscoelastic materials[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].

Regarding aging viscoelastic materials, a closed-form solution has been proposed by Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] to handle the case of spherical inclusions in an aging linear viscoelastic matrix. Moreover, a reliable numerical procedure has been proposed to efficiently estimate the time-dependent strains. This procedure does not require inverting the Laplace-Carson transform and is able to handle any isotropic compliance.

This paper is devoted to validating and extending Sanahuja's method to ellipsoidal inclusions. Yet, the extension is limited to isotropic aging viscoelastic matrices featuring time-independent Poisson's ratio in the sense of Hilton [START_REF] Hilton | The significance of (an)isotropic viscoelastic Poisson ratio stress and time dependencies[END_REF].

This extension may be valuable to study cementitious materials. Indeed, modern formulations of concrete may include aggregates [START_REF] De Larrard | Structures granulaires et formulation des bétons -Concrete Mixture -Proportionning -A scientific approach[END_REF], steel fibers, expanded polystyrene particles [START_REF] Babu | Behaviour of lightweight expanded polystyrene concrete containing silica fume[END_REF][START_REF] Roy | Taking into account the inclusions' size in lightweight concrete compressive strength prediction[END_REF] or wood shavings [START_REF] Bederina | Effect of the addition of wood shavings on thermal conductivity of sand concretes: Experimental study and modelling[END_REF] as inclusions and such inclusions can change the viscoelastic properties of the material [START_REF] Chern | Compressive creep and shrinkage of steel fibre reinforced concrete[END_REF].

A recent numerical study has shown that the size distribution and the shape of aggregates have little effect on the time-dependent strains of concretes made with non-aging cementitious matrices [START_REF] Lavergne | Investigation of the effect of aggregates' morphology on concrete creep properties by numerical simulations[END_REF]. Full 3D numerical simulations and semi-analytical homogenization schemes delivered similar estimates of the time-dependent strains. Yet, this study was limited to polyhedral aggregates with an aspect ratio close to 1 (the aggregates were neither flat nor elongated). Consequently, there is a question left: does the aspect ratio of the aggregates affect the time-dependent strains of concrete made with aging cementitious matrices? To answer this question, the following steps have been carried out:

In section 1, the method of Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] is extended to ellipsoidal inclusions. Eshelby's solution for an isotropic aging viscoelastic matrix featuring a time-independent Poisson's ratio is presented. The strain within the inclusion is still uniform and a time-dependent Eshelby's tensor may be defined. Then, the closed-form of the time-dependent localization tensor is derived using the Volterra operator. Finally, the Mori-Tanaka estimate of the overall viscoelastic behavior is obtained. The more sophisticated Ponte Castañeda-Willis linear estimate, which accounts separately for inclusion shape effects and effects of the spatial distribution of inclusion centers, is formally similarly extended to composites with elastic inclusions embedded in an aging viscoelastic matrix.

In section 2, the numerical procedure described in [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] is used to evaluate the proposed homogenization schemes. The resulting estimates of the overall timedependent strains are compared to existing ones for a fiber-reinforced polymer with a time-shift aging viscoelastic matrix. Then, a complex microstructure featuring 60% of polyhedral elastic aggregates embedded in an aging viscoelastic cementitious matrix is considered. The estimates of the time-dependent strains as evaluated by the method of Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] and by full 3D numerical computations are first compared for spherical inclusions. Finally, the proposed extension to ellipsoidal inclusions of the Mori-Tanaka and Ponte Castañeda-Willis schemes are used to estimate the overall response of concrete-like materials. In particular, it will be shown that the aspect ratio of aggregates used in concrete does not significantly affect their viscoelastic behavior.

. Aging viscoelasticity

The stress tensor σ(t) in a viscoelastic material depends on the history of strain tensor ε(t). If the constitutive law is linear, the Boltzmann superposition principle states that the material properties are defined by a relaxation function (fourth order tensor), g(t, t ), such that:

σ(t) = t -∞ g(t, t )dε(t )
where the integral is a Stieltjes integral. Similarly, the compliance function (fourth order tensor), t(t, t ) is such that:

ε(t) = t -∞ t(t, t )dσ(t )
If the elapsed time since loading is the only relevant parameter, the material is non-aging:

t(t, t ) = ¨(t -t )
However, the assumption of non aging is not made in the following derivations. If the viscoelastic behavior is isotropic, a spherical relaxation function K(t, t ) and a deviatoric relaxation function G(t, t ) are defined, such that:

p(t) = t -∞ 3K(t, t )de(t ) σ d (t) = t -∞ 2G(t, t )dε d (t )
where e(t) = tr(ε(t))/3, p(t) = tr(σ(t))/3. The tensors σ d (t) and ε d (t) are respectively the deviatoric parts of σ(t) and ε(t):

σ d ij (t) = σ ij (t) -p(t)δ ij ε d ij (t) = ε ij (t) -e(t)δ ij
where δ ij is the Kronecker symbol.

The homogenization method of Sanahuja

The homogenization method of Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] operates in the time domain to deal with a composite featuring isotropic aging viscoelastic phases. A spherical inclusion featuring a linear viscoelastic isotropic behavior (K i (t, t ), G i (t, t )) is embedded in an infinite matrix featuring a linear viscoelastic isotropic behavior (K m (t, t ), G m (t, t )) and a strain history E(t) is applied far from the inclusion.

As usually for linear elasticity, the solution of this problem provides the exact solution of the localization problem of a composite made of such inclusions embedded in the matrix, in the so-called "dilute limit", i.e. for volume fractions of inclusions sufficiently low so that mechanical interactions between inclusions can be neglected. For isotropic composites, this localization solution is fully determined when purely spherical or deviatoric overall strains histories E(t) are considered. Sanahuja proved that the strain history in the inclusion is uniform, with value ε(t). In addition, a localization tensor e(t, t ) is defined such that:

ε(t) = t -∞ e(t, τ ) : dE(τ )
To ease the computations, the Volterra operator has been introduced. It reads:

f • g(t, t ) = t -∞ f (t, τ )d τ g(τ, t )
for any scalar functions f and g. The identity element of the Volterra operator H is defined from the Heaviside function:

(t, t ) → H(t -t ) = 1 if(t > t ) 0 if(t < t )
The value of the H function at t = t does not need to be specified. The inverse of f in the sense of the Volterra operator is denoted as

f -1 so that f -1 • f = H.
The spherical part A k (t, t ) of e(t, t ) reads:

A k = (3K i + 4G m ) -1 • (3K m + 4G m )
The closed-form expression of the deviatoric part A g (t, t ) has been computed as well:

A g = H + 2(2H + 3D m ) • (2G i • (2H + 3D m ) + G m • (6H -D m )) -1 (G m -G i )
where

D m = (K m + G m ) -1 • 2 3
G m These expressions provide functionnal relations between time histories of spherical and deviatoric strain prescribed far away from the inclusion and the induced uniform strain history in the latter. They can be used to extend, at least from a formal point of view, to aging viscoleasticity any linear elastic homogenization scheme based on Eshelby's solution, by simply substituting the classical tensor double contraction operations by Volterra operations. The resulting expressions will involve multiple time convolutions as well as Volterra inversions, and will thus be rather involved.

To face this difficulty, numerical procedures to compute the Volterra operator and its inverse proposed in [START_REF] Bažant | Numerical determination of long-range stress history from strain history in concrete[END_REF] are used to turn these formula into a practical tool [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF]. The estimate of the effective behavior of a concrete will be compared to results of 3D numerical simulations in section 2.3.1. It is to be mentioned that the method of Sanahuja produces an estimate of the effective behavior at once: the output is a matrix representing the global linear viscoelastic behavior of the composite material over an initially specified period of time. A single run of the method of Sanahuja handles all loading directions at once while two full 3D numerical computations are needed to simulate hydrostatic and shear creep tests. 4

In its initial formulation, the method of Sanahuja is limited to homogenization models based on the Eshelby solution of an isotropic spherical inclusion embedded in an isotropic matrix. Since solving the Eshelby inhomogeneity problem for ellipsoidal inclusions is likely to be a tedious and complex task because in particular of the multiplicity of the loading cases to be considered to address the induced anisotropy, another viewpoint is introduced in the next section. Eshelby's inhomogeneity problem is solved for ellipsoidal inclusions provided that the linear viscoelastic matrix features a time-independent Poisson's ratio which means that the time dependence of the relaxation function is purely scalar, of the form g(t, t ) = f (t, t )g(0, 0), as illustrated in figure 1

The hypothesis of a time-independent Poisson's ratio

It is well known that the contraction ratio of a viscoelastic materials is a function of time or frequency [START_REF] Lakes | On Poisson's ratio in linearly viscoelastic solids[END_REF][START_REF] Lee | Determination of viscoelastic Poisson's ratio and creep compliance from the indirect tension test[END_REF]. The ratio of the lateral strain to the axial strain during a tensile creep test is called the creep Poisson's ratio as in [START_REF] Kim | Experimental studies on creep of sealed concrete under multiaxial stresses[END_REF], the effective time-dependent Poisson's ratio ν(t, t ) being the ratio of the lateral strain ε lat to the axial strain ε during a tensile stress relaxation test [START_REF] Tschoegl | The Phenomenological Theory of Linear Viscoelastic Behavior, an introduction[END_REF][START_REF] Wineman | Mechanical Response of Polymers, an introduction[END_REF][START_REF] Tschoegl | Poisson's ratio in linear viscoelasticity -a critical review[END_REF][START_REF] Kim | Experimental studies on creep of sealed concrete under multiaxial stresses[END_REF][START_REF] Lakes | On Poisson's ratio in linearly viscoelastic solids[END_REF][START_REF] Lee | Determination of viscoelastic Poisson's ratio and creep compliance from the indirect tension test[END_REF]:

ν(t, t ) = - ε lat ε
The time-dependent Poisson's ratio is consistent with the Poisson's ratio defined in the Laplace-Carson space [START_REF] Tschoegl | Poisson's ratio in linear viscoelasticity -a critical review[END_REF][START_REF] Lakes | On Poisson's ratio in linearly viscoelastic solids[END_REF][START_REF] Li | Computing the time evolution of the apparent viscoelastic/viscoplastic Poisson's ratio of hydrating cement paste[END_REF] or in the frequency domain [START_REF] Read | The Determination of Dynamic Properties of Polymers and Composites[END_REF][START_REF] Pritz | Measurement methods of complex Poisson's ratio of viscoelastic materials[END_REF][START_REF] Tschoegl | Poisson's ratio in linear viscoelasticity -a critical review[END_REF][START_REF] Boiko | Measurement method of complex viscoelastic material properties[END_REF] for non-aging materials and it represents a valid constitutive property [START_REF] Tschoegl | Poisson's ratio in linear viscoelasticity -a critical review[END_REF][START_REF] Li | Computing the time evolution of the apparent viscoelastic/viscoplastic Poisson's ratio of hydrating cement paste[END_REF]. The experimental measurement of a creep Poisson's ratio demands highly accurate measurements be made on the same specimen, at the same time, and under the same conditions of the experimental environment [START_REF] Yee | Dynamic bulk and shear relaxation in glassy polymers. i. experimental techniques and results on PMMA[END_REF][START_REF] Tschoegl | Poisson's ratio in linear viscoelasticity -a critical review[END_REF].

Regarding cementitous materials, the creep Poisson's ratio of concretes [START_REF] Parrott | Lateral strains in hardened cement paste under short-and long-term loading[END_REF][START_REF] Kogan | Creep of concrete under multiaxial compression[END_REF][START_REF] Kim | Experimental studies on creep of sealed concrete under multiaxial stresses[END_REF]] has been measured. The time-dependent Poisson's ratio [START_REF] Kim | Experimental studies on creep of sealed concrete under multiaxial stresses[END_REF] has been estimated and was found to be almost uniform at constant humidity (approximately ν = 0.18) [START_REF] Bažant | Mathematical modeling of creep and shrinkage of concrete[END_REF][START_REF] Kim | Experimental studies on creep of sealed concrete under multiaxial stresses[END_REF]. This fact has been used to model concrete [START_REF] Granger | Comportement différé du béton dans les enceintes de centrales nucléaires : analyse et modélisation[END_REF][START_REF] Bažant | Creep and shrinkage prediction model for analysis and design of concrete structures-model B3[END_REF][START_REF] Bažant | Justification and refinements of model B3 for concrete creep and shrinkage 1. statistics and sensitivity[END_REF]. Grasley and Lange found that the Poisson's ratio of cement paste was increasing for some materials or almost time-independent for others [START_REF] Grasley | Constitutive modeling of the aging viscoelastic properties of Portland cement paste[END_REF][START_REF] Li | Computing the time evolution of the apparent viscoelastic/viscoplastic Poisson's ratio of hydrating cement paste[END_REF]. The accuracy of experimental measurements of the Poisson's ratio may be undermined by additional strains such as autogeneous shrinkage in cementitious materials, which must be accounted for [START_REF] Grasley | Constitutive modeling of the aging viscoelastic properties of Portland cement paste[END_REF]. The viscoelastic Poisson's ratio of a hydrating cement paste has been investigated by full 3D numerical simulations [START_REF] Li | Computing the time evolution of the apparent viscoelastic/viscoplastic Poisson's ratio of hydrating cement paste[END_REF][START_REF] Li | Modeling the apparent and intrinsic viscoelastic relaxation of hydrating cement paste[END_REF] in which the cement paste is described as a composite material. The dissolution of cement grains results in a predicted decreasing viscoelastic Poisson's ratio at early ages, but later-age cement paste is capable of exhibiting increasing, decreasing, or uniform evolution of the Poisson's ratio with time [START_REF] Li | Computing the time evolution of the apparent viscoelastic/viscoplastic Poisson's ratio of hydrating cement paste[END_REF].

The hypothesis of a time-independent Poisson's ratio is nevertheless considered as reasonable for cementitious materials and it is necessary for the extension proposed in the current section. It is formally explained in the following equations.

The aging viscoelastic constitutive phases of the composite materials considered here are supposed to have a time-independent Poisson's ratio ν. The consequences on the relaxation functions K(t, t ) and G(t, t ) are derived by writing that the lateral stress must be null during tensile stress relaxation tests:

0 = σ lat (t) = t -inf ty 2(G(t, t ) + 3K(t, t ))dε lat (t ) + t -inf ty (3K(t, t ) -2G(t, t ))dε(t ) = t -inf ty (-ν2(G(t, t ) + 3K(t, t )) + (3K(t, t ) -2G(t, t )))dε(t )
Since the previous equation is valid for all loading times t and all times t, the spherical relaxation function K(t, t ) varies directly with the deviatoric relaxation function G(t, t ) and:

ν = 3K(t, t ) -2G(t, t ) 2(3K(t, t ) + G(t, t ))
Hence, the relaxation tensor of these materials have the following property [START_REF] Hilton | The significance of (an)isotropic viscoelastic Poisson ratio stress and time dependencies[END_REF]:

g(t, t ) = g(0, 0)f (t, t ) (1)
where g(0, 0) is the instantaneous elastic stiffness at t = 0 and f (t, t ) is a scalar function such that f (0, 0) = 1.

Eshelby's solution for ellipsoidal inclusions

Let us consider here an elastic solid of stiffness g(0, 0) that undergoes a uniform eigenstrain ε * in a bounded domain V of characteristic function χ V (x). Because of the linearity of the equations governing local strain and stress fields, the strain ε(x) at any point x in the body is a linear function of the eigenstrain ε * :

ε(x) = (x) : ε *
This strain field ε(x) is compatible and the corresponding stress field σ(x) is:

σ(x) = g(0, 0) : ((x) -χ V (x)I) : ε *
The stress field σ(x) is self-balanced. Moreover, assuming the elastic body infinite and the domain V of ellipsoidal shape, Eshelby [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] showed that the tensor (x) is uniform on the inclusion:

(x) = x ∈ V
where is Eshelby's tensor [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF]. In the case of an isotropic tensor of module g(0, 0), the tensor only depends on the Poisson's ratio ν of g(0, 0) and on the shape of the ellipsoidal inclusion.

It is now supposed that, in an ellipsoidal domain V of characteristic function χ V (x), an aging viscoelastic solid of relaxation tensor g(0, 0)f (t, t ) undergoes a history of eigenstrain ε * (t). Let us show that the strain field history ε(t, x) = (x) : ε * (t) is the solution to the aging viscoelastic Eshelby's problem. First, it is compatible at all time t. Second, the corresponding stress field history reads:

σ(t, x) = t -∞ g(t, t ) : (dε(t , x) -dε * (t )χ V (x)) σ(t, x) = t -∞ g(t, t ) : ((x) : dε * (t ) -dε * (t )χ V (x))
The property of the aging viscoelastic compliance described in equation 1 is used:

σ(t, x) = t -∞ g(0, 0)f (t, t ) : ((x) : dε * (t ) -dε * (t )χ V (x)) σ(t, x) = g(0, 0) : ((x) -χ V (x)I) : t -∞ f (t, t )dε * (t ) An equivalent eigenstrain history ε * * (t) = t -∞ f (t, t )dε * (t )
is introduced and directly linked to stress field history at the same time:

σ(t, x) = g(0, 0) : ((x) -χ V (x)I) : ε * * (t) Hence, σ(t, x) is self-balanced at all time t. Consequently, ε(t, x) = (x) : ε * (t) is
the solution of the aging viscoelastic Eshelby's problem for this class of materials. The viscoelastic tensor (t, t , x) is similar to the elastic case:

(t, t , x) = H(t -t )(x)
The strain field reads:

ε(t, x) = t -∞ (t, t , x) : dε * (t )
In particular, it must be noticed that the viscoelastic localization tensor is still uniform on the ellipsoid and the viscoelastic Eshelby's tensor (t, t ) = H(t -t ) is defined. For any eigenstrain history ε * (t), the uniform strain field in the inclusion ε(t) reads:

ε(t) = t -∞ (t, t ) : dε * (t ) ε(t) = : ε * (t)
From a more physical point of view, this means that the strain field in this viscoelastic problem is essentially identical to the one of the purely elastic one. The time dependence of this strain field directly follows the time dependence of the prescribed eigenstrain, On the other hand, the stress field is sensitive to the time-dependence of the constitutive relation. This property will be used hereafter to extend standard linear homogenization schemes based on Eshelby's solution for ellipsoidal inclusions to aging visco-elasticity.

Extension of homogenization schemes

Eshelby's reasoning to extend Eshelby's solution to the inhomogeneity problem of an ellipsoidal inclusion having an elastic stiffness differing from the one of the remainder is applied to aging viscoelastic materials in the current section.

Definition of Volterra's tensorial operators and properties

The tensorial Volterra operator described by Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] have interesting properties which ease further computations. Let g a (t, t ) and g b (t, t ) be two relaxation tensors (order 4) and ε(t) be an history of the strain field (tensor of order 2). The tensorial Volterra operator of order 2, noted as: is the operator such that:

g a :ε(t) = t -∞ g a (t, t ) : dε(t )
It must be noticed that:

H(t -t )I:ε = ε
where H(t -t ) is Heaviside's function and I is the unit tensor of order 4. The tensorial Volterra operator of order 4 is such that:

g a :g b (t, t ) = t τ =-∞ g a (t, τ ) : d τ g b (τ, t )
Some properties of these operators are listed in the current section. Their proofs, in Appendix C, rely on the causality principle. Obviously, these operators are bilinear. Moreover, the right and left identity of the tensorial Volterra operator of order 4 is H(t -t )I:

H(t -t )I:g a = g a (2) g a :H(t -t )I = g a (3)
These operators can be associated:

(g a :g b ):ε = g a :(g b :ε) (4)
The tensorial Volterra operator of order 4 is associative:

g a :(g b :g c ) = (g a :g b ):g c (5)
Nevertheless, the Volterra operator is not commutative. The inverse of the relaxation tensor g a is the compliance tensor g -1 a such that:

g -1 a :g a = H(t -t )I
The left and right inverses are equal since the Volterra operator is associative. These notations are applied to the viscoelastic Eshelby's solution. The strain history in the inclusion ε reads:

ε = H:ε *
The stress history in the inclusion reads:

σ = g:(ε -ε * ) 1.4.2

. Extension to the ellipsoidal inclusion

A homogeneous strain field E(t) is added to the viscoelastic Eshelby's solution for a uniform eigenstrain history ε * (t) in an ellipsoidal region of a linear viscoelastic material.

The total strain in the inclusion is ε + E where ε = H:ε * . The stress field in the inclusion is:

σ = g:(ε -ε * + E)
A viscoelastic inhomogeneity of relaxation tensor g i (t, t ), free from eigenstrain, is substituted to the reference material in the ellipsoidal region. To ensure a compatible strain field history at all time, this inclusion feature the same total strain history ε + E. Hence, the stress field in the ellipsoidal region reads:

σ(t) = g i :(ε + E)
To ensure a self-balanced stress field, the stress field in the ellipsoidal inclusion is required to be equal to the one before substitution:

g:(ε -ε * + E) = g i :(ε + E) -g:(ε * ) = (g i -g):(ε + E) -(g i -g):(H:ε * ) -g:(ε * ) = (g i -g):E [-(g i -g):H -g]:(ε * ) = (g i -g):E
This requirement dictates the choice of the eigenstrain history ε * :

ε * = [(g -g i ):H -g] -1 :(g i -g):E
The strain in the inclusion ε l reads:

ε l = H: [(g -g i ):H -g] -1 :(g i -g):E + E = H: [(g -g i ):H -g] -1 :(g i -g) + HI :E
Hence, the expression of a localization tensor for the inhomogeneous inclusion ε l = :E has been derived. It is similar to the elastic case, except for the use of the Volterra operator:

= HI + H:g -1 :(g i -g)

-1 (6)
This result is restricted to a viscoelastic reference material g featuring a time-independent

Poisson's ratio, as described in equation 1. It must be noticed that there is no restriction on g i . Indeed g i can represent any aging linear viscoelastic anisotropic material.

An anisotropic behavior may also be considered for the reference material, as long as the property defined in equation 1 is satisfied. Eshelby's tensor has been explicitly evaluated up to transversely isotropic material [START_REF] Mura | Isotropic inclusions, in: Micromechanics of defects in solids[END_REF] and numerically evaluated in other cases [START_REF] Gavazzi | On the numerical evaluation of Eshelby's tensor and its application to elastoplastic fibrous composites[END_REF].

The particular case of an elastic inclusion g i (t, t ) = g i H(t -t ) leads to further simplifications:

-1 = HI + H:(g -1 :g i H -HI) -1 :ε l = ε l + : (g -1 :(g i : ε l )) -: ε l
The equation -1 :ε l = E where ε l is the unknown is a Volterra integral equation of the first kind.

Homogenization schemes

Mean field homogenization schemes considered in the present study are the Hashin-Shtrikman lower bound [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF] or the Mori-Tanaka scheme [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF], as reconsidered by Benveniste [START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF]. These mean field methods rely on Eshelby's equivalent inclusion theory [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] to estimate the stress concentrations in ellipsoidal inclusions. The elastic Mori-Tanaka estimate g e M T accounts for the volume fraction of inclusions c i and the distribution of orientations of inclusions f (ψ). It is the solution of equation:

c i ψ f (ψ)(g e M T -g i (ψ)) : e (ψ)dψ + (1 -c i )(g e M T -g(0, 0)) = 0 (7)
Here, g i (ψ) is the elastic stiffness of inclusions having orientation ψ ; g(0, 0) is the elastic stiffness of the matrix; e (ψ) is the strain concentration tensor expressing the strain in the inclusions having orientation ψ as a linear function of the strain at infinity, g(0, 0) being the elastic stiffness tensor of the reference material. Tensors g i (ψ) and e (ψ) are computed by rotating g i (0) and e (0) using Bond transformations [START_REF] Bond | The mathematics of the physical properties of crystals[END_REF][START_REF] Auld | Acoustic fields and waves in solids[END_REF].

The equation 7 is extended to the case of viscoelastic materials [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF]:

c i ψ f (ψ)(g M T -g i (ψ)):(ψ)dψ + (1 -c i )(g M T -g) = 0
where (ψ) is derived by rotating (0), which expression is given in equation 6. The relaxation tensor of the overall material g M T is computed as:

g M T = c i ψ f (ψ)g i (ψ):(ψ)dψ + (1 -c i )g : c i ψ f (ψ)(ψ)dψ + (1 -c i )HI -1 (8)
This closed-form formula summarizes the extension of the Mori-Tanaka scheme to aging linear viscoelastic matrices featuring a time-independent Poisson's ratio. It is to be noticed that there is no restriction on the behavior of the inclusions: aging linear viscoelastic anisotropic inclusions can be considered. The alternative Ponte Castañeda-Willis scheme also provides an estimate of the overall stiffness of a composite featuring ellipsoidal inclusions [START_REF] Ponte Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF]. It accounts for inclusion shape and spatial distribution independently. Its extension to linear aging viscoelasticity reads:

g P CW = g + I -c i ψ f (ψ):(ψ)dψ: d -1 : c i ψ f (ψ)(ψ)dψ (9)
where (ψ) = (g i -g) -1 + (ψ) -1 . Hill's tensor (ψ) is derived by rotating (0) = :g -1 and d = d :g -1 where d is Eshelby's tensor corresponding to spherical inclusions. The similarities and differences between the Mori-Tanaka scheme and the Ponte Castañeda-Willis scheme are presented in references [START_REF] Ponte Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF][START_REF] Hu | Some reflections on the Mori-Tanaka and Ponte Castañeda-Willis methods with randomly oriented ellipsoidal inclusions[END_REF][START_REF] Hori | Double-inclusion model and overall moduli of multi-phase composites[END_REF][START_REF] Hu | The connections between the double-inclusion model and the Ponte Castaneda-Willis, Mori-Tanaka, and Kuster-Toksoz models[END_REF].

To turn these closed form formulas (Eq. 8,9) into a practical way to upscale the viscoelastic behavior of composite materials, a procedure to represent strain history and to compute the Volterra operators is defined in the next section.

Applications of the extended Mori-Tanaka scheme

Numerical computation of Volterra's operators

A numerical procedure is now depicted to compute the localization tensor according to equation 6 and to evaluate the overall behavior defined by the extended Mori-Tanaka scheme in equation 8. This procedure is similar to the one of Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF], except that it must be able to deal with tensor histories instead of scalar functions of time.

Following Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] and Bažant [START_REF] Bažant | Numerical determination of long-range stress history from strain history in concrete[END_REF], time is discretized into steps 0 < t 0 < ... < t n and a quadrature rule is used to approximate integrals and solve the Volterra equation.

The strain history ε(t ) and the relaxation function g(t i , t ) are approximated as a piecewise linear functions on [t 0 ; t n ]: for all time t i and at time t ∈ [t j ; t j+1 ], these functions are evaluated as:

ε(t ) = ε(t j ) + t -t j t j+1 -t j (ε(t j+1 ) -ε(t j )) g(t i , t ) = g(t i , t j ) + t -tj tj+1-tj (g(t i , t j+1 ) -g(t i , t j ))
Moreover, it is assumed that ε(t ) = 0 on [0; t 0 [. This difference to the method of Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] is introduced to accurately depict situations where the material is loaded at t 0 , such as relaxation tests or creep tests. The corresponding stress history σ = g:ε is estimated by trapezoidal approximation:

σ(t i ) = g(t i , t 0 )ε(t 0 ) + i-1 j=0 1 2 (g(t i , t j ) + g(t i , t j+1 )) (ε(t j+1 ) -ε(t j ))
Hence, strain history and stress history are represented as vectors ε = [ε(t 0 ), ..., ε(t i ), ..., ε(t n )] and σ = [σ(t 0 ), ..., σ(t i ), ..., σ(t n )]. The relaxation function g is represented as a matrix C such that σ = C.ε. As written by Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF], its generic blocs of size 6 × 6 reads:

2g ij =            2g(t 0 , t 0 ) i = j = 0 g(t i , t 0 ) -g(t i , t 1 ) i > 0, j = 0 g(t i , t j-1 ) -g(t i , t j+1 ) i > 1, 0 < j < i g(t i , t i-1 ) + g(t i , t i ) i > 0, j = i 0 i < j
This matrix is lower triangular per block, as could have been anticipated from the causality principle. An elastic behavior corresponds to a diagonal per block matrix.

Computing and g M T resumes to matrix-matrix products, matrix-vector products and inverting lower triangular per block matrices. It must be noticed that the last step is equivalent to time-stepping if blocked lines are inverted one after another [START_REF] Huet | Adaptation d'un algorithme de Bažant au calcul des multilames viscoélastiques vieillissants[END_REF].

The accuracy of the extended Mori-Tanaka scheme, which uses the numerical procedure described above, is to be checked.

Validation against reference solutions

The method described above must be validated against existing tools whenever possible. The case of aging viscoelastic polymer is interesting due to the fact that the Laplace-Carson method may be used after a change of variable defining an equivalent time. The time-dependent strains estimated by this method are used as references to validate the model described above in case of elongated inclusions.

Aging of polymer by the time-shift method

A composite material made of elastic glass fibers embedded in an aging viscoelastic polymer as in [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] is considered. The compliance of an aging viscoelastic polymer is modeled by the time-shift method and a Kohlrausch-Williams-Watts (KWW) function [START_REF] Kohlrausch | Theorie des elektrischen rückstandes in der leidner flasche[END_REF][START_REF] Williams | Non-symmetrical dielectric relaxation behaviour arising from a simple empirical decay function[END_REF]. It reads:

t(t, t ) = t 0 e t 1+µ -t 1+µ (1+µ)τ µ ref τ 0 m
where t 0 is an isotropic stiffness tensor and m = 1/3 is common to various polymers below the glassy temperature and to metals. If the snapshot condition (t -t ) << t is satisfied, the expression of the compliance resumes to the momentary creep curve:

J(t, t ) ≈ J 0 e t 1+µ (t-t ) τ µ ref τ 0 m
which is similar to the KWW function if t is set. The aging parameters µ and τ ref define an equivalent time λ = t 0 dτ φ(τ ) where φ(t) = (t/τ ref ) -µ . Values of these parameters have been fitted to the experimental creep performed by Read [START_REF] Read | Physical ageing and creep in PVC[END_REF] on Polyvinyl Chloride (PVC). It must be noticed that the ratio between the function φ ZW (t) = 1 + t/b with b = 3000s defined in [START_REF] Zheng | A new constitutive equation for the long-term creep of polymers based on physical aging[END_REF] and the function φ(t) used in the present article is nearly uniform between 1 day and 100 days for PVC, leading to the conclusion that these definitions of the equivalent time are consistent in this range. Yet, the equivalent time defined in [START_REF] Zheng | A new constitutive equation for the long-term creep of polymers based on physical aging[END_REF] can be preferred if extrapolation to long term creep is to be performed.

Furthermore, Kelvin chains are adjusted to the KWW function to ease the forward Laplace-Carson transform. These values from [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] To upscale the viscoelastic response of a non-aging viscoelastic material, the Laplace-Carson transform is combined to mean-field homogenization schemes [START_REF] Hashin | Complex moduli of viscoelastic composites-i. general theory and application to particulate composites[END_REF][START_REF] Lévesque | Numerical inversion of the Laplace-Carson transform applied to homogenization of randomly reinforced linear viscoelastic media[END_REF]. The Laplace-Carson transform turns a non-aging viscoelastic homogenization problem into a set of elastic homogenization problems parametrized by p > 0. The transform of a function g(t) is ĝ(p) = p ∞ 0 g(t)e -pt dt (Appendix A). This transform is still usable on polymer-based composites since aging is defined as an equivalent time.

For a given p, the elastic Mori-Tanaka estimate g p M T accounts for the volume fraction of inclusions c i and the distribution of orientations of inclusions f (ψ). It is the solution of equation:

c i ψ f (ψ)(g p M T -g i (ψ)) : p (ψ)dψ + (1 -c i )(g p M T -g p m ) = 0
Here, g i (ψ) is the elastic stiffness of inclusions having orientation ψ ; g p m is the elastic stiffness of the matrix corresponding to p ; p (ψ) is the strain concentration tensor expressing the strain in the inclusions having orientation ψ as a linear function of the strain at infinity, g p m being the elastic stiffness tensor of the reference material. Tensors g i (ψ) and p (ψ) are computed by rotating g i (0) and p (0) using Bond transformations [START_REF] Bond | The mathematics of the physical properties of crystals[END_REF][START_REF] Auld | Acoustic fields and waves in solids[END_REF]. Formula to compute p (0) in the local reference are recalled in references [START_REF] Mura | Isotropic inclusions, in: Micromechanics of defects in solids[END_REF][START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF][START_REF] Pierard | Micromechanics of elasto-plastic materials reinforced with ellipsoidal inclusions[END_REF].

Results on a fiber reinforced polymer

The effect of large aspect ratios on the time-dependent strain is estimated by the Laplace-Carson method and the extended Mori-Tanaka scheme, the volume fraction of glass fibers being set at 20% and the time of loading at 10 days. A large aspect ratio tends to reduce the time-dependent strains. The estimate of these strains are merely identical for the two methods, which both rely on the Mori-Tanaka model to perform the homogenization step. The relative error between the estimates of strains remain below 1% (Fig. 2). This error arises from the fit of the Kelvin chain and from the Gaver-Stehfest inversion formula for the Laplace-Carson method or from the representation of the compliance as a matrix and subsequent computations for the time-space method.

Since the accuracy of the extended Mori-Tanaka scheme is established, it is applied to aging cementitious materials in the next section. Figure 2: Top: the estimated normalized creep strains are displayed as functions of time elapsed since loading t -t 0 . Time-dependent strains for different aspect ratios a/b, for bulk (left) and shear (right) loadings, are estimated by using the Mori-Tanaka scheme in the Laplace-Carson space and by the Mori-Tanaka scheme in the space domain. Bottom: the relative error between the estimated time-dependent strains is displayed as a function of the time elapsed since loading t -t 0 . The time of loading is 10 days after quenching.

Validation against full 3D numerical simulations 2.3.1. Full 3D numerical simulations

In the frame of periodic homogenization, the determination of the overall viscoelastic behavior of a periodic microstructure can be obtained by solving the following auxiliary problem on the periodic unit cell V .

div σ(x, t) = 0 x ∈ V ε(x, t) = t -∞ t(x, t, t )dσ(x, t ) x ∈ V ε(x, t) = E(t) + ∇ s u(x, t) x ∈ V u(x, t) periodic x ∈ ∂V σ(x, t) • n(x) anti -periodic x ∈ ∂V
Here E(t) is the time-dependent overall strain, u(t, x) is the displacement field in V ,∇ s u(x, t) is its symmetric gradient, ∂V is the boundary of V and n(x) is the outer normal to ∂V .

Actually, E(t) is the volume average of ε(x, t) and we denote by Σ(t) the volume average of σ(x, t). Full 3D numerical simulations have already been performed to upscale mechanical properties of composites in the frame of the periodic homogenization theory. The finite element method or the Fast Fourier Transform (FFT) method [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Šmilauer | Identification of viscoelastic C-S-H behavior in mature cement paste by FFT-based homogenization method[END_REF][START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF][START_REF] Escoda | Modélisation morphologique et micromécanique 3D de matériaux cimentaires[END_REF][START_REF] Sliseris | An accelerated simulation method of medium density wood fiber boards[END_REF] are used to solve elastic problems. The Random Sequential Adsorption algorithm [START_REF] Feder | Random sequential adsorption[END_REF] is used to generate periodic microstructures [START_REF] Bentz | A Hard Core/Soft Shell Microstructural Model for Studying Percolation and Transport in Three-Dimensional Composite Media[END_REF][START_REF] Pan | Analysis of 3d random chopped fiber reinforced composites using FEM and random sequential adsorption[END_REF]. Overlapping between polyhedral inclusions is prevented thanks to the Gilbert-Johnson-Keerthi distance algorithm [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three space[END_REF] as in [START_REF] Lavergne | Investigation of the effect of aggregates' morphology on concrete creep properties by numerical simulations[END_REF]. Convex polyhedral inclusions are defined as the convex hull of 12 points. Polyhedra are then scaled to match the volume of the corresponding inclusion. The resulting microstructure is displayed on figure 3.

The 3D numerical method used in the present article is the one designed by Šmilauer and Bažant [START_REF] Šmilauer | Identification of viscoelastic C-S-H behavior in mature cement paste by FFT-based homogenization method[END_REF] developed for cementitious materials. This method which relies on the exponential algorithm [START_REF] Zienkiewicz | A numerical method of visco-elastic stress analysis[END_REF][START_REF] Taylor | Thermomechanical analysis of viscoelastic solids[END_REF][START_REF] Bažant | Creep and shrinkage law for concrete at variable humidity[END_REF] is a time-iteration procedure to solve the viscoelastic problem for the case of steady loads. It features an integration of the constitutive equations on each time step assuming a uniform stress rate, to enable the time step to grow exponentially when performing a relaxation (or creep) simulation. The exponential algorithm has been recently combined with the FFT algorithm as solver for the unit cell tangent problem [START_REF] Šmilauer | Identification of viscoelastic C-S-H behavior in mature cement paste by FFT-based homogenization method[END_REF][START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF].

The strain ε triggered by a periodic polarization field τ in an homogeneous material of stiffness g * 0 , submitted to the average strain E is given by the Lippman-Schwinger equation:

ε = E - * 0 * τ where * 0 is a Green operator. The convolution of * 0 and τ is computed in the frequency domain thanks to the FFT. The polarity tensor τ is chosen so as to account for the heterogeneity of the considered material g(x):

τ (x) = (g(x) -g * 0 ) : ε(x)
The strain field must satisfy the following equation:

ε = E - * 0 * ((g -g * 0 ) : ε)
This equation is solved by a fixed point algorithm [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]( Appendix B).

Although this 3D numerical method can treat large and complex microstructures, it requires large amount of memory and time. The implementation used in this article is parallel so as to be ran on clusters [START_REF] Lavergne | Investigation of the effect of aggregates' morphology on concrete creep properties by numerical simulations[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].

The case of cementitious materials: the B3 model

The B3 model [START_REF] Bažant | Solidification theory for aging creep[END_REF][START_REF] Bažant | Creep and shrinkage prediction model for analysis and design of concrete structures-model B3[END_REF] At early age, the phenomenon of aging, which is caused by cement hydration and probably also by gradual formation of bonds akin to polymerization, causes major complications for the modeling of concrete creep. The aging of a cement paste is treated as a consequence of volume growth of the load-bearing solidified matter (hydrated cement) whose compliance ¨is non-aging:

¨(t -t ) ∝ ln(1 + λ -n 0 (t -t ) n
) where λ 0 = 1 day. Although a value of n = 0.1 is used at concrete scale [START_REF] Bažant | Solidification theory for aging creep[END_REF][START_REF] Bažant | Creep and shrinkage prediction model for analysis and design of concrete structures-model B3[END_REF], higher values (n = 0, 24) have been used for a cement paste with high water to cement ratio [START_REF] Granger | Effect of composition on basic creep of concrete and cement paste[END_REF]. At the nanometric scale of the calcium-silicate-hydrates (C-S-H), this parameter may be set at n ≈ 0.35 [START_REF] Šmilauer | Identification of viscoelastic C-S-H behavior in mature cement paste by FFT-based homogenization method[END_REF]. The parameter n drives the discrepancy between the elastic response t -t << 1, where ¨(t -t ) ≈ (t -t ) n and the long term trend t -t >> 1, where ¨(t -t ) ≈ n ln(t -t ). The creep tests of Grasley & Lange [START_REF] Grasley | Constitutive modeling of the aging viscoelastic properties of Portland cement paste[END_REF] may be used to set this value.

The early-age compliance t ea (t, t ) of the solidifying cementitious material should satisfy the following conditions:

dt ea dt (t, t ) = 1 v(t) dd t (t -t ) t ea (t, t) = g -1 q1 + 1 v(t) ¨(0)
where v(t) is proportional to the volume fraction of load-bearing material at time t and g q1 is the dynamic part of the elastic stiffness. The inverse of this dynamic Young Modulus is q 1 = 0.6/E 28 where E 28 is the measured elastic stiffness at 28 days. This ratio of 0.6 is consistent with the experimental results of Kawasumi et al. [START_REF] Endo | Determination of the static elastic constant of concrete derived from the elastic constant of cement paste[END_REF] (Tab.2).

According to the B3 model, v(t) reads:

1 v(t) = q 2 λ 0 t + 0.29 * (w/c) 4
where w/c is the water to cement weight ratio. q 2 is related to q 1 and the elastic stiffness at 28 days E 28 . As in [START_REF] Bažant | Creep and shrinkage prediction model for analysis and design of concrete structures-model B3[END_REF], it is assumed that the elastic stiffness E 28 is measured at t -t = 0.01 and t = 28 days. Hence:

q 2 = 0.4 ln(1 + (0.01) n )E 28 1/28 + 0.29(w/c) 4
The Young Modulus increases with time up to a limit:

1 E(t) = 1 E(28)   0.6 + 0.4 1 t + 0.29w/c 4 1 28 + 0.29w/c 4   > 1 E(28)   0.6 + 0.4 0.29w/c 4 1 28 + 0.29w/c 4  
For w/c = 0.5, the B3 model expects the long term Young Modulus to be 30% higher than the Young Modulus at 28 days while experimental results on sealed samples of concretes [START_REF] Granger | Comportement différé du béton dans les enceintes de centrales nucléaires : analyse et modélisation[END_REF] exhibits an increase of about 10% of the Young Modulus between 28 days and 1 year. Moreover, models ACI209, Eurocode2 and the BPEL91 expect respectively the long term Young Modulus to be 8%, 9% and 3% higher than the Young Modulus at 28 days, for a normal cement (R or ASTM I) [START_REF] Berthollet | Contribution à la modélisation du béton vis-à-vis du vieillissement et de la durabilité : interaction des déformations de fluage et du comportement non-linéaire du matériau[END_REF]. Consequently, the function v(t) is modified under the assumption that aging stops once the Young modulus reaches 1.1 times the Young modulus at 28 days and v(t) must remain below a maximum value v max , which reads:

v max = 1 1 1.1 -0.6 E(28)ln(1 + 0.01 n )
The strain rate reads:

εea (t) = g -1 q1 + v(t) -1 ¨(0) σ(t) - t 0 v -1 (t) dd t (t -t ) σ(t )dt
To enable full 3D numerical simulations, the non-aging part of the compliance Φ(t-t ) is considered as a series of Kelvin chains (Tab. 3):

εea (t) = g -1 q1 + v -1 (t)g -1 0 σ(t) + t 0 v -1 (t) n 1 1 τ k e -t-t τ k g -1 k σ(t )dt
The internal variables are:

γ k (t) = t 0 1 τ k e -t-t τ k g -1 k σ(t )dt
The evolution equations are:

εea (t) = g -1 q1 + v -1 (t)g -1 0 σ(t) + v -1 (t) n 1 γ k (t) γk (t) + 1 τ k γ k (t) = g -1 k σ(t)
A irrecoverable part of creep strain ε q4 in cementitious materials is taken account of in the second compliance t q4 (t, t ). Its expression was derived from the theory of the relaxation of microprestress [START_REF] Bažant | Microprestresssolidification theory for concrete creep. i : Aging and drying effects[END_REF]

n 1 1 -e -t-t τ k E -1 k + E -1 0 is fitted according to the non- aging part of the compliance Φ(t -t ) = ln(1 + λ -n 0 (t -t ) n ), where λ 0 = 1 day. cement paste Elastic E 28 13.5 GPa w/c 0.5 q 4 2 10 -6 MPa -1 aggregates Elastic E 0 60 GPa ν 0.25
Table 4: Parameters of the viscoelastic behavior of a cement paste, according to model B3. The Poisson's ratio of the cement paste is assumed to be equal to 0.2.

Pickett effect when the relative humidity changes. This aging compliance features a logarithmic long term trend:

t q4 (t, t ) = ln t t g -1 q4
where the inverse of the Young modulus of g q4 is q 4 = 2.10 -6 MPa -1 for a saturated cement paste [START_REF] Šmilauer | Identification of viscoelastic C-S-H behavior in mature cement paste by FFT-based homogenization method[END_REF]. The corresponding evolution equation is:

εq4 (t) = 1 t g -1 q4 σ(t)
The strain rate ε(t) is the sum of these two strain rates. It results in the following evolution equations:

ε(t) = g -1 q1 + v -1 (t)g -1 0 σ(t) + v -1 (t) n 1 γ k (t) + 1 t g -1 q4 σ(t) γk (t) + 1 τ k γ k (t) = g -1 k σ(t)

Results

The behavior of the matrix is the one of a cement paste (Fig. 4) of water to cement ratio w/c = 0.5 and Young Modulus at 28 days E 28 = 13.5GPa. The inclusions are considered as elastic, as aggregates in concrete (Tab. 4) and the volume fraction of inclusions is 60%. The sieve curve corresponds to the one of a concrete with a small maximal diameter (13.5mm) and a coarse sand (Fig. 3). A 3D microstructure featuring polyhedral inclusions is generated and discretized on a 384×384×384 grid as described in [START_REF] Lavergne | Investigation of the effect of aggregates' morphology on concrete creep properties by numerical simulations[END_REF]. Uniaxial creep tests are simulated, the concrete being loaded at 7, 28 and 365 days. The numerical estimates of the time-dependent strains are compared to the one produced by the time-space method of Sanahuja (Fig. 5). The difference between instantaneous strains is small and the time-space method overestimates the time-dependent strains compared to 3D numerical results. A significant discrepancy exists as the concrete is loaded at early ages (≤7 days). The difference between estimates of time-dependent strains increases with time elapsed since loading but remains acceptable. In conclusion of this comparative study, the method of Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] can be considered as a practical and efficient method to upscale the viscoelastic properties of a concrete featuring round aggregates.

New possibilities: studying influence of aggregates'aspect ratio

The time-space method is now used to study the influence of the aspect ratio of aggregates on the time-dependent strain of concrete, the volume fraction of aggregates being set to 60%. Indicators describing the inclusions' shape such as elongation ratio and flakiness ratio have been defined and measured on real samples of aggregates [START_REF] Kwan | Particle shape analysis of coarse aggregate using digital image processing[END_REF][START_REF] Astm | Standard Practice for Description and Identification of Soils (Visual-Manual Procedure)[END_REF][START_REF] Afnor | Essais pour déterminer les caractéristiques géométriques des granulats -Partie 3 : détermination de la forme des granulats -Coefficient d'aplatissement[END_REF][START_REF] Hu | Shape characterization of concrete aggregate[END_REF]. For an ellipsoidal inclusion with principal axes of lengths a > c > b, the elongation ratio is a/c and the flakiness ratio is b/c. The elongation index and flakiness index correspond to the volume fraction of inclusions featuring respectively an elongation ratio greater than 1.8 and a flakiness ratio lower than 0.6. For normal mix design, the combined (flakiness + elongation) index for coarse aggregates must be limited to 25%, for workability reasons [START_REF] Gambhir | Concrete Technology: Theory and Practice[END_REF]. If the inclusions are spheroids, or ellipsoid of revolution, of semidiameters a, b, b. The aspect ratio is defined as a/b. By using the time-space method, it The normalized creep strains ε x /σ x of a cement paste according to the modified B3 model [START_REF] Bažant | Creep and shrinkage prediction model for analysis and design of concrete structures-model B3[END_REF] are displayed as functions of the time elapsed since loading t -t 0 for different times of loading t 0 . An uniaxial creep test σ x is considered. The Young modulus (t -t 0 ≈ 0.01 days) depends on the time of loading t 0 . The long term trend of creep strains is logarithmic (t -t 0 ≥ 100 days).

is shown in figure 6 that the aspect ratio of aggregates has little effect on the estimated time-dependent strains of concretes, even as the aspect ratio a/b of all aggregates is 0.33 or 3, corresponding to a combined index of 100%. Moreover, full 3D numerical computations are performed on periodic microstructures featuring spherical (a/b = 1) elongated (a/b = 3) or flat (a/b = 0.33) inclusions (Fig. 7). It must be mentioned that a single run of the 3D numerical simulation takes more than one hour using eight nodes of the Athos cluster of EdF Lab [START_REF]Athos EdF R&D[END_REF] while running the time-space method lasts less than five minutes on a common laptop computer. It is shown on figure 8 that the results of the time-space method are consistent with the ones of full 3D numerical simulations. In case of non-spherical inclusions, the estimates of the Ponte Castañeda-Willis scheme are closer to the results of 3D numerical than the ones of Mori-Tanaka. Indeed, both the Ponte Castañeda-Willis scheme and 3D numerical simulation expect the flat inclusions to be slightly more effective than elongated inclusions at reducing creep strains. Yet, the time-dependent strains estimated by 3D numerical simulation are slightly higher than the one of Ponte Castañeda-Willis scheme though the later is a lower bound on the elastic stiffness in the range of elasticity. It may be due to numerical errors since the numerical simulations are performed by using Kelvin chains adjusted to the compliance used by the 20 Ponte Castañeda-Willis scheme. Moreover, in 3D numerical simulations, the behavior of voxels overlapping on two phases is set according to a Reuss bound: the estimated creep strains slightly decreases as the grid size increases [START_REF] Lavergne | Investigation of the effect of aggregates' morphology on concrete creep properties by numerical simulations[END_REF]. Finally, it must be noticed that the considered aspect ratio (a/b = 3 or a/b = 0.33) and volume fraction (c i = 60%) are out of the domain defined by Ponte Castañeda and Willis to enforce the hypothesis of impenetrability of the inclusions [START_REF] Ponte Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF][START_REF] Hu | Some reflections on the Mori-Tanaka and Ponte Castañeda-Willis methods with randomly oriented ellipsoidal inclusions[END_REF]:

c i ≤ a/b if(a/b ≤ 1) c i ≤ 1 a/b 2 if(a/b ≥ 1)
However, the small difference between the estimated time-dependent strains spotted in the case of spherical inclusions is not larger in the case of elongated of flat inclusions. Indeed, both full 3D numerical simulations and homogenization schemes lead to the conclusion that the influence of the aspect ratio of aggregates on the time-depend strains is limited. Consequently, as long as an isotropic distribution of the orientations of aggregates is considered, the aspect ratio of aggregates cannot be held responsible for the discrepancies between the time-dependent strains of concretes of similar formulations.

Discussions

It has been shown that the time-space method is relevant to estimate the timedepend strains of a concrete since these estimates are similar to the one obtained by full 3D numerical simulations (Fig. 5,8). Moreover, the small difference between these 22 estimates and its increase with the time elapsed since loading can be attributed to the low contrast between the elastic stiffness of the aggregates and the tangent stiffness of the cement paste on each time step. Indeed, the contrast of mechanical properties between the phases coupled with the large volume fraction of inclusions may explain such a difference, since it is higher at early age and it increases with the time elapsed since loading. A definition of the contrast between phases featuring isotropic elastic behavior Θ stems from the FFT algorithm [START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF]. It reads:

Θ = max( K max K min , G max G min )
where K min and K min (respectively G min and G max ) are the minimum and maximum of the bulk modulus (respectively shear modulus) in the microstructure. For numerical computations, the contrast of tangent stiffness between the cement paste and the aggregates on the first time step is Θ = 11 at 1 day and Θ = 4.5 at 365 days. 20000 days after loading, it increases to Θ = 14 in both cases. In the range of elasticity, Ghossein and Lévesque [START_REF] Ghossein | A fully automated numerical tool for a comprehensive validation of homogenization models and its application to spherical particles reinforced composites[END_REF][START_REF] Ghossein | A comprehensive validation of analytical homogenization models: The case of ellipsoidal particles reinforced composites[END_REF] have shown that that none of the analytical models they have tested provides accurate estimates for all contrasts and all volume volume fractions of inclusions. Yet, they have checked that the error on the overall elastic strains estimated by the Mori-Tanaka scheme is low if the contrast between the stiffnesses of inclusion and matrix is lower than 20 and for volume fractions of inclusions lower than 50%. Choosing the tangent stiffness to define the contrast between the phases is arguable, since the tangent stiffness depends on the duration of the time step. Yet there is not an obvious definition of the contrast between phases for viscoelastic materials and the duration of time-step can be related to the loading duration and variability, which must be accounted for as viscoelastic materials are considered. The tangent stiffnesses of a series of Kelvin chain, the Generalized Maxwell model and the Zener model are bounded, but it is not the case of the Maxwell model or of the Burger model due to the dashpot in series. In these cases, the tangent stiffness tends to zero for long time steps and the contrast would become infinite if the loading duration were not accounted for. The analytical result presented in this paper is restricted to viscoelastic materials featuring a time-independent Poisson's ratio. Some important materials, such as asphalt concrete [START_REF] Di Benedetto | Three-dimensional lin-ear behavior of bituminous materials: Experiments and modeling[END_REF][START_REF] Nguyen | Prediction of linear viscoelastic behaviour of asphalt mixes from binder properties and reversal[END_REF][START_REF] Kassem | Viscoelastic Poisson's ratio of asphalt mixtures[END_REF][START_REF] Islam | Determining temperature and time dependent Poisson's ratio of asphalt concrete using indirect tension test[END_REF] do not comply with this restriction. Moreover, the effect of bulk viscosity is not considered as rubber [START_REF] Rendek | Amplitude dependence of filler-reinforced rubber: Experiments, constitutive modelling and FEM -implementation[END_REF][START_REF] Fritzen | Reduced basis homogenization of viscoelastic composites[END_REF] is studied and the same assumption is made to model mixed oxides nuclear fuels [START_REF] Siefken | SCDAP/RELAP5/MOD3.3 Code Manual MATPRO A Library of Materials Properties for Light-Water-Accident Analysis[END_REF][START_REF] Largenton | Extension of the nonuniform transformation field analysis to linear viscoelastic composites in the presence of aging and swelling[END_REF], except if hot pressing it to be considered [START_REF] Siefken | SCDAP/RELAP5/MOD3.3 Code Manual MATPRO A Library of Materials Properties for Light-Water-Accident Analysis[END_REF]. The extended Mori-Tanaka scheme proposed in this paper is not able to handle these composite materials due to the hypothesis of a time-independent Poisson's ratio of the matrix. Incremental homogenization approaches [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A timeintegration approach[END_REF][START_REF] Ricaud | Effective properties of linear viscoelastic heterogeneous media: Internal variables formulation and extension to ageing behaviours[END_REF][START_REF] Masson | Incremental homogenization approach for ageing viscoelastic polycrystals[END_REF], full-field computations [START_REF] Li | Modeling the apparent and intrinsic viscoelastic relaxation of hydrating cement paste[END_REF][START_REF] Li | Computing the time evolution of the apparent viscoelastic/viscoplastic Poisson's ratio of hydrating cement paste[END_REF][START_REF] Islam | Determining temperature and time dependent Poisson's ratio of asphalt concrete using indirect tension test[END_REF] or reduced-order models [START_REF] Fritzen | Reduced basis homogenization of viscoelastic composites[END_REF][START_REF] Largenton | Extension of the nonuniform transformation field analysis to linear viscoelastic composites in the presence of aging and swelling[END_REF] proved useful to estimate the overall timedependent strains of such composite materials.

Conclusion

The homogenization method of Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] is extended to the case of elongated or flat inclusions embedded in an aging viscoelastic matrix, as long as the reference material features a relaxation tensor which can be written as g(t, t ) = g(0, 0)f (t, t ). It corresponds to reference materials having a time-independent Poisson's ratio and it is a common assumption due to lack of clear experimental evidence on the evolution of the viscoelastic Poisson's ratio. As for the method of Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF], a practical procedure to derive an estimate of the time-dependent strains of a composite material is proposed.

The method described above has been successfully compared to existing methods, namely a semi-analytical homogenization method using the Laplace-Carson transform for non-aging viscoelastic materials and full 3D numerical simulations. Estimates of the time-dependent strains of a fiber-reinforced polymer and a concrete have been compared and they are found similar. Consequently, these three methods are consistent. The extended method was used to estimate the influence of the aspect ratio of aggregates on the time-dependent strains of concretes: it is shown that the viscoelatic behavior of concrete is not significantly affected by the aspect ratio of aggregates.

The extension presented in this article has many advantages. As in the method of Sanahuja, the output is an estimate of the effective behavior of the material and the computational time is much lower that the one of a single full 3D numerical simulation. Moreover, there is no restriction on the behavior of the inclusions: anisotropic aging viscoelastic materials can be handled as easily as isotropic ones. Regarding the matrix , as long as the reference material features a time-independent Poisson's ratio, any compliance can be incorporated without additional adjustment.

If the viscoelastic Poisson's ratio of the reference material is time-dependent, the method presented in this paper is not valid anymore and extending the method of Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] to the case of elongated or flat inclusions is not straightforward. A perturbation method is currently being investigated to account for a change of Poisson's ratio with time.

Elastic inclusions (volume fraction f i ) are embedded in a viscoelastic matrix modeled by a single Kelvin chain. The relaxation problem reads: div σ(x, t) = 0

x ∈ V σ(x, t) = g i ε(x, t)

x ∈ inclusions σ(x, t) = g m ε(x, t) + τ g m ε(x, t)

x ∈ matrix ε(x, t) = E(t) + ∇ s u(x, t)

x ∈ V u(x, t) periodic x ∈ ∂V σ(x, t) • n(x) anti -periodic x ∈ ∂V

In the Laplace-Carson space, for each p, this set of equation corresponds to the elastic problem: div ( σ(x, p)) = 0 x ∈ V σ(x, p) = g i ε(p)

x ∈ inclusions σ(p) = (1 + pτ )g m ε(x, p)

x ∈ matrix ε(x, p) = E(p) + ∇ s û(x, p)

x ∈ V û(x, p) periodic x ∈ ∂V σ(x, p) • n(x) anti -periodic x ∈ ∂V

The Laplace-Carson transform of the macroscopic stress < σ > (p) is computed by solving the elastic problem and the last stage is inverting this transform. Computing the binomial coefficients requires high precision and the long double type (IEEE 754, decimal on 128 bits) provided it. If M is too low, the formula lacks precision [START_REF] Whitt | A unified framework for numerically inverting Laplace transforms[END_REF]. If M is too large, small errors on ĝ( k ln(2) t ) may trigger large errors on the outcome. M is set to 7.

Therefore, to estimate the response at time t, about 14 elastic computations are required. The Gaver-Stehfest formula does not seem to be practical for FEM since it lacks stability or precision. It is suitable as long as the numerical error in the Laplace-Carson space remains very low. Hence, this formula is useful when semi-analytical estimates by the Mori-Tanaka scheme or the self-consistent scheme are computed in the Laplace-Carson space.

As stated in equation 4, the operator of order 4 and the one of order 2 can be associated:

(g a :g b ):ε = g a :(g b :ε) Let's compute g a :(g b :ε)(t): g a :(g b :ε)(t) = The associativity has been demonstrated.

1 .

 1 Extension of the model of Sanahuja to ellipsoidal inclusions 1.1. Estimating the overall time-dependent strains 1.1.1

Figure 1 :

 1 Figure 1: Left: the inhomogeneity problem solved by Sanahuja[81]. Center and right: the inclusion problem and the inhomogeneity problem solved in the present study.

Figure 3 :

 3 Figure 3: A 4 cm-wide sample of a concrete made of coarse sand (d > 0.5 mm) and small gravel (d < 13.5 mm) built by the RSA algorithm is shown along with its discrezation on a 384 × 384 × 384 grid. The volume fraction of polyhedral inclusions in the cubic cell is 60%.

Figure 4 :

 4 Figure4: The normalized creep strains ε x /σ x of a cement paste according to the modified B3 model[START_REF] Bažant | Creep and shrinkage prediction model for analysis and design of concrete structures-model B3[END_REF] are displayed as functions of the time elapsed since loading t -t 0 for different times of loading t 0 . An uniaxial creep test σ x is considered. The Young modulus (t -t 0 ≈ 0.01 days) depends on the time of loading t 0 . The long term trend of creep strains is logarithmic (t -t 0 ≥ 100 days).

Figure 5 :

 5 Figure5: The normalized creep strains ε x /σ x of a concrete estimated by the time-space method of Sanahuja[START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] and full 3D numerical simulations are displayed as functions of the time elapsed since loading t -t 0 . Considered times of loading are t 0 = 7 days (top left), t 0 = 28 days (top right) and t 0 = 365 days (bottom). The inclusions are assumed to be spherical so as to use the model of Sanahuja. The loading corresponds to an uniaxial creep test.

Figure 6 :

 6 Figure 6: The estimated normalized creep strains ε x /σ x of concrete featuring different aspect ratios of aggregates are displayed as functions of the time elapsed since loading t -t 0 . Considered time of loading are t 0 = 7 days (top left), t 0 = 28 days (top right) and t 0 = 365 days (bottom). In case of uniaxial creep tests, the aspect ratio of aggregates has little effect on the time-dependent strains of concretes estimated by the time-space method.

Figure 7 :

 7 Figure 7: Microstructures featuring spherical (a/b = 1), elongated (a/b = 3) or flat (a/b = 0, 33) aggregates are discretized on 384 × 384 × 384 regular grids to perform full 3D numerical simulations. The volume fraction of aggregates is 60%.

Figure 8 :

 8 Figure8: The normalized creep strains ε x /σ x of concretes estimated by the Mori-Tanaka scheme, the Ponte Castañeda-Willis scheme and full 3D numerical simulations are displayed as functions of the time elapsed since loading t -t 0 . Considered aspect ratios of aggregates for these uniaxial creep tests are a/b = 1 (left), a/b = 3 (center) and a/b = 0.33 (right). The time of loading is t 0 = 28 days. In case of spherical inclusions a/b = 1, the estimate of the Ponte Castañeda-Willis scheme is equal to the one of the Mori-Tanaka scheme.

Appendix A. 2 .

 2 Inverting the Laplace-Carson transformLots of methods are available to invert the Laplace-Carson transform. In the present study, the Gaver-Stehfest formula[START_REF] Stehfest | Algorithm 368: Numerical inversion of Laplace transforms [d5[END_REF] has been used:

  d t g b (t , τ ) : d τ ε(τ ) = t τ =-∞ g a :g b (t, τ ) : d τ ε(τ ) = (g a :g b ):ε(t)As stated in equation 5, the tensorial Volterra operator of order 4 is associative:g a :(g b :g c ) = (g a :g b ):g c First, let's compute g a :(g b :g c )(t, t ) : g a :(g b :g c )(t, t ) = t τ =-∞ g a (t, τ ) : d τ [g b :g c ](τ, t ) = t τ =-∞ g a (t, τ ) : dτ [ τ u=-∞ g b (τ, u) : d u g c (u, t )] = t τ =-∞ g a (t, τ ) : [g b (τ, τ ) : d τ g c (τ, t ) + τ u=-∞ d τ g b (τ, u) : d u g c (u, t )] Due to the causality principle, if u > τ then d τ g b (τ, u) = 0. Hence: g a :(g b :g c )(t, t ) = t τ =-∞ g a (t, τ ) : t u=-∞ d τ g b (τ, u) : d u g c (u, t ) ∞ g a (t, τ ) : d τ g b (τ, u) : d u g c (u, t ) = t u=-∞ [g a :g b (t, u)] : d u g c (u, t )= (g a :g b ):g c (t, t )

  .

		440
		420	t 0 = 0.03 ... 1 ... 10 30 ...	1000 (days)
	) MPa -1 x (10 -6 x /σ ε	340 360 400 380 320		experimental results fitted, KWW
		300		fitted, Kelvin chains
		1	100	10000	1e+06	1e+08
				t-t 0 (s)

Table 1 :

 1 are recalled in table 1 along with properties of glass fibers. Values fitted according to the experimental results of Read et. al.[START_REF] Read | Physical ageing and creep in PVC[END_REF] on rigid PVC, at T=23°C. All Poisson's ratios of hard PVC are assumed to be 0.40. Both the values of the aging parameter µ = -0.958 and the shape parameter m = 0.324 are close to the values identified by Struik or Read et. al.[START_REF] Struik | Physical aging in amorphous polymers and other materials[END_REF][START_REF] Read | Physical ageing and creep in PVC[END_REF]

	µ	τ ref (s)		i	τ i (s)	E i (GPa)
	-0.958	1			0 instantaneous	3.27
	(a) Aging parameters		1 2	4.8 × 10 -6 4.8 × 10 -5	906 387
	E 0 (GPa) τ 0 (s)	m	3	4.8 × 10 -4	198
	3.28	22.5 0.324	4	4.8 × 10 -3	83.9
	(b) KWW function		5 6	4.8 × 10 -2 4.8 × 10 -1	43.0 16.5
	E 0 (GPa)	ν		7	4.8 × 10 0	9.21
	80	0.22		8	4.8 × 10 1	1.23
	(c) Fibers				(d) Kelvin chains

2.2.2. The Mori-Tanaka scheme in the Laplace-Carson space

Table 2 :

 2 has been designed to estimate the mechanical behavior of concrete, including its basic creep. The creep compliance of the B3 model t B3 is the sum of two aging compliances. The first one t ea represents the effect of hydration at early age and the second one t q4 is irrecoverable, while satisfying the logarithmic long term trend of creep strains of concretes. t B3 (t, t ) = t ea (t, t) + t q4 (t, t ) dynamic Young Modulus E d static Young Modulus E s E s /E d According to Kawasumi et al. the static modulus of a cement paste E s may be retrieved from the dynamic Young Modulus E d thanks to the following equation: E s = 0.001628 × (E d /0.0981) 1 .315 where E s and E d are in MPa. The ratio E s /E d is 0.6 at 28 days in the aging viscoelastic model. The ratio being lower for small values of E s means that the static modulus ages more than the dynamic modulus. This dynamic modulus does not age at all in the B3 model[START_REF] Bažant | Solidification theory for aging creep[END_REF][START_REF] Bažant | Creep and shrinkage prediction model for analysis and design of concrete structures-model B3[END_REF] .

	GPa	GPa
	1	0.33	0.30
	5	2.52	0.50
	10	6.27	0.62
	15	10.7	0.71
	20	15.6	0.78
	25	21	0.84

Table 3 :

 3 A Dirichlet series

	, designed to explain long-term aging and model the

  t t =-∞ g a (t, t ) : d t (g b :ε)(t ) =-∞ g a (t, t ) : g b (t , t )dε(t ) + t τ =-∞ d t g b (t , τ ) : dε(τ ) Due to the causality principle, if τ > t then d t g b (t , τ ) = 0. Hence: g b (t , t )dε(t ) + t τ =-∞ d t g b (t , τ ) : dε(τ ) = t τ =-∞ d t g b (t , τ ) : dε(τ )

	=	t t =-∞ g a (t, t ) : d t	t τ =-∞ g b (t , τ ) : dε(τ ) (t )
	t t and : = g

a :(g b :ε)(t) = t t =-∞ g a (t, t ) : t τ =-∞ d t g b (t , τ ) : d τ ε(τ ) = t τ =-∞ t t =-∞ g a (t, t ) :

Appendix A. The correspondence principle and the Laplace-Carson transform

Appendix A.1. The correspondence principle The non-aging linear viscoelastic problem corresponds to elastic problems thanks to the Laplace-Carson transform. The transform of a function g(t) is ĝ(p) = p ∞ 0 g(t)e -pt dt. The transform of its derivative ġ(t) is ĝ(p) = pĝ(p) -p.g(0).

Appendix B. FFT solvers

The basic FFT algorithm (Alg. 1) is the one of Moulinec and Suquet as described in [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. The accelerated FFT algorithm (Alg. 2) is the one of Eyre and Milton [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF]. These algorithms are particular cases of the polarization-based scheme [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF][START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF]. The error on equilibrium error eq and 0 (ξ) are computed as described in reference [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. If the loading is a macroscopic stress Σ instead of a macroscopic strain E, the macroscopic strain E i is modified at each time step as E i = g 0 (Σ -σ i ) + ε i , where σ i denotes the volume average of σ i . In this case, the error on boundary condition is modified as

Algorithm 1 Basic FFT scheme Initial strain field ε 0 (x) and prestress σ 0 (x) are provided Initial strain field ε-1 (ξ) is set while error eq > 10 -7 × error eq,0 do for x ∈ points do

end for σi ← FFT(σ) Compute error on equilibrium error eq for ξ ∈ frequencies do ε(ξ

Appendix C. Definition of Volterra's operators and properties

The properties of Volterra's operators listed in section 1.4.1 are demonstrated in the current appendix.

The left and right identity of the tensorial Volterra operator of order 4 is H(t -t )I, as stated in equation 2 and 3:

Due to the causality principle, if τ < t then g a (τ, t ) = 0. Hence:

Algorithm 2 Accelerated FFT scheme Initial strain field ε 0 (x) and prestress σ 0 are provided Initial strain field ε-1 (ξ) is set while error eq > 10 -7 × error eq,0 or error bc > 10 -7 × error bc,0 or error comp > 10 -7 do