F Lavergne

Karam Sab
email: karam.sab@enpc.fr

Julien Sanahuja

Michel Bornert

Charles Toulemonde

Estimation of creep strain and creep failure of a glass reinforced plastic by semi-analytical methods and 3D numerical simulations

Keywords: creep, polymer, aging, homogenization

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Estimation of creep strain and creep failure of a glass reinforced plastic by semi-analytical methods and 3D numerical simulations.

Introduction

Polyvinyl chloride (PVC) is a material of choice for construction applications, such as pipes for water and gas, house sidings or window frames where long service life is required [1]. The European demand of PVC was 5000 ktons in 2012 [START_REF]PlasticsEurope, Plastics -the facts 2013 an analysis of european latest plastics production, demand and waste data[END_REF]. Many polymers may be part of composite materials, such as fiber-reinforced materials. For instance, glass fibers may be incorporated into a polymer matrix Postprint of Mechanics of Materials doi:10.1016/j.mechmat.2015.06.005 such as PVC to increase stiffness, creep resistance and dimensional stability [START_REF] Rahrig | Glass fiber reinforced vinyl chloride polymer products and process for their preparation[END_REF][START_REF] Kinson | Compositions de polychlorure de vinyle renforcée par des fibres de verre avec stabilité dimensionnelle et résistance à la traction[END_REF][START_REF] O'brien-Bernini | Wet use chopped strand glass as reinforcement in extruded products, wO Patent App[END_REF]. The effect on tensile strength of glass fibers depends on the volume fraction of inclusions, on the fiber length and on fiber orientation. These parameters can be significantly affected by processing operations [START_REF] Tungjitpornkull | Processing technique and fiber orientation angle affecting the mechanical properties of e-glass fiber reinforced wood/PVC composites[END_REF][START_REF] Hohe | Modeling of uncertainties in long fiber reinforced thermoplastics[END_REF].

In case of water pipes, different permanent loads coexist. Cooling after extrusion could trigger internal tensile stress in the range 1.5 to 4.8MPa on the inner diameter, the water pressure induces a permanent hoop stress and additional stresses can occur as a result of non-uniform soil settlement [START_REF] Breen | t Veld, Expected lifetime of existing PVC water systems ; summary[END_REF]. The lifetime of such pipes depends on the operating conditions and on the mechanical properties of the constitutive materials.

Glass reinforced plastic pipes may feature a complex structure, including chopped strand mat layers on the inner side and filament wound layers on the outer side [START_REF] Guedes | On the prediction of long-term creep-failure of GRP pipes in aqueous environment[END_REF][START_REF] Diniz Melo | Mechanical behavior of GRP pressure pipes with addition of quartz sand filler[END_REF][START_REF] Rafiee | Simulation of functional failure in GRP mortar pipes[END_REF]. Short-time hydraulic failure occurs on the inner diameter [START_REF] Diniz Melo | Mechanical behavior of GRP pressure pipes with addition of quartz sand filler[END_REF]. Ring deflection tests may be performed according to standard ISO 9967 to estimate the time-dependent strain of a pipe and study the long term creep failure [START_REF] Afnor | Determination du taux de fluage pour les tubes en matières thermoplastiques[END_REF][START_REF] Guedes | On the prediction of long-term creep-failure of GRP pipes in aqueous environment[END_REF]. In this case, failure mode was always the same, fiber breakage, and always localized in the same region, on the inner diameter [START_REF] Guedes | On the prediction of long-term creep-failure of GRP pipes in aqueous environment[END_REF]. The strain at failure is almost constant, slightly decreasing with the time elapsed since loading. An internal pressure test is defined by the standard ASTM D2992 to estimate the long term static hydrostatic strength of glass-fiber-reinforced pipes [START_REF] Astm | Standard practice for obtaining hydrostatic or pressure design basis for fiberglass (glass-fiber-reinforced thermosetting-resin) pipe and fittings[END_REF]. The time-to-failure is expected to depend on the pressure level, following a power law. Hence, studying the viscoelastic behavior of glass reinforced plastics is required to adjust delayed fracture criterion [START_REF] Guedes | Analysis of a delayed fracture criterion for lifetime prediction of viscoelastic polymer materials[END_REF] and estimate the lifetime and durability of such pipes.

Polymers exhibit a viscoelastic response that depends on the time of loading t . This phenomenon is described as physical aging. It is well known that the viscoelastic properties of polymer glasses are significantly influenced by a physical aging process [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]. This is the reason why the standard ISO 9967 specifies both the temperature (23+ -2°C) and the age of the samples at loading (21+ -2 days), measured since quenching. Experimental evidences on creep properties of various polymer materials below the glassy temperature have been gathered by Struik [START_REF] Struik | Physical aging in amorphous polymers and other materials[END_REF] and further creep tests have been performed on PVC since [START_REF] Read | Physical ageing and creep in PVC[END_REF][START_REF] Zhou | Creep performance of PVC aged at temperature relatively close to glass transition temperature[END_REF]. A PVC quenched from 90°C to 20°C at starting time t 0 = 0 is considered in the present study.

This article is focused on the aging creep and failure of the chopped mat strand layer of the pipe, modeled as a composite material made of E-glass fibers and PVC. Homogenization methods are designed to upscale the viscoelastic behaviors of such materials. Mean-field homogenization schemes such as the one of Mori-Tanaka [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF][START_REF] Wang | The influence of inclusion shape on the overall viscoelastic behavior of composites[END_REF][START_REF] Li | Micromechanical modeling of viscoelastic properties of carbon nanotube-reinforced polymer composites[END_REF] or coupling the single fiber problem and the rule of mixture [START_REF] Hohe | Modeling of uncertainties in long fiber reinforced thermoplastics[END_REF] produce estimates of these elastic properties while taking account of fiber aspect ratio and distribution of orientations. Complete modeling of the failure of fiber-reinforced plastics, including damage of fibers and matrix, have been performed by Sasayama et. al. [START_REF] Sasayama | Prediction of failure properties of injection-molded short glass fiber-reinforced polyamide 6,6[END_REF] and Hashimoto et. al. [START_REF] Hashimoto | Prediction of tensile strength of discontinuous carbon fiber/polypropylene composite with fiber orientation distribution[END_REF]. Some of these micromechanical approaches have recently been able to treat aging linear viscoelastic materials and their accuracy is to be checked by numerical simulations [START_REF] Ricaud | Effective properties of linear viscoelastic heterogeneous media: Internal variables formulation and extension to ageing behaviours[END_REF][START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF][START_REF] Šmilauer | Identification of viscoelastic C-S-H behavior in mature cement paste by FFT-based homogenization method[END_REF].

Mean-field homogenization schemes and 3D numerical simulations of matrixinclusion materials have been compared in the range of elasticity [START_REF] Iii | Stiffness predictions for unidirectional short-fiber composites: Review and evaluation[END_REF][START_REF] Moussaddy | A new definition of the representative volume element in numerical homogenization problems and its application to the performance evaluation of analytical homogenization models[END_REF][START_REF] Ghossein | A comprehensive validation of analytical homogenization models: The case of ellipsoidal particles reinforced composites[END_REF][START_REF] Moumen | Effect of reinforcement shape on physical properties and representative volume element of particles-reinforced composites: Statistical and numerical approaches[END_REF]. The case of elongated or flat inclusions has been explored [START_REF] Iii | Stiffness predictions for unidirectional short-fiber composites: Review and evaluation[END_REF][START_REF] Moussaddy | A new definition of the representative volume element in numerical homogenization problems and its application to the performance evaluation of analytical homogenization models[END_REF][START_REF] Ghossein | A comprehensive validation of analytical homogenization models: The case of ellipsoidal particles reinforced composites[END_REF] with the conclusion that Lielen's model [START_REF] Lielens | Prediction of thermomechanical properties for compression moulded composites[END_REF] was the most accurate one provided that the inclusion is stiffer that the matrix. Such comparisons have also been performed on viscoelastic [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A timeintegration approach[END_REF], plastic [START_REF] Pierard | Micromechanics of elastoplastic materials reinforced with ellipsoidal inclusions[END_REF] or viscoplastic [START_REF] Pierard | Micromechanics of particle-reinforced elastoviscoplastic composites: Finite element simulations versus affine homogenization[END_REF][START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A timeintegration approach[END_REF] matrices to assess the accuracy and performances of different methods. Lahellec and Suquet [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A timeintegration approach[END_REF] have used 2D full finite element simulations to validate a semi-analytical scheme which combines the Hashin-Shtrikman estimate for periodic materials and a time-stepping procedure to compute the in-plane time-dependent strain of aligned circular fibers in a viscoelastic matrix. The range of applications of such models depends on the contrast between phases and volume fraction of inclusions.

The objective of the present article is to estimate the aging viscoelastic behavior and creep failure of a E-glass fiber reinforced PVC. Therefore, the Mori-Tanaka scheme is to be compared to full-field computations. The main features of the model presented in this article are :

1. Both the Mori-Tanaka scheme and the 3D numerical procedure rely on a time-shift procedure to account for the aging creep of PVC. 2. The Mori-Tanaka scheme and a 3D numerical procedure are presented and successfully compared. 3. The influence of volume fraction of fibers and length of fibers are investigated. 4. An estimate of creep failure related to the strength of glass fiber, based on the largest principal stress (Rankine criterion), is proposed.

The behavior of each phase is described and the microstructure is presented in the first section. The methods to obtain the overall property are briefly recalled. The outputs of these methods are compared in the second section and the effect of the length of fibers and of the volume fraction of fibers on the time dependent strain and creep failure are estimated.

Microstructure : geometry and mechanical properties

A short-fiber reinforced plastic

The polymer-based composite that is described in the present article will be a mix of hard PVC and short E-glass fibers. Figure 1a shows the microstructure of the considered material [START_REF] Rahrig | Glass fiber reinforced vinyl chloride polymer products and process for their preparation[END_REF].

In the considered fiber-reinforced material, elastic inclusions are added to the viscoelastic matrix, their elastic stiffness being the one of E-glass fibers : the Young modulus is E = 80GPa and the Poisson's ratio is ν = 0.22 [START_REF] Wallenberger | Glass fibers, asm handbook[END_REF][START_REF] Agy | High strength glass fibers[END_REF][START_REF] Mounier | Evaluation of transverse elastic properties of fibers used in composite materials by laser resonant ultrasound spectroscopy[END_REF]. These inclusions are assumed to be cylindrical chopped fibers, with a diameter of 10µm and a length of 100µm. The distribution of their directions is chosen as isotropic and their volume fraction is 15%. The density of E-glass being about 2.55 and the one of PVC being about 1.35, a weight ratio of 34 phr (part per hundred part of resin) corresponds to a volume fraction of 15%. Above 50 phr (or 20% volume fraction), processing difficulties may appear during the extrusion [START_REF] Rahrig | Glass fiber reinforced vinyl chloride polymer products and process for their preparation[END_REF].

The considered chopped strand mat features shorter fibers than the one used in automotive applications [START_REF] Haque | Polymer/wucs mat for use in automotive applications[END_REF] (11mm to 75mm) or pipes, comparable to fibers used for flooring materials(0.2mm to 1mm) [START_REF] Nakano | Use of polyvinyl chlorine resin sheets as flooring material[END_REF]. The considered microstructure could be similar to the one obtained by the use of chopped strand glass as reinforcement in extruded products, where the integrity of the dry glass fiber strands is broken down by a screw extruder [START_REF] Yuan | The effect of high pressure on the mechanical behavior of short fiber composites[END_REF][START_REF] O'brien-Bernini | Wet use chopped strand glass as reinforcement in extruded products, wO Patent App[END_REF][START_REF] Tungjitpornkull | Sombatsompop, Mechanical characterization of e-chopped strand glass fiber reinforced wood/PVC composites[END_REF][START_REF] Phelps | A model for fiber length attrition in injection-molded long-fiber composites[END_REF].

Failure mode

Experimental tests have shown that the lifetime of a pipe under permanent loading decreases as the level of the applied loading increases [START_REF] Afnor | Determination du taux de fluage pour les tubes en matières thermoplastiques[END_REF][START_REF] Guedes | On the prediction of long-term creep-failure of GRP pipes in aqueous environment[END_REF]. Moreover, the strain at failure decreases slightly with respect to lifetime. The failure mode of a glass fiber-reinforced PVC highly depends on the coupling agent [START_REF] Rahrig | Glass fiber reinforced vinyl chloride polymer products and process for their preparation[END_REF] : the observed failure mode of water pipes was fiber breakage located at points where the largest tensile stress occurred [START_REF] Guedes | On the prediction of long-term creep-failure of GRP pipes in aqueous environment[END_REF]. These experimental evidences will drive the proposed estimates of the lifetime of a pipe using the Mori-Tanaka procedure described in this paper. Indeed, it is assumed that the fibers are perfectly bonded to the matrix until failure occurs.

A fiber is considered as brittle, and a maximum stress criterion is defined to describe its failure [START_REF] Hashimoto | Prediction of tensile strength of discontinuous carbon fiber/polypropylene composite with fiber orientation distribution[END_REF]. The large variability in strength found in brittle fibers [START_REF] Chiao | Stress-rupture of s-glass epoxy multifilament strands[END_REF] is well modeled by Weibull-Poisson statistics and is due to various random flaws on the fiber surface [START_REF] Phoenix | Modeling the statistical lifetime of glass fiber/polymer matrix composites in tension[END_REF]. Moreover, the strength of glass fibers is a function of time when subject to permanent loads [START_REF] Barbero | Phenomenological prediction of tensile strength of e-glass composites from available aging and stress corrosion data[END_REF] : moisture ingress in glass fiber reinforced polymers increases stress corrosion cracking in the fibers [START_REF] Taylor | Mechanism of fracture of glass and similar brittle solids[END_REF][START_REF] Stuart | Dependence of ultimate strength of glass under constant load on temperature, ambient atmosphere, and time[END_REF][START_REF] Charles | Static fatigue of glass. i[END_REF][START_REF] Charles | Static fatigue of glass. ii[END_REF][START_REF] Khennane | Durability of glass polymer composites subject to stress corrosion[END_REF] and thus shortening their lifetime in humid or alkali environments [START_REF] Ren | Investigation on stress-rupture behavior of a chopped-glass-fiber composite for automotive durability design criteria[END_REF][START_REF] Phillips | Prediction of long-term stress-rupture life for glass fibre-reinforced polyester composites in air and in aqueous environments[END_REF]. Nevertheless, for the seek of simplicity, these features will be ignored in the present study and the same uniform-in-time maximum stress criterion σ c will be set for all fibers.

Let σ R (σ) be the largest principal stress of the stress tensor σ. Then the scalar Rankine criterion stipulates that the failure occurs if σ R (σ) > σ c . This criterion coincides with Hashimoto's one [START_REF] Hashimoto | Prediction of tensile strength of discontinuous carbon fiber/polypropylene composite with fiber orientation distribution[END_REF] in case of uniaxial tensile stress in the fiber.

Aging linear viscoelasticity

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

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

J(t, t ) = Φ(t -t )
Non-aging compliance may be approximated by a series of Kelvin chain, which arises from a rheological model made of springs and dashpots (Fig. 2). The advantage of a series of Kelvin chain is that internal variables may be defined, which eases numerical computations. Its compliance writes :

Φ K (t -t ) = n 1 1 -e -t-t τ k C -1 k + C -1 0
where C 0 is the elastic stiffness (order four tensor) and for each 1 ≤ k ≤ n, C k is the stiffness corresponding to the characteristic time τ k . Hence, to define a composite material, it is necessary to describe its microstructure and to depict the behavior of each phase.

Polymers : Temperature shift and Time shift approach

A time-shift approach is designed to model the aging viscoelastic behavior of the matrix. See Grasley & Lange [START_REF] Grasley | Constitutive modeling of the aging viscoelastic properties of portland cement paste[END_REF] for a description of this approach, which has been used to model cement paste. Samples are loaded at different times t elapsed since quenching and the time dependent strains are measured. The obtained strain curves plotted as functions of log(t -t ) are identical up to a shift in the horizontal direction, defining a master curve Φ, so that the aging compliance is written as :

J(t, t ) = Φ(ξ(t) -ξ(t ))
where ξ(t) is a pseudo-time. Struik [START_REF] Struik | Physical aging in amorphous polymers and other materials[END_REF], following Kohlrausch [START_REF] Kohlrausch | Theorie des elektrischen rückstandes in der leidner flasche[END_REF], found the Kohlrausch-Williams-Watts (KWW) [START_REF] Williams | Non-symmetrical dielectric relaxation behaviour arising from a simple empirical decay function[END_REF] function

Φ KW W (t) = J 0 e t τ 0 m
where J 0 is an isotropic stiffness tensor and m = 1/3 is common to various polymers below the glassy temperature and to metals. This function is also valid for PVC at 63 • C [START_REF] Zhou | Creep performance of PVC aged at temperature relatively close to glass transition temperature[END_REF].

The pseudo time ξ(t) is given by [START_REF] Morland | Stress analysis for linear viscoelastic materials with temperature variation[END_REF] :

ξ(t) -ξ(t ) = t t dv A(v)
where A(t) is an activity. For instance [START_REF] Grasley | Constitutive modeling of the aging viscoelastic properties of portland cement paste[END_REF],

A(t) = t t ref -µ T (t) T ref -µ T
where µ is the time shift factor and µ T is a temperature shift factor. A change of variable u = ξ(t), σ * (u) = σ(t) where u is the equivalent time may be performed to turn the aging problem into a non-aging one :

ε(t) = ξ(t) 0 Φ(ξ(t) -v) dσ * du (v)dv 1.5.

Identification of the viscoelastic parameters

The identification of the parameters of the constitutive law is performed in two steps. First, the aging parameter µ KW W and the parameters of the KWW function m, τ 0 and J 0 are fitted at once according to the experimental results. Then, to enable 3D numerical computations, a series of Kelvin chain is fitted according to the KWW function.

To perform the first step, the Levenberg-Marquardt algorithm [START_REF] Levenberg | A method for the solution of certain non-linear problems in least squares[END_REF][START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF] is used, as implemented in the gnuplot software [START_REF] Williams | Gnuplot 4.6: an interactive plotting program[END_REF]. This algorithm minimizes the scalar

χ 2 = i j Φ KW W (ξ(t ij ) -ξ(t i )) -y ij δy ij 2
where y ij is the experimental strain measured at time t ij on the sample loaded at time t i and δy ij is the standard deviation of y ij . δy ij could be interpreted as an estimation of the error on y ij , or as a weight w ij = 1/ δy ij . These weights are chosen so that all curves have the same weight and all decades of a given curve have the same weight for fitting. Let N i be the number of points in a creep curve and N ij be the number of points of that curve between t ij / √ 10 and √ 10t ij . The corresponding weight worths w ij = 1/(N ij .N i ).

The procedure described in the present section was applied to the digitalized experimental results of Struik [START_REF] Struik | Physical aging in amorphous polymers and other materials[END_REF] (Fig. 3)(Tab. 1) and to those of Read et. al. [START_REF] Read | Physical ageing and creep in PVC[END_REF]. The obtained results are presented in figure 4 and table 2. Note that the identified values of µ and m are similar to those identified by Read et. al. [START_REF] Read | Physical ageing and creep in PVC[END_REF].

For the second step, a series of Kelvin chain Φ K is identified to the master curve Φ KW W , keeping the same equivalent time ξ(t). Hence, the aging parameter µ K = µ KW W is left unchanged. As performed in the collocation method [START_REF] Schapery | Approximate methods of transform inversion for viscoelastic stress analysis[END_REF][START_REF] Cost | A multidata method of approximate Laplace transform inversion[END_REF][START_REF] Lévesque | Numerical inversion of the Laplace-Carson transform applied to homogenization of randomly reinforced linear viscoelastic media[END_REF], characteristic times τ i are chosen as terms of a geometric sequence, with one term per decade, and corresponding Young modulus E i are required to be positive. This requirement is sufficient to ensure the thermodynamic correctness of the identified compliance [START_REF] Lévesque | Numerical inversion of the Laplace-Carson transform applied to homogenization of randomly reinforced linear viscoelastic media[END_REF]. The KWW function is sampled at equivalent times in a geometric sequence, between half the smallest characteristic time and twice the largest one, so that each decade has the same weight.

For a regular uniaxial loading, the representation of the KWW function by a series of Kelvin chain is accurate as shown in figure 3 and4. It must be noticed that identifying directly µ and E i on the experimental data with the Levenberg-Marquardt algorithm has been tried unsuccessfully.

Homogenization of aging viscoelastic materials

In an heterogeneous material, the compliance J(x, t, t ) depends on the spacial position x. Mean-field homogenization schemes and numerical simulations are designed to estimate the overall compliance of the composite material.

Micromechanical models

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 in this study since aging is defined as an equivalent time.

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. Weng [START_REF] Weng | The theoretical connection between Mori-Tanaka's theory and the hashinshtrikman-walpole bounds[END_REF], extending the results of Zhao and Tandon [START_REF] Zhao | Elastic moduli for a class of porous materials[END_REF], have proven that the effective moduli of the composite containing either aligned or random oriented, identically shaped ellipsoidal inclusions, as estimated by the Mori-Tanaka scheme, have the same expression as those of the Hashin-Shtrikman-Walpole bounds, only with the latter's comparison material identified as the matrix phase and Eshelby's tensor interpreted according to the appropriate inclusion shape. For a given p, the elastic Mori-Tanaka estimate C 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 (ψ)(C p M T -C i (ψ)) : T p (ψ)dψ + (1 -c i )(C p M T -C p m ) = 0
Here, C i (ψ) is the elastic stiffness of inclusions having orientation ψ ; C p m is the elastic stiffness of the matrix corresponding to p ; T p (ψ) is the strain concentration tensor expressing the strain in the inclusions having orientation ψ as a linear function of the strain at infinity, C p m being the elastic stiffness tensor of the reference material. Tensors C i (ψ) and T p (ψ) are computed by rotating C i (0) and T 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 T p (0) in the local reference are recalled in references [START_REF] Mura | Isotropic inclusions[END_REF][START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF][START_REF] Pierard | Micromechanics of elastoplastic materials reinforced with ellipsoidal inclusions[END_REF].

The average stress in the inclusions of orientation ψ is expressed as a linear function of the overall strain through the localization tensor B p (ψ) given by :

B p (ψ) = C i (ψ) : T p (ψ) : c i ψ f (ψ )T p (ψ )dψ + (1 -c i )1 -1
Once the Laplace-Carson transform is inverted, the Rankine criterion of the average stress in the inclusions of orientation ψ is computed to detect the failure of these inclusions. The Rankine criterion is computed by using the routine dsyev() of the LAPACK package [START_REF] Anderson | Lapack: A portable linear algebra library for high-performance computers[END_REF].

It should be mentioned that many homogenization schemes have recently been adapted by Sanahuja [START_REF] Sanahuja | Effective behaviour of ageing linear viscoelastic composites: Homogenization approach[END_REF] to treat arbitrary aging compliances such as the KWW compliance, or any interpolation of experimental points : instead of using the Laplace-Carson transform, this method operates in the time domain. The expression of the strain localization operator requires computations of Volterra integrals by using trapezoidal rules and inversions of lower triangular per block matrices [START_REF] Bažant | Numerical determination of long-range stress history from strain history in concrete[END_REF][START_REF] Huet | Adaptation d'un algorithme de bazant au calcul des multilames visco-élastiques vieillissants[END_REF]. Sanahuja's method handles the case of matrix-inclusion composites with spherical inclusions and we are currently investigating its extension to matrix-inclusion composites with ellipsoidal inclusions.

3D numerical computations

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 0 J(x, t, t ) dσ dt (x, t )dt 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).

3D numerical computations 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 | Numerical prediction of the stiffness and strength of medium density fiberboards[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](Fig. 1b). 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]. The distribution of orientations of fibers is isotropic : the direction of each fiber is randomly picked on the unit sphere. To test if the distribution of the orientations of fibers is isotropic, a 400µm-wide microstructure featuring 1223 fibers is built and chi square tests are performed. Fiber orientations are binned into 20 sectors of equal angle and the estimated number of fibers in each sector (n f,e = 61, 1) is compared to the observed ones n f,o (θ), θ ∈ 1..20 (Fig. 5). The chi-squared test statistic

χ 2 = 20 θ=1 (n f,o (θ) -n f,e ) 2 /
n f,e is computed and is found to be 9.7, 14.5 and 18.6 depending on the axis chosen to split the sectors (x, y and z). If χ 2 follows a chi-squared distribution of 19 degrees of freedom, there is a 70% probability that χ 2 is lower than 21.7. Hence, the obtained values of χ 2 are not surprising and an isotropic distribution of the orientations of fibers can produce such a set of observations. 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 constant stress rate, to enable the time step to grow exponentially when performing a relaxation (or creep) simulation. This assumption is adapted to treat the case of a time-shift aging compliance based on a series of Kelvin chains, which writes:

J (t, t ) = k≥1 1 -e -t µ+1 -t µ+1 (µ+1)t µ ref τ k C -1 k + C -1 0
The internal variables are:

γ k (t) = t 0 t µ τ k t µ ref e -t µ+1 -t µ+1 (µ+1)t µ ref τ k C -1 k σ(t )dt
and the evolution equations write:

ε(t) = C -1 0 σ(t) + k≥1 γ k γk (t) + t µ τ k t µ ref γ k (t) - µ t γ k (t) = t µ τ k t µ ref C -1 k σ(t)
The constitutive equation is integrated on the time step [t i ; t i+1 ] under the following assumption:

σ(t) = σ(t i ) + ∆σ t µ+1 -t µ+1 i t µ+1 i+1 -t µ+1 i
where µ is the aging parameter (Appendix B). This assumption is unchanged for non-aging material (µ = 0), for which the stress varies linearly on the time step.

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 FFT algorithm requires the microstructure to be discretized on a regular grid. Consequently, issues regarding the automatic generation of high-quality adapted meshes for the finite element method are avoided.

The strain ε triggered by a periodic polarization field τ in an homogeneous material of stiffness C * 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 C(x):

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

ε = E -Γ * 0 * ((C -C * 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 C).

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].

Such 3D computations may deliver more precise estimates of the overall behavior than mean-field methods since they rely on an accurate description of the microstructure. Moreover, a distribution of stress concentrations per phase may be retrieved. Large scale computations are required to match both the need for a Representative Elementary Volume and the need for a precise description of the microstructural details.

A careful assessment of the quality of 3D numerical computations is performed in the next section.

Convergence study Time discretization

To assess the accuracy of the integration of the constitutive law on the time step, the response of the matrix to an uniaxial relaxation test σ num (t) is numerically computed for different rate of growth of the time step b = (t i+1 -t i )/(t it i-1 ) (Fig. 6). The Laplace-Carson transform provides a reference σ LC (t) to compute a relative error (σ num (t) -σ LC (t))/σ LC (t). This relative error is found to be very small, as long as the rate of growth remains limited. A rate of growth of 1.118 is set for further computations. Approximately 180 time steps are needed to compute the response between t 1 -t 0 = 10 -3 s and t -t 0 =50 days, the materials being loaded t 0 = 10 days after quenching.

Space discretization

The FFT algorithm solves the tangent elastic problem on each time step. To perform fast Fourier transform, strain ε(t, x), stress σ(t, x) and internal variables are stored on a regular grid of N × N × N points, where N is called the grid size. The microstructure is also discretized on a regular grid (Fig. 7). Hence, some information about the microstructure is lost, which could trigger an error due to discretization. In order to lower the error, for each voxel, a local volume fraction based on 64 sensing points is computed and a Reuss-like constitutive law is computed and assigned to each voxel [START_REF] Lavergne | Investigation of the effect of aggregates' morphology on concrete creep properties by numerical simulations[END_REF]. Therefore, the microstructure is made of black voxels (pure matrix), white voxels (pure inclusions) and gray voxels (composite).

To estimate the error due to space discretization, a given microstructure, l = 133µm in length, is discretized at different grid step (Fig. 7,9) and the response to a shear creep test are computed. Grid sizes N and numerical performances are provided in table 3: the largest grid is N = 648 and the numerical shear creep test took less than a day thanks to parallel computing. The probability distribution function of the Von Mises stress depicts the magnitude of stress concentrations within the microstructure and it is a matter of concern for further durability assessment, especially if non-linear phenomenon were to be considered. The estimation of this distribution depends on the grid size N and it becomes more accurate as the grid size increases (Fig. 9).

Representative Elementary Volume

It is well-known that the asymptotic overall response should not depend on the generated sample neither on the size of the unit cell l which must be large enough to be representative of the microstructure ( [START_REF] Sab | On the homogenization and the simulation of random materials[END_REF][START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF] among others). A size of the unit cell l = 200µm of two times the length of the fibers is chosen and the unit cell is discretized on a N = 384 grid so that the error due to representativity is of the same magnitude as the error due to discretization. To assess the error due to representativity, different sizes of the unit cell l and nine samples of each size were generated (Fig. 10a). Since the distribution of orientations of fibers is chosen as uniform, the overall behavior is expected to be isotropic. For a given size, numerical hydrostatic and shear creep tests are performed to estimate the average of overall responses and the relative standard deviation of these responses. As l increases, the average of nine overall responses does not change much and the standard deviation of these nine responses decreases(Fig. 10b). The standard deviation of the elastic strain of an hydrostatic test is much lower that the standard deviation of the elastic strain of a shear test, as if the representative elementary volume for the shear test were larger than the one for a hydrostatic test. Two reasons could explain such a discrepancy. On the one hand, the Poisson's ratio of the matrix, the soft phase, is twice as large as the one of the fibers. Hence the ratio of bulk moduli is lower than the ratio of shear moduli: the contrast between phases is larger if shear is considered. On the other hand, when an hydrostatic test is performed, all fibers are acting as reinforcements, independently of their orientation. In the shear test σ xy , fibers aligned along z exhibit lower stress concentrations than fibers in the xy plane. Hence, the overall result might be more sensible to the set of orientations of fibers in the cubic cell.

It should be mentioned that the relative standard deviations are lower than the error due to space discretization (5%): the later may be considered as a bias. The average overall results are stable because the grid size N was proportional to l. Relative standard deviations increase during numerical creep tests, as the contrast between tangent stiffnesses of inclusions and matrix increases. The estimate of the probability distribution function of the Von Mises stress depends slightly on l: it gets smoother as l increases (Fig. 11).

Results and discussions

Overall properties

The results of the Hashin-Shtrikman lower bound in the Laplace-Carson space are first compared to 3D numerical results featuring spherical inclusions or short cylindrical inclusions in figure 12, the volume fraction of inclusions being 20%. The length of short cylindrical inclusions is equal to their diameter. For the hydrostatic creep test, the numerical time-dependent strains are close to the one estimated by the Hashin-Shtrikman bound, while a discrepancy is to be noticed for the shear creep test. This discrepancy may be attributed to the contrast between the mechanical behavior of phases: the elastic bulk modulus of glass fibers is 8.1 times the one of PVC and the shear modulus of glass fiber is 28.3 times the one of PVC. During viscoelastic 3D computations, this contrast changes at each time step: the ratio of bulk moduli changes from 8.1 to 10.9, while the ratio of shear moduli changes from 28.3 to 36.5. For the shear creep test, the timedependent strains depends slightly on the shape of inclusions: short cylindrical inclusions induce a small decrease of the time-dependent strain compared to spherical inclusions.

For 100µm-long fibers, the aspect ratio of fibers induces a larger difference between the behavior estimated by the Hashin-Shtrikman bound and the one estimated by 3D numerical computations. The Hashin-Shtrikman bound does not take account of the shape of inclusions. Hence, the shape of inclusions is a good candidate to explain the difference between the overall strain estimated by the Hashin-Shtrikman bound and the one estimated by numerical simulations. Indeed, the Mori-Tanaka estimate, with an aspect ratio identical to the one of the fibers [START_REF] Diniz Melo | Mechanical behavior of GRP pressure pipes with addition of quartz sand filler[END_REF], is very close to the result of numerical simulations (Fig. 13).

A parametric study of the influence of the volume fraction of inclusions (Fig. 14) and of the aspect ratio (Fig. 13) of fibers is performed. The range of aspect ratio and volume fraction is limited by the Random Sequential Algorithm applied to build the unit cells. Mori-Tanaka estimates are close to numerical results, as long as the aspect ratio is lower than 10 and the volume fraction lower than 20%.

Stress concentrations

The 3D numerical computations produce estimates of stress concentrations (Fig. 15). Hence, it is possible to display the Von Mises stress within the matrix and the Rankine criteria within the fibers. It is shown here that there is little difference between the instantaneous Von Mises stress and the one 50 days after loading (Fig. 16). The Von Mises stress slightly decreases far from the fibers and increases close to the fibers. The fibers already bear much of the loading right after loading time and it increases slightly with time. It is clearly visible on the probability distribution function of the Rankine criteria during an hydrostatic creep test. Since no orientation of fibers is favored by this loading, stress concentrations in fibers trigger a rise on the probability distribution function and this rise shifts toward larger stress concentrations during the creep test. On the contrary, in case of a shear creep test, there is no rise on the probability distribution function and this distribution does not change over time elapsed since loading. A small increase of large stress concentrations in the matrix is visible on the probability distribution function of the Von Mises stress in the matrix.

In case of an hydrostatic loading, the Mori-Tanaka model expects the Rankine criteria to be equal in all fibers. Yet, the numerical computations exhibit a wide distribution of the Rankine criteria in the fibers (Fig. 16). Though the overall strains predicted by these models are similar, stress concentrations may be different.

Prediction of lifetime of the material under tensile stress

A tensile stress σ l is applied to the material and both the overall strain and the Rankine criteria in fibers are estimated. According to the Mori-Tanaka estimate of stress concentrations, the fibers aligned with the direction of loading feature the largest Rankine criterion σ R (t) and they are expected to fail first [START_REF] Hashimoto | Prediction of tensile strength of discontinuous carbon fiber/polypropylene composite with fiber orientation distribution[END_REF]. The Von Mises stress in these fibers increases with time elapsed since loading, which explains the delayed rupture of the material. Since the model remains in the range of linear viscoelasticity, a single run is necessary to define a ratio r(t) = σ R (t)/σ l between the Rankine criterion and the level of loading and this ratio increases with time. A critical stress of fibers σ c is introduced and the lifetime t l under load σ l is such that r(t l )σ l = σ c = r(0)σ 0 , where σ 0 is the tensile strength of the material. The Mori-Tanaka estimate predicts a decrease of the lifetime as the creep load increases (Fig. 17), which is consistent with experiments. The estimated strain at failure is almost uniform in time. The critical stress of fibers σ M T c is set to 2.4GPa to match the experimental instantaneous tensile strength of the material.

3D numerical computations also deliver an estimate of the Rankine criteria in the microstructure. As in [START_REF] Lavergne | Effects of grain size distribution and stress heterogeneity on yield stress of polycrystals: A numerical approach[END_REF], the failure of the specimen is expected to occur when the Rankine criterion is above σ 3D c in a volume fraction of fiber equal to c 3D R . Note that, using the maximum value of the criterion in all the specimen (c 3D R -→ 0) would make the result too volatile and using a large c 3D R would not be realistic. The volume fraction c 3D R being set, σ 3D c is adjusted to match the overall instantaneous tensile strength of the material (Fig. 18). If σ 3D c is set to 2.4GPa, a c 3D R of 0.5% is required to match the instantaneous tensile strength. Yet, for a given creep load, the numerical simulations predict a shorter lifetime than the Mori-Tanaka estimate using the same value σ M T c = 2.4GPa. Indeed, the fact that the Mori-Tanaka estimate is based on the Rankine criterion of the average stress for a given direction is a reason of this discrepancy: the 3D numerical estimate accounts for the heterogeneity of the stress field in the fibers. A larger c 3D R mitigates the effect of stress concentrations: for c 3D R = 2% and σ 3D c = 1.2GPa, the numerical estimate of creep failure is close to one obtained by the Mori-Tanaka scheme, which requires the Rankine criteria of the average stress in fiber parallel to the loading direction to be σ M T c = 2.4GPa. Hence, the Mori-Tanaka model is a practical tool to estimate the creep failure of the composite material.

Aging affects the Mori-Tanaka estimate of creep failure in the same way it changes the creep strain: the composite material is aging viscoelastic and can be modeled by the time-shift method, the µ parameter being the one identified on the creep tests of the polymer matrix. The lifetime of the pipe is expected to increase if it is loaded later than 21 days after quenching.

A parametric study of the effect of the volume fraction of fibers and fiber length on the Mori-Tanaka estimate of creep failure is performed (Fig. 19). The tensile strength of the glass fibers is set at σ M T c = 2.4GPa and the composite material is loaded at 21 days. A raise of the volume fraction of inclusions induces a raise of the tensile strength of the material and its lifetime. If the fibers are more than 1mm long, the fiber length has little influence on the tensile strength and creep failure. On the contrary, fibers smaller than 300µm may enable larger tensile strengths and improve durability. In this case, larger strains may appear and the durability might be limited by the mechanical properties of the polymer or by the quality of bonding between fibers and polymers, which are not taken account of in the present study.

Discussion

The model presented in this article has several limitations compared to other models in the literature. It has been mentioned that it does not account for the statistical variability of the strength of fibers [START_REF] Chiao | Stress-rupture of s-glass epoxy multifilament strands[END_REF][START_REF] Phoenix | Modeling the statistical lifetime of glass fiber/polymer matrix composites in tension[END_REF] and the effect of moisture on the stress corrosion cracking of glass fibers [START_REF] Barbero | Phenomenological prediction of tensile strength of e-glass composites from available aging and stress corrosion data[END_REF][START_REF] Khennane | Durability of glass polymer composites subject to stress corrosion[END_REF] is ignored. Moreover, this model does not describe the progressive rupture of the fibers until the final rupture of the specimen: models have been designed to investigate damage phenomena [START_REF] Sasayama | Prediction of failure properties of injection-molded short glass fiber-reinforced polyamide 6,6[END_REF][START_REF] Hashimoto | Prediction of tensile strength of discontinuous carbon fiber/polypropylene composite with fiber orientation distribution[END_REF], by including damage in the matrix or debonding. Regarding uniaxial composites, the shear lag model [START_REF] Hedgepeth | Stress concentrations in filamentary structures[END_REF][START_REF] Beyerlein | Shear-lag model for failure simulations of unidirectional fiber composites including matrix stiffness[END_REF] may be extended to viscous matrices [START_REF] Blassiau | Three-dimensional analysis of load transfer micromechanisms in fibre/matrix composites[END_REF][START_REF] Kotelnikova-Weiler | Kinetic of fibre ruptures in a UD composite material with a viscoelastic matrix subject to traction [Cinétique de ruptures de fibres dans un matériau composite UD soumis à la traction avec une matrice viscoélastique[END_REF] to investigate the load transfer after the breakage of a fiber and explain delayed failure of the composite under constant load. The present study remains in the range of the linear strain theory, which is not seen as a limitation, since the experimental strain at failure is less than 2% [START_REF] Guedes | On the prediction of long-term creep-failure of GRP pipes in aqueous environment[END_REF]. Numerical strain at failure reaches 7% if 100µm-long fibers are considered. In such cases, it is likely that either damage in the matrix or fiber debonding would have to be considered, as performed in [START_REF] Yang | Multiscale modeling of matrix cracking coupled with interfacial debonding in random glass fiber composites based on volume elements[END_REF].

Regarding the ability of the Mori-Tanaka estimate to depict the effect of short fibers on creep strains, our results is limited to volume fraction of fibers lower than 20%. In this range, the Mori-Tanaka estimate is expected to be close to the one of Lielen [START_REF] Lielens | Prediction of thermomechanical properties for compression moulded composites[END_REF]. Hence, our results comply with existing results in the range of elasticity [START_REF] Iii | Stiffness predictions for unidirectional short-fiber composites: Review and evaluation[END_REF][START_REF] Moussaddy | A new definition of the representative volume element in numerical homogenization problems and its application to the performance evaluation of analytical homogenization models[END_REF][START_REF] Ghossein | A comprehensive validation of analytical homogenization models: The case of ellipsoidal particles reinforced composites[END_REF].

The fiber failure may be defined by using the axial stress [START_REF] Baxevanis | A micromechanically based model for damage-enhanced creep-rupture in continuous fiber-reinforced ceramic matrix composites[END_REF][START_REF] Hashimoto | Prediction of tensile strength of discontinuous carbon fiber/polypropylene composite with fiber orientation distribution[END_REF]. It may be the case in the shear-lag model, where fibers may be considered as springs [START_REF] Beyerlein | Shear-lag model for failure simulations of unidirectional fiber composites including matrix stiffness[END_REF]. The use of the Rankine criteria should deliver comparable results: it has been observed on our 3D computational results that the maximal principal stress is almost parallel to the axis of the fiber during an uniaxial tensile creep test (Fig. 20).

Conclusion

The results of full field computations and mean-field homogenization methods have been compared in the range of viscoelasticity, on a fiber reinforced polymer. If the volume fraction of fibers is lower than 20% and the behavior of the matrix similar to the one of PVC, combining the Laplace-Carson transform on the equivalent time and the Mori-Tanaka scheme was sufficient to produce an accurate estimate of the overall time-dependent strain.

3D numerical simulations produce estimates of stress concentrations. It has been shown that stress concentrations may change during a creep test. In particular, the fibers bear an increasing part of the load: an increase of Rankine criteria in fibers is clearly visible during an hydrostatic creep test. The Von Mises stress may slightly increase in the viscoelastic matrix as well, on the high side of the Von Mises stress probability distribution function, even if the most part of the matrix is relaxing. Such a feature might trigger a delayed damage of the material. As 3D numerical computations and the Mori-Tanaka scheme lead to similar estimate of creep failure based on the Rankine criteria, the Mori-Tanaka model remains an efficient and practical tool to estimate the creep failure of the composite material.

Regarding 3D numerical simulations, further developments are necessary to properly handle the case of long fibers such as the one used in the automotive industry or pipes.

Appendix A. Hashin-Shtrikman bound and Laplace-Carson transform

Appendix A.1. The correspondence principle

The non-aging linear viscoelastic problem corresponds to elastic problems thank 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). 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) = C i ε(x, t) x ∈ inclusions σ(x, t) = C m ε(x, t) + τ C 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) = C i ε(p) x ∈ inclusions σ(p) = (1 + pτ )C 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 Hashin-Shtrikman bound is an analytical model which provides estimates of the macroscopic response of the elastic material.

< σ(p) >= Ĉ HS (p) < ε(p) >= Ĉ HS (p)E
ĈHS-(p) is isotropic and its bulk modulus and shear modulus are:

K m (p) = (1 + pτ )K m µ m (p) = (1 + pτ )µ m K HS-(p) = K m (p) + f i 1 K i -Km(p) + 3(1-f i ) 3Km(p)+4µm(p) µ HS-(p) = µ m (p) + f i 1 µ i -µm(p) + 6(Km(p)+2µm(p))(1-f i ) 5µm(p)(3Km(p)+4µm(p))
The Laplace-Carson transform of the macroscopic stress < σ > (p) is computed and the last stage is inverting this transform.

The internal variables are:

γ k (t) = t 0 t µ τ k t µ ref e -t µ+1 -t µ+1 (µ+1)t µ ref τ k C -1 k σ(t )dt and: ε(t) = C -1 0 σ(t) + k≥1 γ k γk (t) + t µ τ k t µ ref γ k (t) - µ t γ k (t) = t µ τ k t µ ref C -1 k σ(t)
The constitutive law is integrated on the time step [t i , t i+1 ] under the following assumption:

σ(t) = σ(t i ) + ∆σ t µ+1 -t µ+1 i t µ+1 i+1 -t µ+1 i
This choice for σ(t) enables to solve the evolution of internal variables through a variation of parameter. The solution for a constant stress γk (t) +

t µ τ k t µ ref γ k (t) - µ t γ k (t) = 0 writes: γ k (t) = Ct µ e - t µ+1 -t µ+1 i (µ+1)t µ ref τ k
The assumption on σ(t) is used:

Ċ = t µ τ k t µ ref µ+1 t µ+1 i+1 -t µ+1 i e t µ+1 -t µ+1 i (µ+1)t µ ref τ k C -1 k ∆σ C = µ+1 t µ+1 i+1 -t µ+1 i e t µ+1 -t µ+1 i (µ+1)t µ ref τ k C -1 k ∆σ + D Hence: γ k (t) = (µ+1)t µ t µ+1 i+1 -t µ+1 i C -1 k ∆σ + Dt µ e - t µ+1 -t µ+1 i (µ+1)t µ ref τ k γ k (t) = (µ+1)t µ t µ+1 i+1 -t µ+1 i   1 -e - t µ+1 -t µ+1 i (µ+1)t µ ref τ k   C -1 k ∆σ + γ k,t i t µ t µ i e - t µ+1 -t µ+1 i (µ+1)t µ ref τ k γ k,i+1 = (µ+1)t µ i+1 t µ+1 i+1 -t µ+1 i   1 -e - t µ+1 i+1 -t µ+1 i (µ+1)t µ ref τ k   C -1 k ∆σ + γ k,t i t µ i+1 t µ i e - t µ+1 i+1 -t µ+1 i (µ+1)t µ ref τ k
The strain must be computed:

ε(t) = C -1 0 σ(t) + k≥1 γ k (t) ∆ε = t i+1 t i C -1 0 σ(t ) + k≥1 γ k (t ) dt 18
The first part, C -1 0 σ(t ) writes:

∆ε 0 = C -1 0 ∆σ
The part coming from the internal variables writes:

∆ε k =   1 - (µ + 1)t µ ref τ k t µ+1 i+1 -t µ+1 i   1 -e - t µ+1 i+1 -t µ+1 i (µ+1)t µ ref τ k     C -1 k ∆σ+ t µ ref τ k t µ i   1 -e - t µ+1 i+1 -t µ+1 i (µ+1)t µ ref τ k   γ k,t i
The sum of these parts is ∆ε. This defines the tangent elastic problem.

Appendix C. 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 = C 0 (Σ -σ i ) + ε i , where σ i denotes the volume average of σ i . In this case, the error on boundary condition is modified as error bc =< σ i > -Σ [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF].

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 σ i (x)←C(x)ε i (x) + σ 0 (x) end for σi ←FFT(σ) Compute error on equilibrium error eq for ξ ∈ frequencies do ε(ξ) i+1 ←ε(ξ) i -Γ 0 (ξ) : σi (ξ) end for ε(0) i+1 ←E ε i+1 ←FFT -1 (ε i+1 ) end while (a) photomicrograph of a cross section of a fiber-reinforced plastic, the cross section being of a fractured tensile dumbbell [START_REF] Rahrig | Glass fiber reinforced vinyl chloride polymer products and process for their preparation[END_REF]. The distribution of the orientations of fibers does not seem isotropic and this particular sample exhibits fiber-matrix debonding. The failure mode highly depends on the coupling agent [START_REF] Rahrig | Glass fiber reinforced vinyl chloride polymer products and process for their preparation[END_REF]. (b) a 400µm-wide sample of a fiber-reinforced plastic by the RSA algorithm is shown. The volume fraction of fibers in the cubic cell is 15%. Fibers are 10µm in diameter and 100µm in length. t-t 0 (s) pure matrix c i =5% c i =15% c i =20% Figure 14: Overall responses estimated by the Mori-Tanaka scheme (dashed lines) and 3D numerical computations (solid lines) for fiber-reinforced polymers featuring different volume fractions of fibers. A shear creep test σxy = 1 (left) and a hydrostatic creep test tr(σ) = 1 (right) are performed. The Mori-Tanaka estimates match the numerical one. The length of fibers is 100µm, corresponding to an aspect ratio of 10. [START_REF] Struik | Physical aging in amorphous polymers and other materials[END_REF][START_REF] Read | Physical ageing and creep in PVC[END_REF]

[ 1 ]

 1 E. B. Rabinovitch, J. W. Summers, The effect of physical aging on properties of rigid polyvinyl chloride, Journal of Vinyl Technology 14 (3) (1992) 126-130. doi:10.1002/ vnl.730140303. URL http://dx.doi.org/10.1002/vnl.730140303

Figure 1 :

 1 Figure 1: Real and numerical microstructures.

Figure 2 :

 2 Figure 2: A rheological model made of a series of Kelvin chains

Figure 3 :

 3 Figure 3: The experimental results of Struik[START_REF] Struik | Physical aging in amorphous polymers and other materials[END_REF] on hard polyvinyl chloride (PVC) at T=20°C are digitalized. First parameters of aging µ and of the KWW function are fitted according to some of these curves (t0 = 10, 30, 100, 300 and 1000 days). Then a series of Kelvin chains is adjusted to the KWW function to enable 3D numerical computations (blue).

Figure 4 :

 4 Figure 4: The experimental results of Read et. al.[START_REF] Read | Physical ageing and creep in PVC[END_REF] on hard polyvinyl chloride (PVC) at T=23°C are digitalized. First parameters of aging µ and of the KWW function are fitted according to these curves. Then a series of Kelvin chains is adjusted to the KWW function to enable 3D numerical computations (blue).

Figure 5 :

 5 Figure 5: Left : fiber orientations are binned into 20 sectors of equal angle along y axis. Right : the estimated number of fiber in each sector n f,e = 61, 1 is compared to the observed ones n f,o (θ), θ ∈ 1..20. An isotropic distribution of orientation can produce such a set of observations.

Figure 6 :

 6 Figure 6: Relative error due to time stepping. Numerical relaxation tests are performed on the matrix for different rate of growth of the time step b = (ti+1ti)/(titi-1). The relative error decreases as the time steps decrease. Rises correspond to characteristic times of Kelvin chains. The relative error seems discrete due to numbers being stored in files as single-precision floats.

Figure 7 :

 7 Figure 7: The same microstructure of length 133µm is discretized on grids of size N = 32, N = 128 and N = 648. Hence the corresponding total numbers of voxels are respectively 32 3 , 128 3 and 648 3 . Cross sections of cylindrical fibers are ellipses, except at edges.

Figure 8 :Figure 9 :

 89 Figure8: The estimated instantaneous strain εxy and the estimated strain 50 days after loading are plotted as a functions of the grid size N . A grid size corresponds to a ratio between the size of the pixel and the diameter of the fibers.

Figure 10 :Figure 11 :Figure 12 :Figure 13 :

 10111213 Figure 10: The elastic response of unit cells of different sizes l at t0 = 10 days is computed.

Figure 15 :Figure 16 :

 1516 Figure 15: Stress concentrations in a a fiber-reinforced polymer may be estimated by 3D numerical computations. A numerical shear creep test σxy = 1 is performed and the Von Mises stress is displayed at t -t0 = 0 (left) and t -t0 = 50 days (right). It is noticeable that the Von Mises stress is highly variable in the matrix. These stress concentrations do not change much over time: the fibers already bear much of the loading at the time of loading. N = 384, l = 200µm

Figure 17 := 2 .

 172 Figure 17: The Mori-Tanaka estimate of the lifetime versus the creep load. The volume fraction of fibers is 15% and their aspect ratio is 300, corresponding to 3mm-long fibers. The material is loaded at 10 days. The strength of fiber is set at σ M T c = 2.4GPa so that the estimated instantaneous tensile strength is close to the one measured by Guedes et al.[14]. The Zhurkov model of creep failure [103] is a straight line t/τ = e U -γσ adjusted by Guedes et al. on their experimental data (τ =3600s, U=38.1MPa,γ = 0.24.

Figure 18 :

 18 Figure 18: The Mori-Tanaka estimate and numerical estimates of the lifetime versus the creep load. The volume fraction of fibers is 15% and their aspect ratio is 10, corresponding to 100µm-long fibers. The material is loaded at 10 days. The grid size is N = 384 and l = 200µm. The 3D criteria depends on parameters c 3D R and σ 3D c , which can be adjusted to match the Mori-Tanaka estimate.

Figure 19 :

 19 Figure 19: The Mori-Tanaka estimate of the lifetime versus the uniaxial traction creep load. Left: effect of the volume fraction of fibers, the length of fibers being 1mm. Right: Effect of the length of fibers, the volume fraction of fibers being 15%. If fibers are more than 1mm long, the size of the fiber has little effect on tensile strength or creep failure. Fiber lengths lower than 300µm may improve the mechanical performance of the composite, but the polymer may become the limiting factor.

Table 1 :

 1 Values fitted according to the experimental results of Struik[START_REF] Struik | Physical aging in amorphous polymers and other materials[END_REF] on rigid PVC at T=20°C and values for E-glass fibers. All Poisson's ratios of PVC are assumed to be 0.40. The aging coefficient µ is set at -0.99998.

	µ	τ ref (s)		i	τ i (s)	E i (GPa)
	-0.958 (a) Aging parameters 1		0 instantaneous 1 4.8 × 10 -6 2 4.8 × 10 -5	3.27 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

Table 2 :

 2 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.

Appendix A.2. Inverting the Laplace-Carson transform

Lots 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:

and

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. This formula is useful when the Hashin-Shtrikman analytical formula or self-consistent estimate are computed in the Laplace-Carson space.

Appendix B. Integration on the time step for the exponential algorithm

The strain rate writes:

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

Error on boundary conditions: error bc ← < e i+1 > -E Error on compatibility:

end while

Figure 20: The Rankine criteria in fibers as computed by 3D numerical computation of an uniaxial creep test σx = 1MPa. Three slices are displayed. Fibers parallel to axis x feature larger Rankine criteria. Black lines correspond to the direction of the maximal principal stress in fiber. This direction is parallel to the axis of the fiber.

N