Keywords: fractional viscoelasticity, homogenization, incremental variational approach, Effective Internal Variable method, Generalized Standard Materials, FFT-based solver viscoélasticité fractionnaire, homogénéisation, approche incrémentale variationnelle, méthode effective Internal Variable, Matériaux Standards Généralisés, solveur FFT. 1

Parvenue à la conclusion de cette thèse, je suis heureuse d'adresser mes sincères remerciements à toutes celles et ceux qui ont contribué à leur échelle à cet aboutissement.

Je remercie Renald dont j'ai bénéficié pendant trois ans de l'écoute, la patience et l'immense pédagogie. Toujours disponible, il a su m'orienter dans les moments les plus laborieux. Merci pour sa naturelle bienveillance aussi bien scientifique qu'humaine.

Je remercie Noël Lahellec et Léo Morin qui ont accepté de rapporter mes travaux, et dont j'ai grandement apprécié nos échanges. Merci à Carole Nadot-Martin d'avoir présidé mon jury, ainsi que pour l'enthousiasme scientifique qu'elle m'a transmise lors des comités de suivi. Un grand merci à Rémi Cornaggia d'avoir examiné mes travaux, qui m'a fait l'honneur de me suivre aux comités, mais aussi au quotidien au laboratoire. Je remercie Djimédo Kondo dont j'admire la passion pour la mécanique et le talent pour la transmettre.

Je tiens à remercier avec émotion la grande dame qu'était Hélène Dumontet, pour m'avoir accueillie en Master et permise de m'y sentir confiante.

Je pense ensuite au laboratoire, à commencer par le service administratif. Merci à Simona pour sa gentillesse. Je remercie chacune et chacun de mes collègues qui ont déjà réalisé ou réalisent encore leur thèse, pour nos conversations (parfois) scientifiques et nos innombrables pintes. Si l'environnement à ∂'Alembert est si épanouissant, c'est bien grâce à elles et eux. Un remerciement spécial à Samy pour son cordial soutien. J'exprime aussi ma gratitude à Ghita, pilier dans mon parcours académique.

Enfin, je remercie chaleureusement ma famille, ma mère, ma soeur, mon grand frère, mon petit frère jeune chercheur en doctorat, pour leur soutien inconditionnel.

Résumé

Cette thèse porte sur la prédiction des propriétés mécaniques effectives de matériaux hétérogènes composés de constituants viscoélastiques fractionaires, au moyen d'une approche incrémentale variationnelle. Nous appliquons la méthode Effective Internal Variable (EIV) dévelopée par Lahellec and Suquet [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF], particulièrement attrayante pour le traitement de comportements viscoélastiques [START_REF] Tressou | Contribution à l'homogénéisation des milieux viscoélastiques et introduction du couplage avec la température par extensions d'une approche incrémentale directe[END_REF]. Contrairement aux methodes d'homogénéisation communément utilisées qui reposent sur le principle de correspondance et pour lesquelles les fluctuations des champs ne sont pas accessibles, cette approche incrémentale permet de calculer les propriétés effectives dans le domaine direct au moyen des méthodes variationalles [3] et [4] qui prennent en compte les seconds-moments des champs mécaniques. La méthode EIV s'inscrit dans le cadre des Matériaux Standards Généralisés (MSG), dans lequel le comportement des matériaux dissipatifs est décrit par deux potentials thermodynamiques convexes. Nous considérons des constituants viscoélastiques fractionnaires, dont la loi constitutive est décrite par des équations différentielles linéaires avec des dérivées fractionnaires. En accord avec des observations expérimentales, ce formalisme prend en compte des effets de mémoire longue à travers la superposition de plusieurs temps caractéristiques [START_REF] Caputo | Linear Models of Dissipation whose Q is almost Frequency Independent-II[END_REF]. La distribution de ces derniers est donnée explicitement par l'expression du spectre en loi puissance. Les potentiels thermodynamiques des matériaux viscoélastiques fractionnaires sont définis en cohérence avec le cadre des MSG. Cette cohérence s'appuie sur l'interprétation rhéologique de l'élément fractionnaire comme un Maxwell généralisé [START_REF] Lion | On the thermodynamics of fractional damping elements[END_REF]. Ainsi, nous tirons parti de l'extension de la méthode EIV à plusieurs variables internes dévelopée par Tressou et al. [START_REF] Tressou | Application of the variational EIV approach to linear viscoelastic phases governed by several internal variables -Examples with the generalized Maxwell law[END_REF] afin d'homogénéiser des matériaux composites contenant des constituants viscoélastiques fractionnaires. De plus, les temps caractéristiques sont adéquatement choisis à partir de la discrétisation du spectre. Cette discrétisation est réalisée avec la procédure de Papoulia et al. [START_REF] Papoulia | Rheological representation of fractional order viscoelastic material models[END_REF], basée sur une méthode des trapèzes améliorée. Plus précisément, nous appliquons cette méthode à la fonction de Mittag-Leffler impliquée dans la définition des spectres de relaxation. Nous abordons deux problèmes hétérogènes différents au moyen de la méthode EIV. Nous considérons d'abord un composite de type matrice-inclusions sous chargement harmonique, pour lequel nous rencontrons des difficultés numériques. Nous évaluons ensuite la méthode EIV sur un polycristal de glace soumis à un essai de fluage.

The mechanical behaviour of materials is considered viscoelastic if their response due to a loading history is deferred in time. For example, the constant force due to the weight of books put on a shelf will cause the shelf to deform over time. Materials like polymers, ice, rocks, wood, steel alloys or human tissues are characterized as such. They display memory effects, i.e. the mechanical response at a given time depends on the history of loading. Real viscoelastic materials are generally heterogeneous at the micro-scale. For example ice which has a polycrystalline microstructure, or foams which have a porous microstructure. Besides the coupled viscous and elastic effects, the interactions between the constituents give rise to long-memory effects that may not be present locally [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF]. These scale-transition effects can be taken advantage of for industrial purpose. For example, the dissipative acoustic properties of a viscoelastic polymer can be enhanced with the addition of rigid inclusions, as studied in the context of submarine stealth. Therefore, either in the context of earth science, biological materials or industrial purposes, the prediction of the effective mechanical properties of composite materials with viscoelastic constituents is of interest. In order to address this problem, one has to (i) specify the local viscoelastic behaviour, and (ii) apply a homogenization method that can take into account such dissipative constituents.

Local linear viscoelastic behaviours are commonly modelled in the framework of classical viscoelasticity. The constitutive laws are thus expressed as linear differential equations with classical derivative operators. Viscoelastic materials are classified in four types [START_REF] Caputo | Linear models of dissipation in anelastic solids[END_REF], characterized by a one or two characteristic times. However, experimental observations have shown that some materials displayed long-memory effects [START_REF] Nutting | A new general law of deformation[END_REF]. Such behaviours can be accurately described in the formalism of fractional viscoelasticity, i.e. by differential equations with fractional derivatives. These operators account for the superposition of several relaxation processes. Based on physical considerations, Bagley and Torvik [START_REF] Bagley | Fractional calculus -A different approach to the analysis of viscoelastically damped structures[END_REF] successfully showed the natural ocurrence of a fractional behaviour for polymers. Even though the physical origin of such behaviour is not yet entirely understood for other types of viscoelastic materials, fractional behaviours are actually found in other dissipative phenomena such as acoustics and dieletrics [START_REF] Holm | Natural Occurrence of Fractional Derivatives in Physics[END_REF].

The homogenization of composite materials with linear viscoelastic constituents can be assessed by means of the correspondence principle, which is based on the resolution of the heterogeneous problem in the spectral domain. This approach allows to apply the mean field homogenization methods developed in the context of linear elasticity. Al-ternatively, the correspondence principle also establishes the appropriate framework for numerical homogenization methods such as FFT-based solvers. Although robust and reliable, this method can be numerically expensive. In this thesis, we thus rather focus on another category of mean field homogenization methods that rely on incremental variational approaches. More specifically, we make use of the Effective Internal Variable developed by Lahellec and Suquet [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF]. This method requires the thermodynamic formulation of the local constitutive laws in the framework of Generalized Standard Materials [START_REF] Halphen | Sur les Matériaux Standards Généralisés[END_REF]. This PhD thesis aims at applying this incremental variational approach to composite materials with fractional viscoelastic constituents, i.e. with local long-memory effects. To this end, the manuscript is organized as follows :

• Chapter 1 : Constitutive behaviours for fractional viscoelasticity

In this chapter, we start with an overview on the viscoelastic phenomena. We present the formulation of the constitutive laws of linear viscoelasticity. They can be equivalently expressed in terms of differential equations or in integral form. We also present the Generalized Standard Materials (GSM) framework and specify the definition of the free-energy and the dissipation potentials. The constitutive laws derive from these two conservative and dissipative potentials. We then specify these three equivalent forms for the four types of classical viscoelastic materials. After that, the fractional damping element is introduced, and we highlight its rheological interpretation as classical elements through the expression of the fractional derivative operator. Finally, the four types of fractional viscoelastic laws are introduced.

• Chapter 2 : Basics of homogenization and contribution to numerical homogenization Firstly, the essential concepts of linear elastic homogenization are presented. We provide the basics of mean-field approaches, followed by full-field numerical homogenization where we focus on FFT-based methods which iteratively solve the Lippman-Schwinger equation for periodic media. This method presents an inherent drawback due to the Gibbs phenomenon. Several smoothing methods were developed to address this problem. We focus on the automatic data filter [START_REF] Morin | Periodic smoothing splines for FFT-based solvers[END_REF] and investigate a possible combination of this filter on both the rigidity and the compliance tensors for tackling a challenging checkerboard microstructure.

• Chapter 3 : Incremental variational approach for composites with fractional viscoelastic constituents

We extend the presentation of mean field homogenization methods for viscoelastic constituents. They are commonly based on the correspondence principle, where the heterogeneous problem is solved as a symbolic elastic one. This approach is effective and can be successfully applied for composites with fractional visoelastic materials.

We then introduce the alternative incremental variational approach developed by Lahellec and Suquet [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF], namely the Effective Internal Variable method (EIV) for fractional viscoelastic constituents. For this purpose, we start with the formulation of the free-energy and the dissipative potentials of the fractional viscoelastic damp-ing element in the GSM framework. After that, we introduce the EIV method in the general case of several internal variables and isotropic constituents [START_REF] Tressou | Application of the variational EIV approach to linear viscoelastic phases governed by several internal variables -Examples with the generalized Maxwell law[END_REF]. Taking advantage of the exact expression of the fractional relaxation spectrum, the internal variables can be determined as a finite number with appropriate moduli. This is done by using the discretization method developed by Papoulia et al., more specifically we apply it to the Mittag-Leffler function. This chapter thus provides all the tools for the implementation of the method.

• Chapter 4 : EIV estimates of composite materials with fractional viscoelastic constituents This last chapter tackles the homogenization of composites through the EIV method.

We provide insights about our numerical implementation with Python. The numerical validation of the code is then presented. Finally, we tackle two different homogenization problems. First, we consider a fractional Zenerian matrix reinforced with spherical rigid inclusions. We rely on the Generalized Self-Consistent scheme for the resolution of the auxiliary thermoelastic problem and confront the EIV estimates with reference results obtained by [START_REF] Gallican | Homogenization estimates for the effective response of fractional viscoelastic particulate composites[END_REF]. Then, we consider polycrystalline ice subject to a creep loading, following the assumption made in [START_REF] Schapery | Linear Elastic and Viscoelastic Deformation Behaviour of Ice[END_REF] that the local basal slip system of the grains follows a fractional Maxwell behaviour. We present our results and discuss them.

Introduction

Viscoelastic phenomena are observed in a wide variety of materials such as metals, wood, rocks, biological materials or polymers. When subjected to a mechanical loading, part of their response is recoverable, while the other part is defered in time. More precisely, they exhibit the capacity to store a fraction of the mechanical energy and dissipate the rest. The mechanical state at some given time inherits the effects of all previous applied loading, this dependance to the history of loading is referred to as memory effects. For small perturbations, as considered in the context of this work, the stress and strain fields are related through linear differential equations with constant coefficents, or equivalently by convolution integrals with difference kernel. The complete description of the timedependency response is contained in the spectrum, which is not accessible by experiments [START_REF] Tschoegl | The Phenomenological Theory of Linear Viscoelastic Behavior[END_REF]. Viscoelastic behaviours are usually classified in four types [START_REF] Caputo | Linear models of dissipation in anelastic solids[END_REF] for which the spectra are identified by a unique characteristic time. However, the response of real materials displays a continuous spectrum identified by an infinite number of characteristic times that generally can not be exactly determined, but only approximated.

Nonetheless, since the early XX th century, some experimental observations have been accurately described using fractional viscoelastic laws, i.e. differential equations with derivative operators of non-integer order, allowing to take into account long-memory effects, see for example [START_REF] Eldred | Kelvin-Voigt versus fractional derivative model as constitutive relations for viscoelastic materials[END_REF] [START_REF] Henriques | Viscoelastic behavior of polymeric foams: Experiments and modeling[END_REF] [21] [START_REF] Pirk | Modeling viscoelastic damping insertion in lightweight structures with generalized Maxwell and fractional derivative models[END_REF] [START_REF] Baudis | Thorough ultrasonic rheology of soft, visco-elastic materials: Example of crosslinked Polyurethane elastomer[END_REF]. The equivalent integral form exhibits a spectrum displaying an infinite superposition of characteristic times, but for which the distribution follows a power-law function that is explicitely known. Fractional calculus was primarily introduced in viscoelastic laws by Scott-Blair [START_REF] Blair | The Classification of the Rheological Properties of Industrial Materials in the light of Power-law Relations between Stress, Strain and Time[END_REF] [25] [START_REF] Blair | Analytical and Integrative Aspects of the Stress-Strain-Time Problem[END_REF], based on the experimental analysis made by Nutting [START_REF] Nutting | A new general law of deformation[END_REF] and Gemant [START_REF] Gemant | Frictional Phenomena. VIII[END_REF] that the creep function of polymeric materials could be modeled by power-law functions in time. As reported by Stassnie [START_REF] Stiassnie | On the application of fractional calculus for the formulation of viscoelastic models[END_REF], Scott-Blair had foreseen that viscoelastic materials were not described by an addition of elastic springs and Newtonian dashpots, but an intermediate in between. It was only later though that the rigorous mathematical definition of the fractional operators for viscoelasticity was established by Caputo [START_REF] Caputo | Linear models of dissipation whose Q is almost frequency independent[END_REF] [START_REF] Caputo | A new dissipation model based on memory mechanism[END_REF]. The reader may refer to [START_REF] Mainardi | An historical perspective on fractional calculus in linear viscoelasticity[END_REF] for a historical perspective on fractional viscoelasticity, and to [START_REF] Rogosin | George William Scott Blair -the pioneer of factional calculus in rheology[END_REF] for an overview of the pioneer contribution of Scott-Blair on fractional modelling in rheology.

Several attempts were made to understand the physical mechanisms underlying the fractional nature of viscoelastic materials. Bagley and Torvik [START_REF] Bagley | A Theoretical Basis for the Application of Fractional Calculus to Viscoelasticity[END_REF] successfully showed that viscoelastic polymers without crosslinking could be modeled by a fractional constitutive law of order 1/2, that naturally emerged from Rouse's theory on the dynamics of molecular chains. In [START_REF] Heymans | Fractal rheological models and fractional differential equations for viscoelastic behavior[END_REF], Heymans and Bauwens built hierarchical models by combinations of springs and dashpots, giving rise to power-law constitutive behaviours. However, they didn't consider any physical microstructure. The idea of a relation between fractional laws and fractal or hierarchical media continued to be explored on actual materials like human tissues [START_REF] Ionescu | Viscoelasticity and fractal structure in a model of human lungs[END_REF] [START_REF] Deseri | Power-law hereditariness of hierarchical fractal bones[END_REF] or foams [START_REF] Mashayekhi | Fractional viscoelasticity in fractal and non-fractal media: Theory, experimental validation, and uncertainty analysis[END_REF] [START_REF] Somayeh Mashayekhi | A physical interpretation of fractional viscoelasticity based on the fractal structure of media: Theory and experimental validation[END_REF]. In [START_REF] Butera | A physically based connection between fractional calculus and fractal geometry[END_REF], Butera and Di Paola studied the viscous flow of a fluid in a fractal porous medium, establishing a connection between fractal geometry and fractional calculus based on geometrical and physical considerations. On the contrary, considering a 2-phase random checkerboard with fractal patterns, Ostoja-Starzewski and Zhang concluded that the fractional nature of the media was not necessarily due to the microgeometry [START_REF] Ostoja | Does a Fractal Microstructure Require a Fractional Viscoelastic Model?[END_REF]. Actually, power-law models seem to arise naturally in other phenomena with attenuation like dielectrics and acoustics. The physical origin of fractional nature is yet to be investigated, but could be related to multiple mechanisms of attenuation [START_REF] Peter | Linking multiple relaxation, power-law attenuation, and fractional wave equations[END_REF] [START_REF] Holm | Natural Occurrence of Fractional Derivatives in Physics[END_REF]. Precisely in keeping with this physical intuition, the relaxation function of any fractional viscoelastic model can be written as a continuous superposition of relaxation processes, where the repartition of the characteristic times follows a power-law function [START_REF] Tschoegl | The Phenomenological Theory of Linear Viscoelastic Behavior[END_REF].

This chapter aims at providing the fundamentals of the theory of classical viscoelasticity, as well as the basics of fractional calculus. The four main types of viscoelasticity are presented and their rheological illustration is specified. We provide the formulation of the thermodynamic potentials of viscoelastic materials, and specify them for each of the four types of classical viscoelasticity. Finally, the rheological interpretation of the fractional damping element as the superposition of classical viscoelastic assemblies is introduced, which is an important feature in the context of this thesis work.

Viscoelastic phenomena

When a viscoelastic material is subjected to a deformation or a stress loading, the flowing processes inside the body requires a finite time to respond to the loading. A very short time negligible compared with the time of the experiment is associated to a purely viscous materials. If the time is infinite, the material is referred to as purely elastic. Generally, all materials are viscoelastic. However they are identified as such if their characteristic time response is comparable with the time scale of the experiment [START_REF] Tschoegl | The Phenomenological Theory of Linear Viscoelastic Behavior[END_REF]. This section deals with the fundamental experiments allowing to characterize the long-time response of viscoelastic materials.

Transient properties : creep and relaxation

The creep and relaxation experiments consist in applying an instantaneous stress or strain loading respectively at time t = t 0 , which is then held constant. Despite the similarity of these experiments, the feasibility of their execution is not the same. The creep test is always realizable, while the relaxation test requires the material to be instantaneously deformed. Besides, they display different characteristic times, the relaxation phenomena occuring faster than the creep phenomena [START_REF] Salencon | Viscoélasticité pour le calcul de structure[END_REF].

Creep

Let's consider a uniaxial shear stress loading of amplitude σ 0 :

σ(t) = σ 0 H(t -t 0 ), (1.1)
Viscoelastic phenomena H(t) is the Heavyside function defined by : At t = t 0 , the material displays an instantaneous strain response, followed by a monotonous increase, see Fig. 1.1. The strain response is given by the relation :

H(t) = 0 if t < 0, 1 if t ≥ 0. ( 1
ε(t) = M(t)σ 0 , ( 1.3)
where M(t) is the creep function of the material.

Relaxation

Analogously, the relaxation test is carried out by applying the uniaxial strain loading of amplitude ε 0 :

ε(t) = ε 0 H(t -t 0 ). (1.4)
The stress response is composed of an instantaneous jump σ 0 at t = t 0 , then followed by a monotonous decrease, see Fig. 1.2. The stress response is governed by the relaxation function L(t) through the relation :

σ(t) = L(t)σ 0 . ( 1

Dynamic mechanical analysis

The dynamic mechanical analysis (DMA) method allows to characterize frequencydependent behaviors. For a harmonic strain (or stress) loading with pulsation ω, the stress (or strain) field exhibits an oscillating response with same pulsation, but with a delay represented by a phase δ as illustrated in Fig. 1.3, associated to the dissipation of the energy. Note that the term "dynamic" refers to the oscillating loading, and is not attributed to any inertia effects.

Let's consider the uniaxial harmonic shear strain loading of pulsation ω :

ε(t) = ε 0 sin(ωt).
(1.6)

The delayed stress response is expressed :

σ(t) = σ 0 sin(ωt + δ).
(1.7)

Figure 1.3 : Dynamic mechanical analysis. The stress response (dashed lines) is delayed by phase δ with respect to the harmonic strain loading (solid line).

Writing the fields in complex notation ε(t) = ε 0 e iωt and σ(t) = σ 0 e iωt+δ , the dynamic stiffness L * (iω) is a complex function of ω :

σ 0 ε -1 0 e iδ = L * (iω) = L ′ (ω) + iL ′′ (ω) (1.8)
The real part L ′ (ω) refers to the storage part of the behaviour, while the imaginary part L ′′ (ω) is associated to the dissipative part, and they both are positive definite. The tangeant of the phase tan(δ) referred to as the loss factor represents the dissipated energy over a period of oscillation, and is expressed as :

tan(δ) = L ′′ (ω) L ′ (ω) .
(1.9)

Non-ageing time-dependent behaviour

When a viscoelastic body is submitted to a mechanical history loading, at each time step the stress depends on the strain history. The stress can thus be expressed as a functional of the strain :

σ(t) = G    t ε(u) -∞    , (1.10)
where t denotes the present time, and u the past time from -∞ to t. Analogously, the strain can be expressed as a functional of the stress history :

ε(t) = J    t σ(u) -∞    .
(1.11)

In this work, we consider materials without aging effects, meaning their properties only depend on mechanical loadings, and not other loadings such as change of temperature, hygrometry or crystallization [START_REF] Salencon | Viscoélasticité pour le calcul de structure[END_REF]. If we consider two identical strain loading histories shifted in time by u :

ε u (τ ) = ε(τ -u), ∀τ , ∀u (1.12)
the stress responses will be identical and shifted in time by the translation u as well :

G    t ε(τ ) -∞    = G    t -u ε(τ -u) -∞    , ∀ε, ∀u. (1.13)
Naturally, the analogous reasoning also stands for a stress loading :

J    t σ(τ ) -∞    = J    t -u σ(τ -u) -∞    , ∀σ, ∀u. (1.14)

Linear constitutive law

In this section, the two equivalent representations of linear viscoelastic constitutive laws are presented, namely the integral representation resulting from Boltzmann's superposition principle, and the differential representation. The Laplace-Carson transforms are also introduced, as the expression of the constitutive laws in the Laplace-Carson domain may be convenient in the manipulation of the integro-differential equations. Furthermore, it allows to highlight the link between the two representations.

Boltzmann's superposition principle

Let's consider an arbitrary strain loading on the path u ∈ [t 0 , t], where t is the present time, and with discontinuous jumps at times t i . The instant t 0 is arbitrarily chosen such that the mechanical fields are null for t < t 0 . Thus we may fix t 0 = 0 in the following : u ∈ [0, t]. The history of the strain loading can be seen as the continuous superposition of singular impulses dε(u). We express the field in the general tensorial form :

ε(t) = t 0 H(t -u) : dε(u) + i H(t -t i ) ε i (1.15)
The stress is thus the superposition of the responses of infinitesimal relaxation tests, written as the following Stieltjes integral :

σ(t) = t 0 L(t -u) : dε(u) + i L(t -t i ) : ε i , (1.16)
that can be rewritten in the succint form :

σ(t) = t 0 L(t -u) : ε(u) du = (L ⊛ ε)(t), (1.17)
where ⊛ denotes the convolution Stieltjes product. Given that the material is at rest at t < t 0 , and applying an integration by part to (1.17), the Boltzmann's formula is obtained [START_REF] Christensen | Theory of viscoelasticity: an introduction[END_REF] :

σ(t) = L(0) : ε(t) + t 0 ε(t -u) dL(u) du : du (1.18)
Similarly, a stress loading with the same form as Eq. (1.15) results in the strain response :

ε(t) = M(0) : σ(t) + t 0 σ(t -u) : dM du (u) du (1.19)

Differential representation

A purely elastic behaviour is expressed typically with the constitutive equation of the Hookean type σ = L : ε, and a purely viscous material with Newton's law σ(t) = η : ε(t), where η is the viscous tensor. In this section, let's consider classical viscoelasticity, i.e. differential operators of integer order only. The extension to fractional operators will be presented in section 1.4.2. Therefore, the general relation between the stress and strain fields is expressed by the following differential equation :

M k=0 p k : d k σ(t) dt k = N k=0 q k : d k ε(t) dt k (1.20)
For small perturbations, the tensorial coefficients are considered constant. We define the linear differential operators :

P (D) = M k=0 p k D k , Q(D) = K k=0 q k D k , (1.21)
where the argument D denotes the operator d./ dt. In order to clarify the equivalence between this differential constitutive law and the integral representation, it is convenient to manipulate these integro-differential operators by using Laplace-Carson transforms (Appendix A). The Laplace transform L of the function f (t) is defined by :

L {f (t)} = ∞ 0 f (t)e -pt dt, (1.22)
where p is the complex Laplace variable. The Laplace-Carson transform f (p) of the function f (t) reads :

f (p) = pL {f (t)} (1.23)
The Laplace-Carson transform of Eq. (1.20) thus reads :

P (p) : σ(p) = Q(p) : ε(p) (1.24)
The same operation is applied to the integral form Eq. (1.17 with I the 4 th -order identity tensor. The viscoelastic compliance is straightforwardly obtained :

M(p) = Q(p) -1 : P (p) (1.28)
The viscoelastic stiffness L(t) and compliance M(t) tensors are thus obtained by inversion of Eqs. (1.26) and (1.28). These rational algebraic functions display singularities in the complex plane in the form of poles, but without essential singularities or branch points [START_REF] Tschoegl | The Phenomenological Theory of Linear Viscoelastic Behavior[END_REF].

Thermodynamics-based formulation of viscoelastic materials

Generalized Standard Materials

The thermodynamic formulation of a viscoelastic material in terms of a potential and a dissipation energy was firstly derived by Biot [START_REF] Biot | Theory of Stress-Strain Relations in Anisotropic Viscoelasticity and Relaxation Phenomena[END_REF]. Moreau then generalized this approach for other dissipative laws [START_REF] Jacques | Sur les lois de frottement, de plasticité et de viscosité[END_REF] (see also [START_REF] Jacques | On Unilateral Constraints, Friction and Plasticity[END_REF]), which constituted the base for introducing the Generalized Standard Materials (GSM) defined by Halphen and Nguyen [START_REF] Halphen | Sur les Matériaux Standards Généralisés[END_REF]. Under isothermal conditions, the behaviour of GSM can be entirely described by means of two potentials : the free-energy w(ε, α) and the dissipation potential φ( ε, α). They are convex functions of the state variables ε and α. In general, the latter refers to any variable that can not be recovered, here the viscous strain. The free-energy potential allows to provide the state laws defining the reversible forces :

σ = ∂w ∂ε (ε, α), A rev = ∂w ∂α (ε, α). (1.29)
The intrisic dissipation reads :

d 1 = σ irr : ε + A irr : α ≥ 0 (1.30)
where one has A irr = -A rev = A = -∂w/∂α. Moreover, we will see thereafter that some viscoelastic models are defined such that ε is a non-dissipative state variable. This leads to the nullity of the irreversible force, such that : σ irr = σσ rev = 0. In this case, the intrisic dissipation and the state laws can be reexpressed as :

d 1 = A : α ≥ 0, (1.31)
and

σ = ∂w ∂ε (ε, α), A = - ∂w ∂α (ε, α). (1.32)
These equations allow to characterize the reversible state of the system, but do not provide informations on the evolution of the dissipative variables. The dissipative potential φ( ε, α) is thus introduced in order to derive such evolution law of the internal variable :

A = ∂φ ∂ α ( ε, α), σ irr = ∂φ ∂ ε ( ε, α) (1.33)
where the differentiable dissipation potential φ must be convex with respect to ( ε, α) and finds its minimum value at ( ε = 0, α = 0). Finally, the constitutive law of GSM is given by :

σ = ∂w ∂ε (ε, α), ∂w ∂α (ε, α) + ∂φ ∂ α ( ε, α) = 0. (1.34)

Thermodynamics and dissipation potentials for linear viscoelasticity

A general formulation of classical linear viscoelastic models can be obtained in the framework of GSM, where the internal variable α refer to the viscous strain (or stress). The free-energy potential is expressed as :

w(ε, α) = 1 2 ε : L εε : ε + α : L αε : ε + 1 2 α : L αα : α (1.35)
where the 4 th -order tensors L εε , L εα and L αα are defined depending on the viscoelastic model. The dissipation potential is also a convex function of the state variables. For generality, let us consider the dual dissipation potential φ * (A), given by the Legendre transform 1 of φ with respect to the internal variable α. The dual potential is expressed as :

φ * (A) = 1 2 A : D : A (1.36)
where D is a symmetrical 4 th -order tensor.

Classical linear viscoelasticity

The general linear viscoelastic constitutive law for a loading on path u ∈ [0, t], where t is the present time, takes the form :

σ(t) = t 0 L(t -u) : dε du (u) du, (1.37)
where the relaxation function can be expressed in its general form [START_REF] Biot | Theory of Stress-Strain Relations in Anisotropic Viscoelasticity and Relaxation Phenomena[END_REF] :

L(t) = L er + δ(t)L vg + ∞ 0 G(z)e -t/z dz (1.38)
δ(t) is the Dirac distribution, the subscripts . er and . vg denote the relaxed elastic modulus and the glassy viscous modulus respectively. The term "relaxed" refers to the mechanical behaviour at long times, while "glassy" indicates the instantaneous response. Finally, G(z) is the relaxation spectrum, that represents the distribution of moduli that weights the different time ranges characterizing the relaxation processes.

This constitutive law is equivalently given by :

ε(t) = t 0 M(t -u) : dσ du (u) du, (1.39)
where the creep function takes the form :

1 The Legendre transform f * (x * ) of the differentiable and convex function f (x) is given by the convex

function f * (x * ) = xx * -f (x), with x * = df(x) dx
1 Constitutive behaviours for fractional viscoelasticity 25

M(t) = M eg + tM vr + ∞ 0 J(z) 1 -e -t/z dz (1.40)
The moduli M eg and M vr respectively indicate the glassy elastic and the viscous relaxed compliances. Dually to G(z), the quantity J(z) represents the retardation spectrum.

The four types of classical viscoelastic models

Linear viscoelastic models can be classified in four types [START_REF] Caputo | Linear models of dissipation in anelastic solids[END_REF] in accordance with the corresponding behaviour at short and long times. In classical viscoelasticity, the asymptotic states can be either purely elastic or purely viscous. The rheological one dimensional representation of the four models consists in assemblies of springs and dashpot, providing intuition on the construction of the constitutive laws, see

1 η 0 0 Kelvin-Voigt 0 0 E η Zener 1 K + E 0 K 0 Anti-Zener 0 1 η + 1 η 2 0 η 2
It is noted that the rheological assembly of the Zener model, illustrated in Table 1.1, can also be represented as a Kelvin-Voigt model with a spring added in series. Analogously, the anti-Zener model is equivalently represented with a Maxwell unit in parallel with a dashpot. The constitutive laws in the general three dimensional case are presented in the following.

Maxwell model

The constitutive law of the Mawxellian constituent reads :

M vr : σ(t) + M eg : σ(t) = ε(t) (1.41)
The relaxation and the retardance functions read :

L(t) = L eg e -t/τσ and M(t) = M eg + M vr t (1.42)
The Maxwell constituent exhibits an instantaneous elastic response, and a purely viscous behaviour at long times, with an elastic-to-viscous transient regime. Its viscoelastic stiffness tensor L(t) displays an attenuating transient regime from one to another asymptotic state characterized by relaxation times which are the eigenvalues of the tensor M eg : M -1 vr .

These constitutive laws can be derived from the free-energy and the dissipation potentials, defined by :

         w(ε, α) = 1 2 (ε -α) : L eg : (ε -α) φ( α) = 1 2 α : M -1 vr : α (1.43)

Kelvin-Voigt model

The constitutive law of the Kelvin-Voigt constituent follows the linear differential equation :

σ(t) = L er : ε(t) + L vg : ε(t) (1.44)
and the viscoelastic stiffness and compliance tensors are respectively given by :

L(t) = L er + L vg δ(t) and M(t) = M er 1 -e -t/τε (1.45)
In contrast to the asymptotic behaviours of the Maxwell constituent, the Kelvin-Voigt displays an instantaneous purely viscous response, and an elastic behaviour at long times. The relaxation times are null, while the retardation times are given by the eigenvalues of the tensor L vg : L -1 er . The above constitutive laws derive from the following potentials :

         w(ε) = 1 2 ε : L er : ε φ( ε) = 1 2 ε : L vg : ε (1.46)
Note that the only state variable is the strain, such that the dissipation potential only depends on the strain rate ε.

1 Constitutive behaviours for fractional viscoelasticity

27

Zener model

The Zener's constitutive differential law reads : (1.47) and its viscoelastic stiffness and compliance tensors are :

σ(t) + L v : L eg -L er -1 : σ(t) = L er : ε(t) + L eg : L v : L eg -L er -1 : ε(t),
L(t) = L er + Ge -t/τσ and M(t) = M eg + J 1 -e -t/τε (1.48)
As reported in Table 1.1, the Zener constituent exhibits an elastic response at both short and long times. The viscous effects are displayed in the transient regime by the viscous tensor L v . Thus the behaviour of the Zener constituent is identified by a transition from the instantaneous elastic state to the viscoelastic transient regime modulated by the characteristic time τ σ -like the Maxwell constituentfollowed by a transition from the transient regime to the relaxed elastic state of characteristic time τ ε -similarly to the Kelvin-Voigt constituent. The relaxation and retardation spectra are respectively defined by the quantities G = L eg -L er and J = M er -M eg . The corresponding potentials are given by :

         w(ε) = 1 2 (ε -α) : M -1 eg : (ε -α) + 1 2 α : (L er -M -1 eg ) : α φ( α) = 1 2 α : L vg : α (1.49)

Anti-Zener model

The constitutive differential law of the classical anti-Zener constituent reads :

σ(t) + L vr -L vg : M e : σ(t) = L vr : ε(t) + L vr -L vg : M e : L vg : ε(t) (1.50)
The stiffness and compliance tensors are given by :

L(t) = L vg δ(t) + Ge -t/τσ and M(t) = M vr t + J 1 -e -t/τε (1.51)
The asymptotic responses at both short and long times are purely viscous, with a viscoelastic transient regime in between. The conservative part is characterized by the elastic tensor M e .

In contrast to the Zener model, anti-Zener constituents exhibit a viscous-to-viscoelastic and a viscoelastic-to-viscous transitions. The retardance times are given by the eigenvalues of the tensor L vr -L vg : M e : L vg : (L vr ) -1 , and the relaxation times, by the eigenvalues of L vr -L vg : M e . The corresponding potentials are expressed as :

         w(ε) = 1 2 ε : M -1 e : ε φ( ε, α) = 1 2 ε : L vg : ε + 1 2 α : M -1 vr : α (1.52)

Generalized models

As mentioned earlier, most viscoelastic materials exhibit a behaviour with multiple relaxation or retardation times. By extending the above models, they can be derived into two forms : the generalized Kelvin-Voigt and the generalized Maxwell constituents, illustrated in Figs 1.4 and 1.5. The rheological representation of the generalized Maxwell model consists in the superposition of K Maxwell branches in parallel, while the generalized Kelvin-Voigt model is interpreted as the superposition of K Kelvin-Voigt units in series. They display the same type of asymptotic behaviours at short and long times as their 1-unit description.

Generalized Maxwell constituent

Each Maxwellian unit j straighforwardly follows the constitutive law Eq. (1.87) :

M j v : σ j (t) + M j e : σj (t) = εj (t) , ∀j ∈ [1; K] (1.53)
The rheological representation Fig. 1.4 highlights that the total strain doesn't depend on the branch : ε(t) = εj (t) ∀j ∈ [1; K]. By linearity, the superposition principle allows to determine the relaxation function. Considering the uniaxial strain loading ε(t) = ε 0 H(t), the constitutive law can be written :

σ(t) = K j=1
L j e e -t/τ j σ : ε 0 , (1.54) where the characteristic times τ j σ are the eigenvalues of the tensors M j e : (M j v ) -1 . The viscoelastic stiffness tensor is directly deduced :

L(t) = K j=1
L j e e -t/τ j σ .

(1.55)

The elastic glassy modulus is given by

L eg = L(t = 0) = K j=1 L j e .
The corresponding potentials are directly obtained from the superposition of the potentials of each individual branch :

1 Constitutive behaviours for fractional viscoelasticity

29                w(ε, α 1 , α 2 , ...) = K j=1 1 2 (ε -α j ) : L j e : (ε -α j ) φ( α1 , α2 , ..) = K j=1 1 2 αj : M j v -1
: αj

(1.56)

Generalized Kelvin-Voigt constituent

Analogously, the constitutive law of the generalized Kelvin-Voigt model is given by the superposition of the K evolution laws :

σ j (t) = L j er : ε j (t) + L j vg : εj (t) , ∀j ∈ [1; K] (1.57)
Similarly using the superposition principle, the retardation compliance tensor can be expressed by :

M(t) = K j=1 M j e 1 -e -t/τ j ε , (1.58)
with the retardation times being the eigenvalues of the tensors L j vg : (L j er ) -1 , and the elastic relaxed compliance is given by M er = K j=1

M j e . The corresponding potentials are significantly different from those of the generalized Maxwell. Indeed, the only state variable is the strain, such that the dissipation potential only depends on the strain rate :

               w(ε) = K j=1 1 2 ε : L j er : ε φ( ε) = K j=1 1 2 ε : L j vg : ε (1.59)

Collocation method

Viscoelastic materials are typically described using the Generalized Maxwell model, where the material parameters are obtained through either relaxation experiments or DMA. The fitting of the experimental data is usually done by means of the collocation method, developed by Schapery [START_REF] Schapery | Methods of transform inversion for viscoelastic stress analysis[END_REF]. This method allows to identify the corresponding moduli, requiring an a priori choice of the relaxation times. The finite series of characteristic times is usually taken as a logarithmic distribution on a chosen interval.

Let's consider the viscoelastic stiffness tensor :

Classical linear viscoelasticity

L(t) = ∞ 0 G(z)e -t/z dz (1.60)
Discretizing the spectrum in a sum of Dirac, G(z) ≈ j G j δ(z -z j ), the relaxation function is approximated by the expression :

L * (t) = K j=1
G j e -t/z j (1.61)

In order to determine the coefficients G j , and given an a priori choice of the parameters K and the characteristic times z j , the method consists in minimizing the quadratic error :

E = +∞ 0 [L(t) -L * (t)] 2 dt (1.62)
This minimization with respect of the coefficients G j (∂E/∂G = 0) implies the following condition in the Laplace-Carson domain :

L(p j ) = L * (p j ) , ∀p j = 1 z j , i = 1, ..., K (1.63)
The exact relaxation function (or experimental data) and its approximation are imposed to be equal at some specific collocation points p j in the Laplace-Carson domain, which correspond to the inverse characteristic times in the direct domain. Note that, following Rekik and Brenner [START_REF] Rekik | Optimization of the collocation inversion method for the linear viscoelastic homogenization[END_REF], the method can be improved by also minimizing the error with respect to the characteristic times z j (∂E/∂p = 0), adding the condition :

d L(p) dp p=p j = d L * (p) dp p=p j , ∀p j = 1 z j , i = 1, ..., K (1.64)
Eqs. (1.63) and (1.64) constitute the following system to solve :

L(p j ) = K k=1 G k 1 p k + p j and d L(p) dp p=p j = - K j=1 G k 1 (p k + p j ) 2 (1.65)
Three constraints must be introduced in order to solve the system (1.65). First, the coefficients G j must be positive definite. The collocation times must be in a chosen interval [z min , z max ]. Finally, at t = 0, the viscoelastic modulus must be equal to the instantaneous elastic modulus. The minimization constrained problem to be solved is thus expressed by :

               min G j ,p j   K j=1 [ L(p j ) - L * (p j )] 2 + K j=1 d L(p) dp p=p j - d L(p) dp p=p j 2   G j ≥ 0; z j ∈ [z min , z max ] and K j=1 G j = L eg (1.66)
1 Constitutive behaviours for fractional viscoelasticity

31

The system can similarly be derived for the viscoelastic compliance tensor, but without the constraint at t = 0. This improved collocation method provides accurate approximations for a low number of parameters.

As mentionned earlier, the collocation times are determined empirically, without theoretical criterion. In [START_REF] Jalocha | Revisiting the identification of generalized Maxwell models from experimental results[END_REF], Jalocha et. al. proposed a method to improve this choice depending on the conducted experiment. This strategy for DMA provided successfull results in [START_REF] Suarez-Afanador | Effective thermo-viscoelastic behavior of short fiber reinforced thermorheologically simple polymers: An application to high temperature fiber reinforced additive manufacturing[END_REF] for the identification of the effective properties of a polymer composite. The relaxation times are obtained from the inverse values of the zeros of the relaxation functions in the complex plane. Despite this improvement, the collocation method displays inherent issues. The inversion of the matrix required in the resolution of the linear system can generate negative spectrum moduli. It is thus worth mentionning an alternative algorithm proposed by [START_REF] Emri | Determination of mechanical spectra from experimental responses[END_REF], where the characteristic times are optimally predetermined and the moduli are iteratively found by minimization of the quadratic error in the direct domain.

Fractional viscoelasticity

The introduction of fractional formalism into viscoelastic constitutive laws allows to take into account long-memory effects with a low number of parameters. As mentioned in the introduction, Bagley and Torvik [START_REF] Bagley | Fractional calculus -A different approach to the analysis of viscoelastically damped structures[END_REF] made a significant step in the physical justification of such behaviours for polymers without crosslinking. The origin of the fractional nature of viscoelastic materials is still investigated, for example in [START_REF] Butera | A physically based connection between fractional calculus and fractal geometry[END_REF], [START_REF] Ostoja | Does a Fractal Microstructure Require a Fractional Viscoelastic Model?[END_REF] or [START_REF] Mashayekhi | Fractional viscoelasticity in fractal and non-fractal media: Theory, experimental validation, and uncertainty analysis[END_REF]. Nevertheless, as pointed out by Holm et al. [START_REF] Holm | Natural Occurrence of Fractional Derivatives in Physics[END_REF], this formalism actually occurs naturally in several physical dissipative phenomena.

The differential equations describing viscoelastic behaviours imply derivative operators with non-integer order γ such that 0 < γ < 1. Compared to classical assemblies, it corresponds to the substitution of the classical dashpot with a fractional damping element, for which the stress is proportional to the fractional derivative of the strain field [START_REF] Koeller | Applications of Fractional Calculus to the Theory of Viscoelasticity[END_REF] (Scott-Blair model [START_REF] Scott Blair | The Estimation of Firmness in Soft Materials[END_REF]). In the equivalent integral form, the relaxation or retardance functions exhibit a power-law spectrum that is explicitely known [START_REF] Mainardi | Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models[END_REF], and thus doesn't require to be approximated. The first part of this section deals with the definition of the fractional integro-differential operators. The Scott-Blair model is then presented and introduced into the four types of viscoelasticity.

Essentials of fractional integro-differential operators

Let's consider causal functions f (t) (all outputs depend on past inputs) that can be either real or complex, with real argument t > 0. In the following, the fractional integral operator is defined, from which the expression of the fractional derivative can be obtained [START_REF] Mainardi | Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models[END_REF].

Fractional integral operator

The Cauchy's formula for repeated integration of order n reads :

0 I n t f (t) := 1 (n -1)! t 0 (t -z) n-1 f (z) dz , n ∈ N (1.67)
and vanishes at t = 0. Eq. (1.67) can be generalized to the real positive non-integer order γ by introducing the Gamma function Γ(n) = (n -1)! . In doing so, the Riemann-Liouville fractional integral is defined :

0 I γ t f (t) := 1 Γ(γ) t 0 (t -z) γ-1 f (z) dz , t > 0 and γ ∈ R + (1.68)
The expression of the Gamma function for non integers reads :

Γ(z) := ∞ 0 e -u u z-1 du , Re(z) > 0 (1.69)

Fractional derivative operator

The derivative operator of order n is noted D n = d n / dt n , and the symbol • denotes the composition between operators. The fractional derivative of order γ can be defined by applying m -γ times the fractional integral operator (1.68), then taking m times the derivative, with m ∈ N and such that m -1 < γ < m. This composition defines the Caputo fractional derivative [START_REF] Caputo | Linear Models of Dissipation whose Q is almost Frequency Independent-II[END_REF] :

D γ f (t) := 0 I m-γ t • D m t f (t) , with m -1 < γ < m (1.70)
yielding the expression :

D γ f (t) :=            1 Γ(m -γ) t 0 1 (t -z) γ+1-m f (m) (z) dz m -1 < γ < m d m dt m f (t) γ = m (1.71)
It is noted that the two possible orders of composition are not equivalent :

D m • 0 I m-γ t ̸ = 0 I m-γ t • D m (1.72)
The other composition defines the fractional Riemann-Liouville derivative. The fractional Caputo derivative actually brings a sort of regularization to its Riemann-Liouville counterpart at t = 0. Based on arguments on the initial conditions in the Laplace domain, Mainardi [START_REF] Mainardi | Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models[END_REF] justifies that the Caputo fractional derivative is fully compatible with the classical viscoelastic constitutive laws, therefore he requires to keep definition (1.71) for viscoelastic laws.

Fractional viscoelastic constitutive law

Scott-Blair model

Based on the empirical power-law functions of Nutting [START_REF] Nutting | A new general law of deformation[END_REF] and Gemant [START_REF] Gemant | Frictional Phenomena. VIII[END_REF] for describing the creep function of polymeric materials, Scott-Blair introduced the fractional derivative operator D γ into viscoelastic constitutive laws [START_REF] Scott Blair | The Estimation of Firmness in Soft Materials[END_REF] with order 0 < γ < 1, thus corresponding to m = 1 in Eq. (1.71). The Scott-Blair model reads :

σ(t) = L f : D γ ε(t) (1.73)
with L f the viscoelastic stifness tensor of dimension Pa.s γ . When γ = 0, Eq. (1.73) reduces to a purely elastic material, while the case γ = 1 corresponds to a Newtonian viscous fluid. The Scott-Blair model thus describes an intermediate mechanical behavior between elasticity and viscosity. The long-memory effect is clearly seen by expliciting the fractional derivative operator in the constitutive law :

σ(t) = L f : 1 Γ(1 -γ) t 0 1 (t -z) γ dε(z) dz dz (1.74)
Contrary to the classical Newtonian law where σ(t) depends on the instantaneous rate of deformation loading, the Scott-Blair model refers to a material for which the stress field at the present time t is related to the whole history of loading.

The equivalent integral form of Eq. (1.73) can be assessed by making use of the Laplace transforms [START_REF] Tschoegl | The Phenomenological Theory of Linear Viscoelastic Behavior[END_REF]. In this domain, the Scott-Blair model is written :

σ(p) = L f : p γ ε(p) (1.75)
The relaxation function in the Laplace domain reads :

L(p) = L f p γ (1.76)
Using the identities in Table A.1, and with the relation

L(p) = M(p) -1
, the fractional relaxation and the retardance functions can be expressed in the direct domain :

L(t) = L f 1 Γ(1 -γ) 1 t γ (1.77) M(t) = M f 1 Γ(1 + γ) t γ (1.78)
The expressions (1.77) and (1.78) highlight the connection between power-law functions and the fractional derivative. Using the definition of the Gamma function (1.69), and applying the change of variable u = t/z, they can be put in their integral form :

           L(t) = ∞ 0 G(z)e -t/z dz G(z) = L f 1 Γ(1 -γ)Γ(γ) 1 z 1+γ (1.79) and            M(t) = ∞ 0 J(z) 1 -e -t/z dz J(z) = M f 1 Γ(1 + γ)Γ(1 -γ) γ z 1-γ (1.80)
It must be emphasized that Eqs. (1.79) and (1.80) are similar to the continuous form of the classical generalized models (1.58) and (1.55). The fractional relaxation and creep functions can thus be respectively interpreted as a continuous superposition of Maxwell units in parallel and Kelvin-Voigt units in series. However, the expression of the weighting functions G(z) and J(z) are explicitely known, the repartition of the characteristic times following a power-law function.

The Scott-Blair element is usually illustrated with a diamond [START_REF] Ronald | A Generalized Derivative Model for an Elastomer Damper[END_REF] (see Fig. 1.6) and is referred to as the spring-pot or fractional dashpot [START_REF] Koeller | Applications of Fractional Calculus to the Theory of Viscoelasticity[END_REF]. Unable to reproduce itself the asymptotic experimental observed states, this element is incorporated into the four types of viscoelastic models in place of the classical dashpot. The four types of fractional viscoelastic models are derived in the next section.

Fractional viscoelastic constitutive law

The fractional viscoelastic constitutive law for a loading path u ∈ [0, t] reads :

σ(t) = t 0 L(t -u) : dε(u) du du (1.81)
where L(t) is the fractional viscoelastic stiffness tensor :

L(t) = L er + L fg t γ Γ(1 -γ) + ∞ 0 G(z)E γ [-(t/z) γ ] dz (1.82)
1 Constitutive behaviours for fractional viscoelasticity

35

where L er is the elastic relaxed modulus, L fg is the glassy fractional viscous tensor, and E γ (z) is the Mittag-Leffler function defined by :

E γ (z) = ∞ k=0 z k Γ(1 + γk) , γ > 0, z ∈ C (1.83)
where the gamma function Γ is an extension of the factorial function for complex numbers.

The classical exponential function is thus a particular case of the Mittag-Leffler function.

The constitutive law is dually expressed by :

ε(t) = ∞ 0 M(t -u) : dσ(u) du du (1.84)
with the fractional compliance viscoelastic tensor :

M(t) = M eg + M fr t γ Γ(1 + γ) + ∞ 0 J(z) (1 -E γ [-(t/z) γ ]) dz (1.85)
where M fr is the fractional viscous tensor at long times.

The equivalent fractional differential form reads :

I + N k=1 p k D γ k : σ(t) = M + M k=1 q k D γ k : ε(t), γ k = k + γ -1 (1. 86
)
where I is the 4 th -order identity tensor. The equivalence between Eqs. (1.81) or (1.84), and Eq. (1.86) can be derived in the same way as in section 1.2.2. Still following [START_REF] Caputo | Linear models of dissipation in anelastic solids[END_REF], the four types of fractional viscoelastic materials are presented in the following.

The four types of fractional viscoelasticity

The extension of the four classical models in fractional viscoelasticity doesn't modify the elastic asymptotic behaviours. However, the viscous asymptotic states and the transient regimes are governed by the fractional viscoelastic contribution. The asymptotic states for the one-dimensional representation of the fractional models are reported in Table 1.2. Similarly to the classical models, the fractional Zener constituent can be equivalently represented with a fractional Kelvin-Voigt in series with a spring. The fractional anti-Zener can also be illustrated with a fractional Maxwell constituent with a fractional dashpot in parallel. In the following, the fractional types of viscoelasticity are compared with the classical models. It is specified that the Mittag-Leffler function can be implemented following the algorithm 4 (see Appendix B) from [START_REF] Gorenflo | Computation of the Mittag-Leffler function[END_REF]. Fractional Maxwell 1 E

1 E 1 τ γ 0 0 Fractional Kelvin-Voigt 0 0 E E 1 τ γ Fractional Zener 1 K + E 0 K 0 Fractional Anti-Zener 0 1 E 1 τ γ 1 + 1 E 2 τ γ 2 0 E 2 τ γ 2

Fractional Maxwell model

The fractional Mawxellian constituent follows the fractional differential equation :

M fr : σ(t) + M eg : D γ σ(t) = D γ ε(t) (1.87)
The relaxation and the retardance functions read :

     L(t) = L eg E γ [-(t/τ σ ) γ ] M(t) = M eg + M fr t γ (1.88)
Similarly to the classical model, the fractional Maxwell constituent displays an instantaneous elastic modulus L eg , but then exhibits a fractional viscous transient regime, followed by a fractional viscous behaviour at long times. The comparison between the classical retardance function (γ = 1) and the fractional one is drawn in Fig.

Fractional Kelvin-Voigt model

The constitutive law of the fractional Kelvin-Voigt reads :

σ(t) = L er : ε(t) + L fg : D γ ε(t) (1.89)
with :

     L(t) = L er + L fg t -γ Γ(1 -γ) M(t) = M er (1 -E γ [-(t/τ ε ) γ ]) (1.90)
The classical relaxation function of the Kelvin-Voigt constituent displays a singular viscous response at t = 0, which is not reported in Fig. 1.10. The fractional model also exhibits an instantaneous (fractional) viscous response, but it is followed by a smooth transient regime, highlighting the long memory effects features of the fractional operator.

Fractional Zener model

The fractional Zener's constitutive differential law reads : (1.91) and its viscoelastic stiffness and compliance tensors :

σ(t)+L f : L eg -L er -1 : D γ σ(t) = L er : ε(t)+L eg : L f : L eg -L er -1 : D γ ε(t),
     L(t) = L er + GE γ [-(t/τ σ ) γ ] M(t) = M eg + J (1 -E γ [-(t/τ ε ) γ ]) (1.92)
The fractional Zener constituents exhibits the same asymptotic elastic states as the classical Zener model, but displays a different transient fractional viscoelastic regime, see Figs.

Fractional anti-Zener model

The differential law of the classical fractional anti-Zener constituent reads :

σ(t) + L fr -L fg : M e : D γ σ(t) = L fr : D γ ε(t) + L fr -L fg : M e : L fg : D 2γ ε(t) (1.93) and :          L(t) = L fg t -γ Γ(1 -γ) + GE γ [-(t/τ σ ) γ ] M(t) = M fr t γ Γ(1 + γ) + J (1 -E γ [-(t/τ ε ) γ ]) (1.94)
1 Constitutive behaviours for fractional viscoelasticity

39

The asymptotic states display purely fractional viscoelastic behaviours, while the transient regime is governed by the "transient" elastic tensor M e .

t

M (t) = 0.2 = 0.5 = 0.7 = 1

Conclusion

In this chapter, we provided the main features of viscoelasticity phenomena, as well as the fundamentals of the theory of classical viscoelasticity [START_REF] Tschoegl | The Phenomenological Theory of Linear Viscoelastic Behavior[END_REF]. The four types of classical viscoelasticity [START_REF] Caputo | Linear models of dissipation in anelastic solids[END_REF] are generally not accurate enough to describe materials displaying longmemory effects. Use can thus be made of generalized rheological models, but the spectrum can only be approximated [47] [48] [49] [START_REF] Emri | Determination of mechanical spectra from experimental responses[END_REF]. Some materials have been empirically described by means of fractional constitutive laws with success [START_REF] Henriques | Viscoelastic behavior of polymeric foams: Experiments and modeling[END_REF] [START_REF] Baudis | Thorough ultrasonic rheology of soft, visco-elastic materials: Example of crosslinked Polyurethane elastomer[END_REF]. They exhibit power-law spectra that can be exactly expressed [START_REF] Mainardi | Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models[END_REF]. The relation between power law spectra and fractional calculus can be highlighted by means of Laplace transforms (or Fourier transforms). It must be also pointed out that fractional viscoelastic behaviours could be rheologically interpreted as an infinite continuous superposition of classical rheological elements. Finally, let us indicate that, based on their rheological interpretation, fractional viscoelastic models can be also put in the framework of Generalized Standard Materials. This will be explicited in chapter 3.

Introduction

Homogenization of composite materials can generally be addressed with two main types of methods. Mean field methods rely on a statistical description of the microstructure. They provide analytical or semi-analytical estimates, improving the understanding of links between the local microstructure and the properties at the macroscopic scale [START_REF] Milton | The Theory of Composites[END_REF]. Homogenization schemes were initially developed in linear elasticity, from the pioneer work of Eshelby on an inclusion embedded in an infinite media [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF]. His key idea was the base for tackling other types of microstructures, such as composite sphere assemblages [START_REF] Hashin | On some Variational Principles in Anisotropic and Nonhomogeneous Elasticity[END_REF], matrix-inclusions [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF], polycrystals [START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF]... Otherwise, one maye use full-field computations which are based on robust numerical methods, for example the Finite Element Method (FEM) [START_REF] Yvonnet | Computational Homogenization of Heterogeneous Materials With Finite Elements[END_REF]. This method involves the generation of meshes and the inversion of large matrices, causing high computational costs for complex microstructures. Alternatively, Moulinec and Suquet [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] developed a computational homogenization method based on Fast-Fourier Transforms (FFT), which operates on regular grids of voxels, allowing the use of digital images of microstructures. This approach is built upon the iterative integration of the Lippman-Schwinger equation solved in the Fourier domain. Initially based on a fixed-point algorithm, this approach has been extend by several authors who proposed other schemes to improve its convergence (see for instance [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF], [START_REF] Michel | A Computational Method Based on Augmented Lagrangians and Fast Fourier Transforms for Composites with High Contrast[END_REF], [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF] and [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF]). The method was originally proposed for mechanical loadings, but it was also applicable to a wide variety of problems, such as conductivity [START_REF] Willot | Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF], diffusivity [START_REF] To | An FFT method for the computation of thermal diffusivity of porous periodic media[END_REF], coupled problems like piezoelectricity [START_REF] Brenner | Numerical computation of the response of piezoelectric composites using Fourier transform[END_REF] or damage and fracture mechanics [START_REF] Chen | Fast Fourier transform Solver for Damage Modeling of Composite Materials[END_REF] [START_REF] Zhu | An incremental-iterative method for modeling damage evolution in voxel-based microstructure models[END_REF]. The interested reader may refer to the recent review by Schneider [START_REF] Schneider | A review of nonlinear FFT-based computational homogenization methods[END_REF].

However, regardless of the nature of the considered problem, this method presents inherent drawbacks. The discontinuities present in the microstructure lead to the Gibbs phenomenon, resulting in spurious oscillations of the mechanical fields. Several methods were developed in order to circumvent this limitation. One strategy consists in modifying the discrete Green's operator, classically obtained from the continuous pseudo-spectral definition. Instead, these methods apply a finite difference discretization of partial derivatives, leading to a smooth Green's opeator and making it consistent with the underlying voxel grid [START_REF] Willot | Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF]. Besides removing the spurious oscillations, this method generally improves the convergence rate of the algorithm. Another approach consists in considering composite voxels on the interfaces, and homogenize them in order to define a new smooth microstructure [START_REF] Gélébart | Filtering material properties to improve FFTbased methods for numerical homogenization[END_REF]. This method requires the identification of the interfaces and the definition of the normal vector for each composite voxel. Although efficient for removing the spurious oscillations, it adds complexity to the computation. It is also worth mentioning the recent method proposed by [START_REF] Zecevic | New large-strain FFT-based formulation and its application to model strain localization in nanometallic laminates and other strongly anisotropic crystalline materials[END_REF] for tackling large-strain problems with FFT-based solvers. Three grids of material points are considered : (i) a reference regular grid, (ii) an initial irregular grid representing a stress-free displacement field applied to the first grid, and (iii) a current irregular grid which allows to consider a better discretization of curved interfaces. Coupled to the modified discrete Green operator procedure, this approach allows to both consider large-strain problems and remove spurious oscillations.

An alternative approach was proposed by Morin et. al [START_REF] Morin | Periodic smoothing splines for FFT-based solvers[END_REF]. These authors developed a periodic data filtering method that automatically smoothes discontinuities, where the width of the interface is modulated by a unique parameter. It does not modify the FFT-based algorithm, but only constitutes an additional step applied to the material properties. In this chapter, after providing the basics of mean field homogenization in linear elasticity, we then introduce the FFT-based algorithm. Finally, the latter mentioned smoothing filter is presented and challenged on a complex microstructure for which we then investigate a weighted combination of this filter.

Basics of homogenization in linear elasticity

The micromechanical approach relies on an upscaling procedure which provides links between the heterogeneous microstructure and the effective properties. The principle is to replace a heterogeneous material with a fictitious one -the Homogeneous Equivalent Media (HEM) -with homogeneous properties that are equivalent to that of the real material at the macroscopic scale. Given the characteristic size d of the heterogeneities, the characteristic scale ℓ of the microstructure must be carefully chosen. Denoting L the dimension of the structure, the following conditions of scale separation must be respected :

-ℓ ≪ L : this condition ensures that the structure is described as a continous medium.

-ℓ ≫ d : this allows to assign the homogeneous properties of the HEM on each point of the structure.

This separation of scale defines the Representative Volume Element (RVE) that statistically represents the heterogeneous material. The RVE will be homogenized, and it can be considered as a material point of the macroscopic material. The micromechanical approach in the framework of linear elasticity relies on three steps :

(i) The representation aims at describing the spatial distribution of the constituents and their mechanical behaviour; (ii) The objective of the localization step is then to establish the relation between the local and the macroscopic stress and strain fields; (iii) Finally, the homogenization procedure determines the effective properties of the HEM.

Micromechanical approach

Representation

The Representative Volume Element of the microstructure must be chosen such that it is statistically uniform (invariant by translation) besides ensuring the separation of scale L ≫ ℓ ≫ d. For a periodic microstructure, this choice would straightforwardly reduce to the periodic element of volume. Let us condider a general random material, the domain of the RVE is denoted Ω and is composed of N phases of volume Ω (r) , r = 1, ..., N . The spatial distribution of phase r is defined by the the indicator function χ (r) (x) :

χ (r) (x) = 1 if x ∈ Ω (r) 0 otherwise. ∀x ∈ Ω (2.1)
Assuming that the constituents are uniform with stiffness L (r) and perfectly bonded on the interfaces, the heterogeneous stiffness tensor is defined as :

L(x) = N r=1 L (r) χ (r) (x) (2.2)
As mentioned before, the microstructure of a periodic media is entirely defined by the unit cell. Regarding random composite materials, the characterization of the microstructure requires a statistical approach. The indicator function indeed depends on the realization α ∈ S, where the space S contains all the micro-configurations of the material. The probability to find phase r for a fixed x is associated to the probability density of the realization p(α) such that :

p (r) (x) := S χ (r) (x, α)p(α) dα, (2.3)
and p (rs) (x, x ′ ) is the probability to simultaneously find phase r at x and phase s at x ′ :

p (rs) (x, x ′ ) := S χ (r) (x, α)χ (s) (x ′ , α)p(α) dα. (2.4)
As the microstructure is statistically uniform, the probability (2.3) is constant, and (2.4) depends only on h = xx ′ . Moreover, the considered composites are assumed to be statistically isotropic, thus p (rs) actually depends on h = |h|. Finally, with ergodicity assumption for the materials of interest, the averages over Ω converge to the averages over S. The spatial averages over Ω and Ω (r) are denoted ⟨.⟩ and ⟨.⟩ (r) respectively. The volume fraction of phase r in the RVE is thus defined by :

c r = p (r) (x) = χ (r) (x) (2.5)
and the second-order correlation function becomes :

p (rs) (h) = χ (r) (x)χ (r) (x + h) (2.6)
The covariance functions are finally given by :

C rs (h) = p (rs) (h) -c r c s (2.7)

Localization

The RVE being well defined, the localization step then aims at establishing the linear relation between the local stress (or strain) fiels, and their macroscopic counterpart :

ε(x) = A(x) : E and σ(x) = B(x) : Σ, (2.8)
where A(x) is the localization tensor and B(x) the concentration tensor. As the microstructure can not be entirely described, the local fields can not be characterized pointwisely either. The averages of the local fields are introduced and are defined by :

⟨ε(x)⟩ = 1 Ω Ω ε(x) dx = E and ⟨σ(x)⟩ = 1 Ω Ω σ(x) dx = Σ (2.9)
The micromechanical approach consists in considering the averages over each volume phase :

ε (r) = ⟨ε⟩ (r) = ⟨A(x)⟩ (r) : E = A (r) : E σ (r) = ⟨σ⟩ (r) = ⟨B(x)⟩ (r) : Σ = B (r) : Σ (2.10)
These relations imply :

⟨A⟩ = I and ⟨B⟩ = I (2.11)
where I is the 4 th -order identity tensor. The determination of A (r) and B (r) is addressed by means of homogenization schemes, depending on the microstructural morpohology.

Homogenization

The homogenized properties are obtained by the combination of the localization equations (2.8), the local constitutive laws and the averages of the local fields. Expliciting the local constitutive laws into (2.9), the effective stiffness tensor L is expressed by : L = ⟨L : A⟩ (2.12)

and the effective compliance tensor M :

M = ⟨M : B⟩ (2.13)

Eshelby problem

The Eshelby inclusion problem has been the base for the development of a variety of homogenization schemes. Eshelby considered the problem of an infinite medium of which a small region is submitted to a stress-free strain (the homogeneous inclusion problem). The inclusion problem can be extended to a heterogeneous inclusion problem, that is a material inhomogeneity [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF].

The homogeneous inclusion problem

Let us consider a linear elastic infinite homogeneous media Ω with a stiffness tensor L, initially in an undeformed configuration. A subdomain I ∈ Ω is submitted to a stress-free strain. The eigenstrain ε t displayed by the material in I is constrained by the surrounding medium, resulting in the real compatible deformation ε I . This accommodation also causes a perturbation of Ω \ I in the neighborhood of the ellipsoidal inhomogeneity. Considering small perturbations, the deformation of the surrounding medium is considered to be null far from I. By linearity, the free deformation and the real deformation in I are related through the 4 th -order Eshelby tensor S E :

ε I = S E : ε t (2.14)
The elastic local constitutive law in Ω \ I reads σ(x) = L : ε(x), while the expression in the inclusion is given by σ(x) = L : (ε(x) -ε t ), that can be rewritten :

σ = L : ε + τ (2.15)
where τ = -L : ε I is the polarization symmetric 2 nd -order tensor. The equations of the problem yield :

             σ(x) = L(x) : ε(x), ∀x ∈ Ω \ I and σ(x) = L : ε(x) + τ , ∀x ∈ I div σ(x) = 0, curl( t curl ε(x)) = 0, ∀x ∈ I ∪ Ω \ I ε = 1 2 ∇u + t ∇u , ∀x ∈ I ∪ Ω \ I (2.16)
The deformation in I can be related to the polarization field, making use of the 4 th -order Hill tensor rather than the Eshleby tensor :

ε I = -P : τ , P = S E : L -1 (2.17)
where P is symmetrical and positive-definite [75] [76]. For an ellipsoidal inclusion, the Hill tensor can be computed. Homogeneous to an compliance tensor, this tensor depends on the elasticity of the infinite medium, the shape of the ellipsoid inclusion and its orientation. Its expression reads :

P = 1 4π|υ| ||x||=1 Γ(x)||υ -1 .x|| -3 dS x (2.18)
where

Γ(x) = [x ⊗ κ -1 ⊗ x]
is the Green operator and κ = x.L.x is the acoustic Christoffel tensor. The 2 nd -order tensor υ characterizes the geometry of the ellipsoid I, where the origin is taken in the center of the inclusion.

We consider the problem (2.16), where the infinite medium is submitted to a homogeneous loading u = ε.x far from the inclusion. The solution of this infinite loading problem is straightforwardly found by the following linear superposition :

ε I = -P : τ + ε (2.19)

The heterogeneous inclusion problem

Let us consider the same linear elastic infinite medium Ω of stiffness tensor L embedding an ellipsoidal region I, initially in an undeformed configuration. However, the inclusion is now made of a different material whose stiffness is L I . The inclusion is submitted to the polarization tensor τ , which leads to the local problem, in which the constitutive law in the inclusion is governed by L I :

             σ(x) = L(x) : ε(x), ∀x ∈ Ω \ I and σ(x) = L I : ε(x) + τ , ∀x ∈ I div σ(x) = 0, curl( t curl ε(x)) = 0, ∀x ∈ I ∪ Ω \ I ε = 1 2 ∇u + t ∇u , ∀x ∈ I ∪ Ω \ I (2.20)
The problem actually reduces to (2.16) by reexpressing the constitutive law in the inclusion : σ(x) = L : ε + τ ′ , where τ ′ = τ + (L I -L) : ε. According to the first Eshelby problem, the deformation in the inclusion reads :

ε I = -P : τ + (L I -L) : ε I (2.21)
The Hill tensor is invertible, so this equation can be rewritten as :

P -1 -L + L I : ε I = -τ (2.22)
where the interaction tensor introduced by Hill [77] appears :

L * = P -1 -L = L : (S E ) -1 -I (2.23)
The interaction tensor is positive-definite and homogeneous to an elasticity tensor. This quantity is related to the stress applied by the infinite medium on the inclusion, in reaction to the homogeneous strain of the inclusion. L * is also invertible, so that the deformation in the inclusion can finally be expressed :

ε I = -L * + L I -1 : τ (2.24)
The Eshelby problem thus consists in considering the material heterogeneities as a polarization field. This idea led to the development of homogenization schemes for various microstructures. First, the introduction of very few other ellipsoidal inclusions in the Eshelby problem led to the dilute scheme, where the interactions between the material heterogeneities are weak. Mori and Tanaka [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] then considered several inclusions embedded in a matrix.

Benveniste [START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF] extended the latter work for the homogenization of polycrystalline materials. The self-consistent scheme introduced by Christensen [START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF] consists in considering the homogenized properties as the reference medium for each phase, and is thus well-suited for polycristalline microstructures. Luciano and Barbero [START_REF] Luciano | Formulas for the stiffness of composites with periodic microstructure[END_REF] obtained closed-form expressions for the effective properties of periodic matrix-fiber microstructures. An alternative of homogenization schemes is the use of variational principles which allow to establish bounds on the effective properties.

Variational principles and inequalities

Indeed, it can be useful to obtain bounds for the effective properties of a composite material in terms of an energy minimization problem. Let us consider the elastic medium Ω submitted to the displacement field u = ε.x on the contour ∂Ω. We denote u * the trial displacement field such that u * ∈ U(ε), where U(ε) refers to the set of kinematically admissible displacement fields leading to the uniform strain field ε on the contour ∂Ω. Based on the minimization of the quadratic error ⟨(ε -ε * (u * )) : L : (εε * (u * ))⟩ ≥ 0, and where σ * is divergence-free, the following inequalities can be obtained :

   ⟨(ε * (u * ) -M : σ * ) : L : (ε * (u * ) : L : ε * (u * ))⟩ ≥ 0 2⟨σ * ⟩ : ε -⟨σ * : M : σ * ⟩ ≤ ⟨ε * (u * ) : L : ε * (u * )⟩ (2.25)
The bounds reduce to equalities when the trial fields coincide with the exact fields. The variational representation is derived by considering the quadratic energies of the medium. Let us denote the convex strain and stress elastic energy densities w (ε * (u * )) and v(σ * ) respectively :

     w (ε * (u * )) = 1 2 ε * (u * ) : L : ε * (u * ) v(σ * ) = 2⟨σ * ⟩ : ε -⟨σ * : M : σ * ⟩ (2.26)
The effective strain energy can be expressed by the following minimization problem :

w(ε) = 1 2 ε : L : ε = min u * ∈U (ε) ⟨w (ε * (u * ))⟩ (2.27)
and the dual complementary energy principle gives :

ṽ(σ) = 1 2 σ : M : σ = min σ * ,div(σ * )=0 ⟨v (ε * (u * ))⟩ (2.28)
For any admissible choice of trial stress and strain fields, the following bounds hold :

max σ * ,div(σ * )=0 ⟨v (σ * )⟩ ≤ 1 2 ε : L : ε ≤ min u * ∈U (ε) ⟨w (ε * (u * ))⟩ (2.29)

Voigt and Reuss bounds

The general bounds (2.29) hold for any kinematically admissible displacement trial field ε * and statically admissible trial field σ * . The most basic assumption for the upper bound is to consider the homogeneous strain field equal to the macroscopic strain, that is ε * (x) = ε. The upper bound is then given by :

ε : L : ε ≤ ⟨ε : L : ε⟩ (2.30)
The effective stiffness tensor thus reduces to the Voigt estimate L V [START_REF] Voigt | Lehrbuch der kristallphysik[END_REF] :

L V = ⟨L(x)⟩ (2.31)
Dually, the Reuss bound consists in the simplest assumption for the stress field, that is the homogeneous field σ * (x) = σ = ⟨M⟩ -1 : ε. The derived lower bound in (2.29) can be expressed :

σ : L-1 : σ ≤ σ : L -1 : σ (2.32)
and leads to :

L R = ⟨L -1 (x)⟩ -1 (2.33)
The general bounds (2.29) can thus reduce to the Voigt-Reuss ones :

ε : ⟨L -1 ⟩ -1 : ε ≤ ε : L : ε ≤ ε : ⟨L⟩ : ε (2.34)
In the above bounds, only the elastic properties and the volume fractions are required, thus providing wide bounds. The estimations can be improved by taking into account the heterogeneity of the fields. They can be improved by means of finer variational formulation.

Hashin and Shtrikman bounds

Hashin and Shtrikman [START_REF] Hashin | On some Variational Principles in Anisotropic and Nonhomogeneous Elasticity[END_REF] introduced a variational principle that admits heterogeneous stress and strain fields. Their formulation is based on the substitution of the real problem with a homogeneous reference medium of stiffness tensor L 0 subjected to the same loading as the initial problem. This reference thermoelastic problem involves the polarization tensor τ * defined through :

σ * = L 0 : ε * + τ * (2.35)
The polarization field is assumed to be of the form τ * = (L -L 0 ) : η * , where η * is a 2 nd -order tensor with the same dimension unit as the deformation. However, it is not required to respect the compatibility conditions. Based on this polarization formulation, the energy and complementary energy principles provide the following inequalities :

-∆Ψ(η * ) + H 0 (η * ) ≤ 1 2 ε : L : ε ≤ H 0 (η * ) + ∆Φ(η * ) (2.36)
where H 0 is the Hashin-Shtrikman functional that reads :

2H 0 (η * ) = ε : L 0 : ε + η * : (L -L 0 ) : [ε + ε * (η * ) -η * ] (2.37)
and ∆Ψ and ∆Φ are quadratic quantities defined by :

2∆Ψ(η * ) = [ε * (η * ) -η * ] : L 0 : (M -M 0 ) : L 0 : [ε * (η * ) -η * ] (2.38)
and

2∆Φ(η) = [ε * (η * ) -η * ] : (L -L 0 ) : [ε * (η * ) -η * ] (2.39)
The formulation (2.36) allows to consider the Hashin-Shtrikman functional (2.37) as the upper or lower bound of the effective elastic energy, depending on the choice of L 0 . It will indeed induce the sign of the quadratic quantities ∆Ψ and ∆Φ, providing an indication of the choice of the reference medium. More precisely :

H 0 ≥ 1 2 ε : L : ε, L 0 ≥ L(x) ∀x ∈ Ω (2.40)
and

H 0 ≤ 1 2 ε : L : ε, L 0 ≤ L(x) ∀x ∈ Ω (2.41)
∆Ψ and ∆Φ are actually exactly null if the trial polarization field is taken as the exact polarization field τ * = τ = (L -L 0 ) : ε, or equivalently if η * is taken equal to the exact strain field ε. Besides, the Hashin-Shtrikman functional reduces to the potential energy.

The choice of the stiffness of the reference medium thus provides the variational principles :

w(ε) = min ⟨η * ⟩ H 0 (η * ), L 0 ≥ L(x) ∀x ∈ Ω (2.42) and w(ε) = max ⟨η * ⟩ H 0 (η * ), L 0 ≤ L(x) ∀x ∈ Ω (2.43)
or more generally, for any choice of L 0 :

w(ε) = stat ⟨η * ⟩ H 0 (η * ), ∀L 0 (2.44)
In practice, the homogeneous reference medium L 0 can be chosen as the lowest and the highest stiffness tensors in a quadratic sense of the local phases. Hashin and Shtrikman adressed the problem of a composite made of two isotropic phases and obtained bounds for the bulk and shear moduli, making the assumption that the polarization field was uniform within each phase [START_REF] Hashin | On some Variational Principles in Anisotropic and Nonhomogeneous Elasticity[END_REF]. They pointed out that the lower bound of the bulk modulus was the exact result for a composite sphere assemblage.

Periodic elastic cell problem

Let us consider a periodic cell. The heterogeneous problem reads :

                   σ(x) = L(x) : ε(x) ∀x ∈ Ω div σ (x) = 0, curl( t curl ε(x)) = 0 ∀x ∈ Ω ε(x) = 1 2 ∇u(x) + t ∇u(x) ∀x ∈ Ω ε(x) = E , ε(x) # , σ(x).n -# (2.45)
The symbol # denotes that ε(x) is a periodic field, and -# that σ(x).n is anti-periodic.

Introducing the fictitious linear elastic homogeneous media L 0 , the constitutive law can be rewritten as σ(x) = L 0 : ε(x) + τ (x), and where the polarization field is defined by τ (x) = L(x) -L 0 : ε(x). The local strain is splitted into its average value and a fluctuation term :

ε(x) = E + ε * (x), (2.46)
or equivalently :

u(x) = E.x + u * (x) (2.47)
By substituting τ (x) with its expression, the equilibrium condition can be rewritten :

div σ(x) = div L 0 : ε(x) + div τ (x) = 0 (2.48)
The solution of this linear problem is given by the periodic Lippman-Scwhinger equation in direct and in Fourier space, where the hat symbol denotes the Fourier transform (alternatively denoted FT ) :

     ε(x) = E -Γ 0 * τ (x) ∀x ∈ Ω ε(ξ) = - Γ0 (ξ) : τ (ξ) ∀ξ ̸ = 0, ε(0) = E (2.49)
where Γ 0 is the Green operator in bold symbol, not to be confused with the classical gamma function Γ. Its analytical expression for an isotropic elastic constituent in Fourier space is given by :

Γ0 khij = - 1 4 (ξ h N ki ξ j + ξ h N kj ξ i + ξ k N hi ξ j + ξ k N hj ξ i ) (2.50) with N 0 ki = κ 0 ik -1 = L 0 ijkl ξ l ξ j -1
, κ 0 ik being the acoustic tensor.

Contribution to numerical homogenization : smoothing filter for FFT-based methods

The basic-scheme [START_REF] Moulinec | A Fast Numerical Method for Computing the Linear and Nonlinear Mechanical Properties of Composites[END_REF] consists in iteratively solving the Lippman-Scwhinger equation with a fixed-point algorithm. It can be summarized by the following equation, where the field at iterate i + 1 is the updated value from iteration i :

εi+1 (ξ) = - Γ0 (ξ) : FT L -L 0 : ε i (ξ), ∀ξ ̸ = 0, εi+1 (0) = E (2.51)
When the material is isotropic, the condition of convergence can be reduced to the following choice of the homogeneous bulk and shear moduli :

κ 0 = 1 2 (κ min + κ max ) , µ 0 = 1 2 (µ min + µ max ) (2.52)
where the indexes . min and . max refer respectively to the minimum and maximum values of the moduli among the phases. Despite the reliability of the basic scheme for linear and non-linear problems, its convergence can be slow.

Accelerated algorithm

Several schemes were then proposed to improve the convergence speed of the initial basic scheme developped by [START_REF] Moulinec | A Fast Numerical Method for Computing the Linear and Nonlinear Mechanical Properties of Composites[END_REF], e.g. [START_REF] Müller | Mathematical vs. Experimental Stress Analysis of Inhomogeneities in Solids[END_REF], [START_REF] Vinogradov | An accelerated FFT algorithm for thermoelastic and non-linear composites[END_REF], [START_REF] Zeman | Accelerating a FFTbased solver for numerical homogenization of periodic media by conjugate gradients[END_REF]. Moulinec and Silva [START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF] compared three of them : the Eyre-Milton [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF], the augmented-Lagrangian [START_REF] Michel | A Computational Method Based on Augmented Lagrangians and Fast Fourier Transforms for Composites with High Contrast[END_REF], and the polarization schemes [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF]. They showed that the two formers are specific cases of the latter one. They rewrote the algorithm such that the similarities are highlighted. The polarization algorithm [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF] corresponds to the augmented-Lagrangian method if a = b = 1, and reduces to the Eyre-Milton if a = b = 2. Contrary to the basic scheme, the considered fileds e i and s i , respectively refering to the strain and stress field, are not required to ensure the compatibility and the equilibrium conditions respectively, before convergence. The algorithm is presented in 1, where the Fast Fourier Transform operator is denoted by FF T . [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF] rewritten following Moulinec and Silva [START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF] Require: e i , s i , a, b, L 0 , ϵ the tolerance while div s i > ϵ and curl( t curl e i ) > ϵ do

Algorithm 1 Polarization algorithm

s i = s i + (1 -b)L 0 : e i (x).s i (x) e i (x) = a s i (x) -b L 0 : e i (x) êi (ξ) = FF T e i (x) êi (ξ) = - Γ0 (ξ) : êi ∀ξ ̸ = 0, êi (0) = b E e i (x) = FF T -1 êi (ξ) e i+1 (x) = (L(x) + L 0 ) -1 : (s i (x) + L 0 : e i (x)) s i+1 (x) = L(x) : e i+1 (x

) end while

As shown by [START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF], the best convergence rate is provided for a = b = 2, i.e. for the Eyre-Milton scheme. For isotropic constituents that are neither voids nor infinitely rigid, the elastic moduli should be chosen as :

κ 0 = √ κ min κ max , µ 0 = √ µ min µ max (2.53)
Therefore, in the following the FFT-based method is computed using the Eyre-Milton scheme, i.e. a = b = 2.

Smoothing splines for FFT-based solvers : data filtering

We present the automatic periodic smoothing filter developed by Morin et. al. [START_REF] Morin | Periodic smoothing splines for FFT-based solvers[END_REF], that smoothes the discontinuous fields responsible for the oscillations. Based on polynomial smoothing, this approach has the advantage not to modify the FFT-based algorithm and only constitutes an additional step before (reducing the discontinuities on the material properties) or after (removing the oscillations on the resulting fields). Following [START_REF] Morin | Periodic smoothing splines for FFT-based solvers[END_REF], it is preferred to apply it on the microstructure as an input before the computation, adding a smooth interface between the phases that can actually be present on real composite materials. Besides, this filter is both easy to implement and numerically inexpensive.

After presenting the method, the filter will be investigated on elastic 2D microstructures for which the analytical solutions are known, allowing to compare the FFT-based computation to reference results. The oscillations issues appearing regardless of the constitutive law, this choice is not reductive. We first present the periodic smoothing filter, as well as the example given in the original article. A new example on a challenging microstructure is then considered, which allows to confront the limits of the filter.

Data filtering

Let us consider the periodic three-dimensional function f defined on the interval [0,

T 1 ] × [0, T 2 ] × [0, T 3 ] discretized with N 1 × N 2 × N 3 voxels.
The volume is uniformly gridded with spatial scales : ∆x = T 1 /N 1, ∆y = T 2 /N 2 and ∆z = T 3 /N 3 . We denote f jkl = f (x j , y k , z l ), with x j = j∆x, y k = k∆y, z l = l∆z and j, k, l = 0, ..., N -1. The smooth estimate f of f is given by the minimization of the following functional [START_REF] Craven | Smoothing noisy data with spline functions[END_REF] :

1 N 1 N 2 N 3 N 1 -1 j=0 N 2 -1 k=0 N 3 -1 l=0 f jkl -fjkl 2 + s T 1 T 2 T 3 T 1 0 T 2 0 T 3 0 f (m) (u, v, w) du dv dw (2.54)
and where the superscipt (m) corresponds to the m th derivative. The solution of Eq. (2.54) is a polynomial smoothing spline of degree 2m -1. The smooth parameter s, which must be chosen, controls the fidelity to the data, measured by the first term, and the roughness of the solution, given by the second term. The polynomial solution is given by the following expression, where • denotes a pointwise product :

f = FT -1 (Υ 3 (s) • F T (f ))
(2.55)

Υ jkl = 1 1 + s ξ 2m j + υ 2m k + ζ 2m l (2.56)
with ξ, υ, ζ the coordinates in the Fourier domain. The filter preserves the mean value of the real function f . In the following, the degree of the polynomial is restricted to m = 1, as a higher order could lead to unphysical solutions such as negative material properties [86] [15].

Checkerboard-like microstructure

We consider the incompressible periodic matrix-fibers microstructure that was considered by Morin et. al. in [START_REF] Morin | Periodic smoothing splines for FFT-based solvers[END_REF]. The fibers are oriented along the z direction and have a squared section as illustrated in Fig. 2.1, and each phase is isotropic, such that the material displays a transversely isotropic symmetry. The periodic unit cell is submitted to the out-of-plane shear loading strain ε = ε 0 (e 1 ⊗ e 3 + e 3 ⊗ e 1 ). The problem can thus be reduced to a two-dimensional problem in the symetry plane. The exact solution of this microstructure is given by Obnosov [START_REF] Yurii | Periodic Heterogeneous Structures: New Explicit Solutions and Effective Characteristics of Refraction of an Imposed Field[END_REF] (see Appendix C) for a conductivity problem. Its extension to elasticity is given in Morin et al. [START_REF] Morin | Periodic smoothing splines for FFT-based solvers[END_REF]. The dimensions of the unit cell are

L 1 = L 2 = 4 with x ∈ [-L 1 /2, L 1 /2] and y ∈ [-L 2 /2, L 2 /2]. The squared inclusion is of dimension L 1 /2 × L 2 /
2 and is denoted as phase (2) while phase (1) denotes the matrix. The elastic tensor of phase i, with i = 1, 2, reads :

L (i) = 2µ (i) 1 0 0 1 (2.57)
And we define the contrast c = µ (2) µ (1) . The checkerboard microstructure is illustrated in Fig. 2.5, where Ω 1 refers to the grain of which the local coordinate system corresponds to the macroscopic one, and Ω 2 refers to the grain which is rotated with the angle 90°. The dotted line at y = L 2 /2 specifies where the local fields will be plotted in the following. In the macroscopic system of coordinates, the stiffness tensors of orientations ( 1) and ( 2) read :

L (1) = 2 µ (a) 0 0 µ (b) , L (2) = 2 µ (b) 0 0 µ (a)
(2.59)

And the effective stifness tensor reads [START_REF] Craster | Four-Phase Checkerboard Composites[END_REF] :

L = 2μ exact 1 0 0 1 , μexact = µ (a) µ (b) (2.60)
Figure 2.5 : Checkerboard microstructure

The analytical solution of this problem is given by [START_REF] Craster | Four-Phase Checkerboard Composites[END_REF] (see Appendix C). We investigate the effect of the filter (2.55) for several material contrasts c = µ (b) /µ (a) with values 2, 10 and 100. The moduli µ 2323 (x) of the exact and smoothed microstructure are represented in Figs 2.6a, 2.6b and 2.6c.

Effective properties

Setting the smooth parameter with the value s = 5 × 10 -3 , and considering the constrasts c = 2, 10 and 100, we compare the performance of the filter on the effective stiffness tensor with the unfiltered simulation and the exact result. On Tables 2.1, 2.2 and 2.3, the effective shear modulus is given, the second column gives the relative error ϵ with respect to the exact modulus such that ϵ = |μ exactμs |/μ exact . The last column give the values of the off diagonal moduli of the effective stiffness tensor. For a small contrast c = 2, the smoothing filter does not improve the effective properties and even deteriorates the accuracy of the unfiltered effective stiffness, as reported in 2.1. This effect is accentuated for a higher but moderate contrast c = 10, see in Table 2.2. For a high contrast c = 100, the filtered computation is entirely unable to give accurate results. Although the effective shear modulus obtained by the unfiltered simulation s = 0 is acceptable, the off diagonal terms are not null, showing the limits of the FFT computation on this challenging microstructure.

Local fields

The visualization of the fields in Figs. 2.8 and 2.9 shows that the smoothing filter on the mechanical properties allows to circumvent the spurious oscillations on the results. Figs. 2.7 show the component ε 13 along the x-direction at y = L 2 /2 (dashed line in Fig. 2.5) and for the contrast c = 100. This graph indicates that the rough singularities in the corners of the grains are expected to display abrupt jumps according to the exact fields. These are not recovered by the FFT simulations, even though the mean values of the fields may lead to the correct effective shear modulus depending on the contrast. 2 Homogenization in elasticity and numerical contribution 61

Proposed weighted smoothing filter

Weighted smoothing filter

The filter applied to the rigid tensor overestimates the effective properties, especially for high contrasts. We thus investigate a combination of the filter on both rigidity and compliance properties Ľ and ( M) -1 , as a way to soften the effective estimates. We introduce the parameter h, referred to as the weighting parameter in the following. The combined filter reads :

Ľ = h Ľ + (1 -h)( M) -1 , 0 ≤ h ≤ 1 (2.61)
As seen in the previous section for the checkerboard microstructure, the effective shear modulus is overestimated, and this effect is emphasized as the contrast increases. The introduction of the smoothed compliance tensor will calculate lower stiffness properties and lead to lower effective shear moduli. The effect of the weighted smoother is firstly investigated on the checkerboard-like microstructure, then on the checkerboard polycrystal, in order to seek for an "optimal" couple (s, h).

Smoothed checkerboard-like microstructure

The effect of the weighting parameter h on the material properties results in lowering the value of the local shear modulus on the interfaces, as seen in Fig. 2.10. As well as in the previous section, the stress component σ 13 (x) is shown in the top right corner of the squared inclusion in Fig. 2.11. The weighted smoother with parameters s = 0.001 and h = 0.8 is closer to the analytical result in comparison with the other configuration. We continue the investigation of the filter on the checkerboard microstructure. We fix the smooth parameter to the value s = 5 × 10 -3 . We seek for the best weighting parameter h that ensures the following conditions :

Effective properties for the weighted filter

Off diagonal terms

The first remarkable point about the effective stiffness tensor obtained from the double filtered microstructures is that the off diagonal terms are systematically null (numerical zeros, around 10 -10 Pa), regardless of the value of h.

Effective shear modulus

Let us now look into the value of the effective shear modulus. Figs. 2.14a, 2.14b and 2.14c show that the introduction of the smoothing compliance leads to a lower estimate of the effective shear modulus compared to h = 1, as expected by its effect on the microstructure. For the contrasts c = 2, 10 and 100, the accuracy of the shear modulus is obtained for the weighting parameter h = 0.4, 0.

Local fields for the weighted filter

The component ε 13 (x) is plotted for c = 2 in Figs. 2.15. We observe a softening of the local response, especially on the center corner. Despite the low contrast c, the FFT computation does not provide accurate local results. Besides, there is almost no difference between the multiple parameters h.

For the large contrast c = 100, the results are shown in Figs. 2.16. The close-up shows that the high singularity at x = 0 results in a collapse for the FFT fields, whereas the analytical solution expects a jump. This effect is due to the fact that the coordinate (x = 0, y = L 2 /2) represents the singular point of the microstructure that joins the 4 squared single crystals. By smoothing the microstructure, the shear modulus µ(x) have a wells shape in this region, producing the collapse effect on the resulting fields. The filter indeed removes the oscillations, however no significant difference is observed between the different values of h.

Conclusion

In this chapter, after presenting the fundamentals of homogenization in linear elasticity, we introduced the FFT-based solver for iteratively solving the periodic Lippman-Schwinger equation in elasticity. We also presented the robust method developed by [START_REF] Morin | Periodic smoothing splines for FFT-based solvers[END_REF] for smoothing the discontinuous interface between the constituents, obviating the resulting spurious oscillations of the FFT-based calculations. Based on polynomial smoothing splines, it consists in minimizing a functional that balances the fidelity to the data and the smoothness of the approximate function, where the smooth parameter s adjusts the width of the interface. Contrary to the other smoothing strategies such as [START_REF] Willot | Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF] and [START_REF] Gélébart | Filtering material properties to improve FFTbased methods for numerical homogenization[END_REF], the filtering data is automatic and does not modify the FFT-based algorithm. We applied the filter on a challenging microstructure with periodic singular points that joins 4 squared section crystals. For high contrasts, the effective properties are not accurate. Indeed, the effective shear modulus is overestimated and the off-diagonal terms are not null as should be expected analytically. We thus proposed a composed filter that weights the contribution of the smooth stiffness tensor and the inverse of the smooth compliance tensor. This combination results in the lowering of the stiffness modulus on the smooth interfaces. The goal was to seek for an "optimal" weighting parameter 0 ≤ h ≤ 1 for s fixed, that could improve the accuracy both on the effective properties and on the local fields. The effective shear modulus is indeed lowered and h can be set such that it is accurate. However, it has to be chosen under 0.5, significantly modifying the real microstructure. Besides, there is no theoretical discrimination criteria for this choice. However, for any h < 1, the off diagonal terms are systematically null. Locally, the fields remain similar regardless of the value of the weighted filter. Nonetheless, this microstructure is particularly challenging, and the effect on the squared inclusions embedded in a matrix showed an actual improvement, although h cannot be rigorously chosen.

Introduction

Unlike linear elasticity, the effective properties of heterogeneous media composed of viscoelastic constituents display mechanical features that are not present at the local scale. Indeed, the mixing of viscoelastic constituents with short memory (discrete spectrum with a finite number of relaxation times) gives rise to long-memory effects at the macroscopic scale (infinite relaxation times) [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF]. This can be evidenced by means of the correspondence principle, which consists in considering the Laplace-Carson transform of the viscoelastic equations. In the spectral domain, these equations take the form of a symbolic elastic problem that can be solved with the classical homogenization methods [START_REF] Hashin | Viscoelastic Behavior of Heterogeneous Media[END_REF]. The viscoelastic effective moduli are obtained by the inverse Laplace transform; however the inversion can generally not be evaluated in closed-form but only numerically. This is usually approximated by means of the collocation method, allowing to estimate the effective properties with a finite number of characteristic times [START_REF] Schapery | Methods of transform inversion for viscoelastic stress analysis[END_REF] [90] [48] [START_REF] Ricaud | Effective properties of linear viscoelastic heterogeneous media: Internal variables formulation and extension to ageing behaviours[END_REF]. The extension to fractional visoelasticity with this method is possible, see for example [START_REF] Dinzart | Self-consistent approach of the constitutive law of a two-phase visco-elastic material described by fractional derivative models[END_REF] for polycrystalline materials, and [START_REF] Gallican | Homogenization estimates for polymer-based viscoelastic composite materials[END_REF] for a viscoelastic matrix reinforced with spherical elastic particles. However, the homogenization methods relying on the correspondence principle can not provide the second-moments of the fields. The inaccessibility of the intraphase fluctuations is due to the impossibility of evaluating their inverse Laplace transforms.

Alternatively, for a possible characterization of the fields statistics, one may use incremental methods applicable in the direct domain. These methods rely on variational approaches. Specifically, we take an interest on the Effective Internal Variable (EIV) method developed by Lahellec and Suquet in 2007 [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF]. The EIV method consists in solving iteratively the heterogeneous viscoelastic problem in the direct domain by considering a thermoelastic problem at each time step. This approach is built upon the thermodynamic formulation of materials with dissipative and conservative effects, referred to as Generalized Standard Materials (GSM), initially introduced by Halphen and Nguyen in 1975 [START_REF] Halphen | Sur les Matériaux Standards Généralisés[END_REF]. At each time step of the considered loading, these two potentials are derived into a unique incremental potential that thus takes into account both conservative and dissipative effects. This formulation allows to make use of the variational procedures developed by Ponte Castañeda in 1991 [3] and 2002 [4]. This method provides accurate results, see for example [START_REF] Boudet | An incremental variational formulation for the prediction of the effective work-hardening behavior and field statistics of elasto-(visco)plastic composites[END_REF] where the authors consider work-hardening behaviour on elasto-viscoplastic composites. The incremental method was also extended into an alternative strain-rate formulation [START_REF] Lahellec | Effective response and field statistics in elastoplastic and elasto-viscoplastic composites under radial and non-radial loadings[END_REF] which has been applied to linear viscoelastic polycrystals [START_REF] Badulescu | Field statistics in linear viscoelastic composites and polycrystals[END_REF]. However, the EIV method presents some inherent limitations for specific cases. The mathematical structure of the EIV equations was explicited in 2020 by Idiart, Lahellec and Suquet in companion papers [START_REF] Martín | Model reduction by meanfield homogenization in viscoelastic composites. I. Primal theory[END_REF], [START_REF] Martín | Model reduction by meanfield homogenization in viscoelastic composites. II. Application to rigidly reinforced solids[END_REF] and [START_REF] Lahellec | Model reduction by meanfield homogenization in viscoelastic composites. III. Dual theory[END_REF], where they showed that non-radial and harmonic loadings could contrevene fundamental hypotheses of the method.

For the purpose of incorporating fractional viscoelastic constituents into the EIV method, the following question arises : how can the the thermodynamic potentials be derived ? The key idea relies on the rheological interpretation of the fractional dashpot as a superposition of classical Maxwell elements (see Chapter 1). The mechanical state of the material can thus be described by means of multiple viscous strains, or equivalently internal variables. As a result, there exist thermodynamic potentials in the framework of the GSM. They are directly obtained from the linear superposition of those of the individual branches. Lion stressed out that this formulation was naturally consistent with the second principle of thermodynamics [START_REF] Lion | On the thermodynamics of fractional damping elements[END_REF]. Hence, the extension of the EIV method for composites with fractional viscoelastic constituents implies the consideration of multiple internal variables for one constituent. Precisely, Tressou addressed the extension of the EIV method to generalized Maxwell constituents, i.e. for several internal variables [START_REF] Tressou | Contribution à l'homogénéisation des milieux viscoélastiques et introduction du couplage avec la température par extensions d'une approche incrémentale directe[END_REF]. In [START_REF] Tressou | Application of the variational EIV approach to linear viscoelastic phases governed by several internal variables -Examples with the generalized Maxwell law[END_REF], Tressou et al. obtained accurate estimates considering rheological assemblies with two and four internal variables, and for different contrasts of characteristic times. In that respect, a discrete rheological description of fractional constituents is required. This can be obtained by discretizing the spectrum of which the expression is explicitely known. More precisely, the four fractional assemblies of viscoelasticity [START_REF] Caputo | Linear models of dissipation in anelastic solids[END_REF] can involve the Mittag-Leffler function. The spectrum thus takes the form of a rational fraction involving polynomials of non-integer order. In this study, we make use of the discretization procedure proposed by Papoulia et al. [START_REF] Papoulia | Rheological representation of fractional order viscoelastic material models[END_REF].

The outline of the present chapter is as followed. For introductory purpose, we first present the commonly used correspondence principle for the application of homogenization methods for viscoelasticity. Then, we introduce the EIV method for multiple internal variables, as developed by Tressou et al. [START_REF] Tressou | Application of the variational EIV approach to linear viscoelastic phases governed by several internal variables -Examples with the generalized Maxwell law[END_REF] for Generalized Maxwell constituents. After that, the discretization of integral viscoelastic laws will be dealt with. We will primarily introduce in detail the method of Papoulia et al. [START_REF] Papoulia | Rheological representation of fractional order viscoelastic material models[END_REF] and the application of their method. We also compare their approach with the classical collocation method. Finally, we present a new application of the discretization procedure on the integral Mittag-Leffler function. This result will allow to extend the EIV formulation for considering any fractional relaxation function in the EIV method.

Homogenization in viscoelasticity

This section deals with homogenization methods for material composites with linear viscoelastic constituents. Generally, these problems are tackled by means of the correspondence principle [START_REF] Mandel | Cours de Mécanique des Milieux Continus[END_REF] involving the transformation of the local problem into the Laplace-Carson domain. This method allowed [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF] to highlight that the local coupling of elastic and viscous effects and the spatial interactions within the different phases gave rise to long-memory effects at the macroscopic scale. Alternatively, the use of a time-integration approach allows to circumvent the computational cost of integral approaches [START_REF] Ricaud | Effective properties of linear viscoelastic heterogeneous media: Internal variables formulation and extension to ageing behaviours[END_REF].

Correspondence principle

Let us consider the heterogeneous domain Ω composed of N phases characterized by the local viscoelastic stiffness tensors L (r) (t) with r = 1, ..., N . The medium is submitted to a uniform strain loading ε(t) = E(t) from t = 0 to t = T . The local problem reads :

                         σ(x, t) = t 0 L(x, t -ξ) : ε(x, ξ) dξ ∀(x, t) ∈ Ω × [0; T ] div σ(x, t) = 0, curl( t curl ε(x, t)) = 0 ∀(x, t) ∈ Ω × [0; T ] ε(x, t) = 1 2 ∇u(x, t) + t ∇u(x, t) ∀(x, t) ∈ Ω × [0; T ] ε(t) = E(t) ∀t ∈ [0; T ] (3.1)
The effective constitutive law is given by the volume average of the local law, yielding the effective viscoelastic tensor [START_REF] Hashin | Viscoelastic Behavior of Heterogeneous Media[END_REF] :

σ(t) = t 0 L(t -ξ) : dε dξ (ξ) dξ, where L(t) = t 0 L(x, t -ξ) : dA dξ (x, ξ) dξ (3.2)
Homogenization methods for elasticity can not be directly applied to solve this problem because of the integral form. One can use the correspondence principle [START_REF] Mandel | Cours de Mécanique des Milieux Continus[END_REF], which consists in transforming the viscoelastic equations (3.1) into the Laplace-Carson domain [START_REF] Hashin | Viscoelastic Behavior of Heterogeneous Media[END_REF], where they read : The use of Laplace-Carson transform leads to a symbolic elastic problem (3.3). For a fixed value of the variable p, classical homogenization methods can be applied. In the case of a harmonic strain loading E(t) = Êe iωt with i 2 = -1, the complex stress and strain fields σ(x, iω) and ε(x, iω) are related through the constitutive law :

                   σ(x, p) = L(x, p) : ε(x, p) ∀x ∈ Ω div σ(x, p) = 0, curl( t curl ε(x, p)) = 0 ∀x ∈ Ω ε(x, p) = 1 2 ∇û(x, p) + t ∇û(x, p) ∀x ∈ Ω ε(p) = Ê(p) , ( 3
σ(x, iω) = L(x, iω) : ε(x, iω), with L(x, iω) = L ′ (x, iω) + iL ′′ (x, iω) (3.5)
The local problem in the quasi-static regime reads :

                   σ(x, iω) = L(x, iω) : ε(x, iω) ∀(x, ω) ∈ Ω × [0; +∞] div σ(x, iω) = 0, curl( t curl ε(x, iω)) = 0 ∀(x, ω) ∈ Ω × [0; +∞] ε(x, iω) = 1 2 ∇û(x, iω) + t ∇û(x, iω) ∀(x, ω) ∈ Ω × [0; +∞] ε(iω) = E(iω) , (3.6)
3 EIV approach for fractional viscoelastic constituents

71

The effective viscoelastic tensor is obtained by the average stress over the whole domain :

σ(iω) = L(iω) : Ê(iω), where L(iω) = L′ (ω) + i L′′ (ω) (3.7)
Regarding the upscaling procedure, Suquet [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF] applied the correspondence principle to a composite material made of viscoelastic constituents.This author considered local Maxwellian constitutive laws, which display short-memory effect as they are defined with a unique characteristic time. He showed that the macroscopic behaviour displayed longmemory effects through an additional integral term in the expression of the effective viscoelastic tensor that is not present at the local scale. In his PhD thesis [START_REF] Gallican | Homogenization estimates for polymer-based viscoelastic composite materials[END_REF], Gallican explicited the same procedure to a mixture of Zener constituents, resulting in the same conclusion of macroscopic long-memory effects.

The inversion of the effective properties in the direct domain can not be obtained exactly, and in general requires a numerical approximation by means for example of the collocation method [START_REF] Schapery | Methods of transform inversion for viscoelastic stress analysis[END_REF] [START_REF] Rekik | Optimization of the collocation inversion method for the linear viscoelastic homogenization[END_REF]. Alternatively, the problem can be formulated by means of internal variables, as proposed by Ricaud and Masson [START_REF] Ricaud | Effective properties of linear viscoelastic heterogeneous media: Internal variables formulation and extension to ageing behaviours[END_REF]. They took advantage of the Prony-Dirichlet series expansion used for the approximation of the inverse effective moduli, by noting that this was equivalent to considering internal variables. Indeed, each term of the approximated relaxation function series is associated to a characteristic time τ k and an internal variable α τ k . The evolution law of the latters follow a Maxwell differential equation. The integral constitutive law thus reduces to :

σ(t) = L : ε(t) - K k=1 α τ k (t) (3.8)
with K the finite number of terms in the Prony series. The homogenization problem can be reduced to the resolution of an implicit system of linear equations that turns out to be exact for some specific two-phase microstructures. This formulation was applied to an ageing porous viscoelastic matrix reinforced with elastic inclusions [START_REF] Ricaud | Effective properties of linear viscoelastic heterogeneous media: Internal variables formulation and extension to ageing behaviours[END_REF], and to an ageing polycrystalline microstructure, making use of the self-consistent approach [START_REF] Masson | Incremental homogenization approach for ageing viscoelastic polycrystals[END_REF].

FFT-based homogenization of viscoelastic composites

The extension of FFT-based methodspresented in the Chapter 2 in the context of elasticityto viscoelastic materials is precisely done by means of the correspondence principle. This is done by considering the Laplace-Carson variable as purely imaginary : p = iω. Let us consider the harmonic strain loading with pulsation ω : ε(t) = εe iωt . The quasi-static local viscoelastic problem in the Laplace-Carson domain reads :

                   σ(x, iω) = L(x, iω) : ε(x, iω) ∀(x, ω) ∈ Ω × [0; +∞] div σ(x, iω) = 0, curl( t curl ε(x, iω)) = 0 ∀(x, ω) ∈ Ω × [0; +∞] ε(x, iω) = 1 2 ∇û(x, iω) + t ∇û(x, iω) ∀(x, ω) ∈ Ω × [0; +∞] ε(iω) = E(iω) , ε # , σ.n # (3.9)
The macroscopic response is given by the relation :

⟨ σ⟩(iω) = L(iω) : ⟨ε⟩(iω) (3.10)
where the real part of the effective complex viscoelastic tensor refers to the elastic response of the RVE, and the imaginary part to the dissipative effects :

L(iω) = L′ (ω) + i L′′ (ω) (3.11)
The symbolic elastic problem (3.9) can be solved in the same way as in linear elasticity, but with considering complex fields. The reference homogeneous tensor L 0 is required to be real. For isotropic constituents that are neither voids nor rigid, the reference moduli are chosen such as [START_REF] Figliuzzi | Modelling the microstructure and the viscoelastic behaviour of carbon 4 Bibliography 133 black filled rubber materials from 3D simulations[END_REF] :

κ 0 = κ ′ min κ ′ max , µ 0 = µ ′ min µ ′ max (3.12)
Regarding fractional viscoelasticity, the correspondence principle is applied without difficulty. For example, the FFT-based solver has been applied in [START_REF] Gallican | Homogenization estimates for the effective response of fractional viscoelastic particulate composites[END_REF] for a fractional Zener matrix reinforced with spherical elastic inclusions, and for a polycrystalline microstructure with fractional Maxwellian slip systems [START_REF] Gallican | Homogenization estimates for polymer-based viscoelastic composite materials[END_REF].

Variational principle for viscoelasticity

As seen in section 2.1.3, variational principles are used in elasticity to construct bounds on the effective properties of heterogeneous materials. The bounds generally describe convex domains containing the effective moduli depending on the volume fraction. Regarding viscoelastic materials, Cherkaev and Gibiansky [START_REF] Cherkaev | Variational principles for complex conductivity, viscoelasticity, and similar problems in media with complex moduli[END_REF] derived variational principles in the spectral domain, leading to a saddle-point min-max problem.

Let us consider the viscoelastic body Ω submitted to a uniform harmonic strain loading of pulsation ω. The harmonic local stress and strain fields σ s and ε s are given by :

ε s (x, t) = Re[ε(x)e -iwt ] = ε ′ (x)cos(wt) -iε ′′ (x)sin(wt) σ s (x, t) = Re[σ(x)e -iwt ] = σ ′ (x)cos(wt) -iσ ′′ (x)sin(wt) (3.13)
where ε and σ are related through the complex viscoelastic tensor : L(x, iω) = L ′ (x, ω) + iL ′′ (x, ω). The constitutive law can be expressed in terms of the following system of real equations :

σ ′ (x) σ ′′ (x) = -L ′′ (x) L ′ (x) L ′ (x) L ′′ (x) ε ′′ (x) ε ′ (x) (3.14)
The real fields σ ′ and σ ′′ are divergence free, and ε ′ and ε ′′ are compatible fields. This system is associated to the quadratic complex quantity :

W e (x, ω) = σ ′ (x) : ε ′ (x) -σ ′′ (x) : ε ′′ (x) + i σ ′′ (x) : ε ′ (x) + σ ′ (x) : ε ′′ (x) (3.15)
and its effective quantity can be expressed by the two following saddle-point (or min-max) variational principles :

Re We (ω, ε ′ , ε ′′ ) = min ⟨ε ′ ⟩=ε ′ max ⟨ε ′′ ⟩=ε ′′ σ ′ (x) : ε ′ (x) -σ ′′ (x) : ε ′′ (x) (3.16) Im We (ω, ε ′ , ε ′′ ) = min ⟨ε ′ ⟩=ε ′ max ⟨ε ′′ ⟩=ε ′′ σ ′′ (x) : ε ′ (x) + σ ′ (x) : ε ′′ (x) (3.17)
This quadratic saddle-shaped "energy" associated to the system (3.14) is not positivedefinite, thus does not have physical interpretation. Cherkaev and Gibiansky [START_REF] Cherkaev | Variational principles for complex conductivity, viscoelasticity, and similar problems in media with complex moduli[END_REF] make use of the Legendre transform to convert Eq. (3.14) into :

j(x) = D(x)e(x) , with j = ε ′′ (x) σ ′′ (x) , e(x) = -σ ′ (x) ε ′ (x) , ( 3.18)
and where :

D(x) = (L ′′ (x)) -1 (L ′′ (x)) -1 : L ′ (x) L ′ (x) : (L ′′ (x)) -1 L ′′ (x) + L ′ (x) : (L ′′ (x)) -1 : L ′ (x) (3.19)
The convex quadratic quantity W = e.De associated to (3.19) is an actual energy and is naturally interpreted as the dissipated energy averaged over a period of oscillation :

W = ε ′′ : L ′′ : ε ′′ + ε ′ : L ′ : ε ′ ≥ 0 (3.20)
According to the positivity of W, L ′′ must be positive semi-definite, so D is positive semi-definite too. The formulation of the effective constitutive law σ = L : ε can be expressed in the same form as Eq. (3.18) :

j = De, (3.21)
where :

D = ( L′′ ) -1 ( L′′ ) -1 : L′ L′ : ( L′′ ) -1 L′′ + L′ : ( L′′ ) -1 : L′ . (3.22)
The effective energy is given by W = e. De and can be expressed through the following variational principle :

W = min ⟨e ′′ ⟩=e ′′ min ⟨σ ′′ ⟩=σ ′′ ⟨W⟩ (3.23)
The formulation of D as a positive semi-definite quantity related to the convex energy W allows to derive the Voigt and Reuss bounds (the following inequality must be understood in a quadratic sense) :

⟨D -1 (x)⟩ -1 ≤ D ≤ ⟨D(x)⟩ (3.24)
Milton and Berryman [START_REF] Milton | On the effective viscoelastic moduli of two-phase media. II. Rigorous bounds on the complex shear modulus in three dimensions[END_REF] also established the generalization of the Hashin-Shtrikman bounds in the case of viscoelasticity. Based on these variational principles in terms of the convex energy W, [START_REF] Milton | On the effective viscoelastic moduli of two-phase media. I. Rigorous bounds on the complex bulk modulus[END_REF], [START_REF] Milton | On the effective viscoelastic moduli of two-phase media. II. Rigorous bounds on the complex shear modulus in three dimensions[END_REF] and [START_REF] Gibiansky | On the effective viscoelastic moduli of two-phase media. III. Rigorous bounds on the complex shear modulus in two dimensions[END_REF] developed bounds of viscoelastic moduli into the complex plane for a fixed volume fraction. In particular in [START_REF] Milton | On the effective viscoelastic moduli of two-phase media. II. Rigorous bounds on the complex shear modulus in three dimensions[END_REF], Milton and Berryman proposed an algorithm to determine the bounding region in the complex plane of the effective shear modulus of a two-phase material. They showed that the bounding region took the form of a lens-shaped region composed of two circular arcs. More recently, Kern et al. [START_REF] Kern | Tight Bounds on the Effective Complex Permittivity of Isotropic Composites and Related Problems[END_REF] found a tighter lower bound by using variational methods. They were also able to show that the higher bound was almost optimal by funding that some specific points of this higher circular arc were exactly attained by certain microstructures, namely the coated sphere assemblage, the coated cylindrical assemblage, and laminates.

Homogenization with fractional viscoelastic constituents

In [START_REF] Dinzart | Self-consistent approach of the constitutive law of a two-phase visco-elastic material described by fractional derivative models[END_REF], Dinzart and Lipiński address the homogenization of a two-phase polymeric material made of Zener constituents. Based on the correspondence principle, the effective modulus is calculated by means of the self-consistent scheme for a two-phase composite with equal volume fractions [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF]. The mean-field estimate is then compared to experimental results for a mixture of polymers (the polymethyl methacrylate and the styrene acrylonytrile copolymer) obtained through the Dynamical Mechanical Thermal Analysis. Identifying the composite behaviour as a fractional Zener material as well, they successfully showed a good agreement between the homogenization estimate and the experimental loss factors.

More recently, Gallican and Brenner also made use of the correspondence principle for the homogenization of a mixture of fractional Zener constituents [START_REF] Gallican | Homogenization estimates for the effective response of fractional viscoelastic particulate composites[END_REF]. Nonetheless, they adopt a different strategy where advantage has been taken of the purely elastic regimes at short and long times. These authors thus reduced the heterogeneous problem to the resolution of two elastic homogenization problems. The exact asymptotic estimates at long and short times provide exact relations on the effective relaxation spectrum that governs the transient regime. Let us recall the expression of the local fractional Zener constitutive law in phase r :

σ(x, t) + L * : D γ σ(x, t) = L (r) er : ε(x, t) + L (r) eg : L * : D γ ε(x, t), ( 3.25)
where L * = L (r) f : (L (r) eg -L (r) er ) -1 of which the eigenvalues provide the characteristic times (τ (r) ) γ of dimension s γ . The relaxation function of phase (r) reads :

L (r) (t) = L (r) er + G (r) E γ -(t/τ (r) ) γ , with G (r) = L (r) eg -L (r) er (3.26)
while the effective relaxation function in the direct and spectral domain is :

               L(t) = Ler + +∞ 0 G(z)E γ [-(t/z) γ ] dz L(iω) = Ler + +∞ 0 (iωz) γ 1 + (iωz) γ G(z) dz (3.27)
and the decomposition of Lγ (iω) into its real and imaginary terms gives :

               L′ (ω γ ) = Ler + +∞ 0 1 q (iω) γ cos πγ 2 + (iω) γ G(z) dz L′′ (ω γ ) = +∞ 0 1 q (iω) γ sin πγ 2 G(z) dz (3.28)
with q = 1 + 2(ωτ ) γ cos πγ 2 + (ωτ ) γ . These authors approximate the effective spectrum G(z) as a sum of Dirac functions, so the effective relaxation function is considered as given by the series :

L(t) = Ler + k G k E γ [-(t/τ k )] (3.29)
Based on the stationary principles for complex viscoelastic materials (3.16) and (3.17), Gallican and Brenner obtained an approximation of the effective spectrum through the set of 3 equations :

                       K k=1 G k = Leg -Ler K k=1 τ γ k ε : G k : ε = N r=1 c r L (r) f :: ⟨ε (r) ⊗ ε (r) ⟩ (r) K k=1 1 τ γ k ε : G k : ε = N r=1 c r G (r) : (L (r) f ) -1 : G (r) :: ⟨ε (r) ⊗ ε (r) ⟩ (r) (3.30)
where the right-hand terms for the second and third equations are obtained by the exact asymptotic relations developed in [START_REF] Gallican | Exact asymptotic relations for the effective response of linear viscoelastic heterogeneous media[END_REF]. Furthermore, Gallican and Brenner [START_REF] Gallican | Homogenization estimates for polymer-based viscoelastic composite materials[END_REF] solved the system (3.30) for a purely viscoelastic deviatoric and isotropic problem, reducing the determination of G to the shear modulus scalar, so the minimal required number of terms is K = 2. Gallican and Brenner investigated three homogenization schemes : the Mori-Tanaka, the self-consistent and the Generalized Self-Consistent (GSC) estimates. These authors compared these different schemes with full-field computation as reference results, and obtained excellent accuracy with the GSC estimate.

Variational approach : Effective Internal Variable

Although homogenization methods relying on the correspondence principle are efficient, they present two main inconvienients. As mentioned in the introduction, the inverse Laplace transform of the second-moments of the fields can not be obtained. Furthermore, the macroscopic constitutive relations obtained through the correspondence principle display an integral kernel due to the long-memory effects. In the context of structural computations, as mentioned by Badulescu et al. [START_REF] Badulescu | Field statistics in linear viscoelastic composites and polycrystals[END_REF], the use of integral kernels requires the storage of the fields on the whole loading path. Alternatively and to avoid such limitation, one may use incremental variational approaches. This class of methods is based on the step-by-step homogenization of inelastic or non-linear heterogeneous problems through the use of variational principles that take into account the intraphase fluctuations. At each increment, the unknown fields depend only on the solution fields from the previous time step rather than the whole history. Specifically, we get interested in the Effective Internal Variable (EIV) approach developped by Lahellec and Suquet [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF]. Initially introduced for linear viscoelasticity, this method relies on the derivation of a unique incremental potential from the two potentials governing the local behaviour in the GSM framework, presented in Chapter 1. This procedure allows to make use of the variational principles of Ponte Castañeda [3] and [4]. The heterogeneous inelastic fields are approached by an "effective" internal variable. Therefore, the incremental EIV method relies on the introduction of an auxiliary thermoelastic problem that can be solved by means of the classical homogenization methods developed for elasticity.

It is worth mentioning that the same authors later proposed an incremental Rate Variational Procedure (RVP) based on the strain rates rather than strains [START_REF] Lahellec | Effective response and field statistics in elastoplastic and elasto-viscoplastic composites under radial and non-radial loadings[END_REF], allowing to consider local threshold. They successfully obtained an accurate macroscopic response to a tension-compression loading, recovering the Bauschinger effect. Badulescu et al. [START_REF] Badulescu | Field statistics in linear viscoelastic composites and polycrystals[END_REF] then applied the RVP approach to polycrystalline materials and obtained accurate predictions of their effective response to a monotonic loading history, but not for cyclic loadings. Boudet et al. [START_REF] Boudet | An incremental variational formulation for the prediction of the effective work-hardening behavior and field statistics of elasto-(visco)plastic composites[END_REF] also tackled the homogenization of composites with local threshold by extending the EIV to a generalized incremental variational formulation (GIV) that incorporates the elastic unloading. Nevertheless, in the case of linear viscoelasti cconstituents, the efficiency of the EIV model has been shown by Tressou et al. in [110], where they successfully applied it for five different homogenization schemes, coupled to thermic effects [START_REF] Tressou | Contribution à l'homogénéisation des milieux viscoélastiques et introduction du couplage avec la température par extensions d'une approche incrémentale directe[END_REF]. They also extend this method to the case of local phases governed by a generalized Maxwell behaviour, implying multiple internal variables [START_REF] Tressou | Application of the variational EIV approach to linear viscoelastic phases governed by several internal variables -Examples with the generalized Maxwell law[END_REF]. For these reasons, we make use of the EIV method for considering fractional viscoelastic constituents.

Energy and potentials of the fractional damping element

This section presents the thermodynamic formulation of the fractional viscoelastic damping element governed by the Scott-Blair model. For the purpose of the implementation in the EIV method, which considers strain loadings, we focus on the relaxation function. The definition of the potentials relies on the rheological interpretation of the fractional dashpot's relaxation function as the superposition of classical Maxwell elements. By linearity, this representation directly leads to the existence of thermodynamic potentials in the GSM framework. Adopting a discrete point of view, let us consider the superposition of K Maxwell elements, as rheologically represented in Fig. 3.1. The k th Maxwell branch is characterized by the 4 th -order elastic and viscous tensors L k and L k v tensors, with k = 1, ..., K and is associated to the viscous strain α k . The thermodynamic potentials of branch k yield :

w k (ε, α k ) = 1 2 (ε -α k ) : L k : (ε -α k ) φ k ( αk ) = 1 2 αk : L k v : αk (3.31)
where the eigenvalues of the tensors (L k e ) -1 : L k v give the characteristic times. The thermodynamic potentials of the Generalized Maxwell constituent are straighforwardly defined as the superposition of those of each individual element :

w(ε, α 1 , α 2 , ...) = K k=1 w k (ε, α k ) φ( α1 , α2 , ...) = K k=1 φ k (ε, α k ) (3.32)
The constitutive law is given by derivation of the free-energy density, meanwhile the evolution law of each internal variable is given by Biot's equation such that :

σ = ∂w ∂ε (ε, α 1 , α 2 , ...), ∂w ∂α k (ε, α 1 , α 2 , ...) + ∂φ ∂ αk ( α1 , α2 , ...) = 0 ∀ k = 1, ..., K (3.33)
The choice of the number of branches K will be addressed in the third section of this chapter.

Presentation of the EIV method

Incremental variational formulation of the heterogeneous problem

Let us consider the heterogeneous body Ω composed of N phases. The thermodynamic potentials are uniform per phase and are distributed with respect to the indicator function :

w(x, ε, α 1 , α 2 , ...) = N r=1 χ (r) (x)w (r) (ε, α 1 , α 2 , ...) φ(x, ε, α1 , α2 , ...) = N r=1 χ (r) (x)φ (r) ( α1 , α2 , ...) (3.34)
The material is subjected to a strain loading path E(t) on the interval of time t ∈ [0, T ]. The corresponding heterogeneous problem reads :

                                     σ(x, t) = ∂w ∂ε (x, ε, α 1 , α 2 , ...), ∀(x, t) ∈ Ω × [0, T ] ∂w ∂α 1 (x, ε, α 1 , α 2 , ...) + ∂φ ∂ α1 (x, α1 , α2 , ...) = 0, ∀(x, t) ∈ Ω × [0, T ] ∂w ∂α 2 (x, ε, α 1 , α 2 , ...) + ∂φ ∂ α2 (x, α1 , α2 , ...) = 0, ∀(x, t) ∈ Ω × [0, T ] ... div(σ(x, t)) = 0, ∀(x, t) ∈ Ω × [0, T ] ⟨ε(x, t)⟩ = E(t), ∀t ∈ [0, T ] (3.35)
The implicit Euler-scheme is applied to (3.35). The interval of time is divided with time step ∆t = t n+1 -t n . The fields σ n , ε n and α 1,n , α 2,n , ... are supposed to be known at t n . The unknown fields σ n+1 , ε n+1 and α 1,n+1 , α 2,n+1 , ... are obtained by solving the following discrete problem, where the subscript . n+1 is omitted for conciseness :

                                     σ = ∂w ∂ε (ε, α 1 , α 2, , ...), ∀x ∈ Ω ∂w ∂α 1 (ε, α 1 , α 2 , ...) + ∂φ ∂ α1 α 1 -α 1,n ∆t , α 2 -α 2,n ∆t , ... = 0, ∀x ∈ Ω ∂w ∂α 2 (ε, α 1 , α 2 , ...) + ∂φ ∂ α2 α 1 -α 1,n ∆t , α 2 -α 2,n ∆t , ... = 0, ∀x ∈ Ω ... div(σ) = 0, ∀x ∈ Ω ⟨ε⟩ = E, (3.36)
This problem constitutes the Euler-Lagrange equations of the following variational problem :

E(E) = Inf ⟨ε⟩=E Inf α 1 ,α 2 ,... ⟨E(ε, α 1 , , α 2 , ...)⟩, (3.37)
where the incremental potential is given by :

E(ε, α 1 , α 2 , ...) = w(ε, α 1 , , α 2 , ...) + ∆tφ α 1 -α 1,n ∆t , α 2 -α 2,n ∆t , ... . (3.38)
The macroscopic stress at time t n+1 is obtained by derivation of Ẽ :

Σ = ∂ Ẽ(E) ∂E (3.39)
which, obviously, defines the searched constitutive law.

Introduction of the "effective" internal variables and application of the variational method

As the internal variable fields α k,n (x) are not uniform within the phases, the polarization tensor is not piecewise uniform either. As a consequence, the problem (3.36) can not be solved with classical linear homogenization methods. The principle of the present method is to replace the fields α n,k (x) with the auxiliary uniform fields represented by the "effective" internal variables α (r) k,n . The actual effective incremental potential Ẽ is then approximated through a variational method involving the resolution of an auxiliary thermoelastic problem. The solution can obviously not coincide pointwisely with the real solution, though it is required that the first and second moments of these two fields are the same.

Following the variational principle established by Ponte Castañeda [START_REF] Castañeda | New variational principles in plasticity and their application to composite materials[END_REF], the piecewise uniform potential is introduced in the form :

E 0 (x, ε, α 1 , α 2 , ...) = N r=1 χ (r) (x)E (r) 0 (ε, α 1 , α 2 , ...), ( 3.40)
with :

E (r) 0 (ε, α 1 , α 2 , ...) = w (r) (ε, α 1 , α 2 , ...) + ∆tφ (r)   α 1 -α (r) 1,n ∆t , α 2 -α (r) 2,n ∆t , ...   (3.41)
and with :

φ (r)   α 1 -α (r) 1,n ∆t , α 2 -α (r) 2,n ∆t , ...   = 1 2 (α 1 -α (r) 1,n ) : L (r) 0,vα 1 : (α 1 -α (r) 1,n ) + 1 2 (α 2 -α (r) 2,n ) : L (r) 0,vα 2 : (α 2 -α (r) 2,n ) + ... (3.42)
The viscous reference tensors L (r) 0,vα k are uniform per phase. Adding and substracting the auxiliary incremental potential E 0 , the variational problem reads :

Ẽ(E) = Inf ⟨ε⟩=E Inf α 1 ,α 2 ,... ⟨E 0 (ε, α 1 , α 2 , ...) + ∆E(α 1 , α 2 , ...)⟩ , ( 3.43)
where ∆E = E -E 0 . The effective incremental potential is obtained through the variational principle [4] and yields :

Ẽ(E) ≈ Inf ⟨ε⟩=E Inf α 1 ,α 2 ,... ⟨E 0 (ε, α 1 , α 2 , ...)⟩ + Stat α * 1 ,α * 2 ,... ∆E(α 1 , α 2 , ...) (3.44)
The stationarity over the last term provides the equation :

L (r) vα k : α k -α k,n ∆t = L (r) 0,vα k :   α k -α (r) k,n ∆t   (3.45)
The adimensional 4 th -order tensors θ (r) k are introduced :

θ (r) k = L (r) vα k -1 : L (r) 0,vα k (3.46)
They are taken such that each tensor is defined by a single variable θ (r)

k , i.e. θ (r) k = θ (r)
k I, with I the identity tensor. It follows that the internal variables can be expressed as :

α k (x) = α k,n -θ (r) k : α (r) k,n 1 -θ (r) k (3.47)
Replacing α k into Eq. (3.44) by its expression (3.47), the effective incremental potential can be rewritten :

Ẽ(E) ≈ Inf ⟨ε⟩=E Inf α 1 ,α 2 ,... ⟨E 0 ⟩ + N r=1 c r K k=1 θ (r) k 2∆t(θ (r) k -1) (α k,n -α (r) k,n ) : L (r) vg : (α k,n -α (r) k,n ) (3.48) The variables α (r) k,n and θ (r)
k still need to be determined such that they give the optimum value for Ẽ(E). Stationarity of the incremental potential over these variables give the K following systems of 2 non-linear equations :

                       θ (r) k = 1 ± ⟨(α k,n -α (r) k,n ) : (α k,n -α (r) k,n )⟩ (r) ⟨(α k -α (r) k,n ) : (α k -α (r) k,n )⟩ (r) α (r) k,n = ⟨α k,n ⟩ (r) + ⟨α k ⟩ (r) θ (r) k (3.49)
It is noted in [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF] that the plus sign in (3.49) 1 overestimates the energy. It is thus required to keep the minus sign.

Auxiliary thermoelastic heterogeneous problem

The infimum problem (3.48) gives the expression of the internal variables :

α k (x) = θ (r) k ∆t L (r) vα k + L (r) α k α k -1 : L εα k : ε(x) + θ (r) k ∆t L (r) vα k : α (r) k,n (3.50)
Note that, for a superposition of Maxwell elements, it is clear that the coupled elastic tensors are all identical :

L (r) εε = L (r) α k α k = L (r) εα k .
The insertion of (3.50) into the expression of Ẽ0 leads to the definition of a thermoelastic comparison composite :

Ẽ0 (E) = Inf ⟨ε⟩=E ⟨w ∆ (ε)⟩ , w (r) ∆ (ε) = 1 2 ε : L (r) ∆ : ε + τ (r) ∆ : ε + f (r) ∆ (3.51)
The quantity ⟨w ∆ (ε)⟩ = Inf ∆ are uniform in phase r and are expressed by :

                                 L (r) ∆ = L (r) εε - K k=1 L (r) εα k : θ (r) k ∆t L vα k + L (r) α k α k -1 : L (r) εα k τ (r) ∆ = - K k=1 L (r) εα k : θ (r) k ∆t L vα k + L (r) α k α k -1 : θ (r) k ∆t L (r) vα k : α (r) k,n f (r) ∆ = K k=1 θ (r) k 2∆t L (r) vα k : α (r) k,n : θ (r) k ∆t L vα k + L (r) α k α k -1 : L (r) α k α k : α (r) k,n (3.52)
The resolution of the variational problem requires, through the expressions of θ (r) k and α (r) k,n , the determination of the first and second moments of the internal variable α k (x), which itself requires the dermination of the first and second moments of ε(x). The latters are obtained by the resolution of the auxiliary thermoelastic problem through the averaged strain-localization tensor A (r) and its thermal part a (r) :

⟨ε(x)⟩ (r) = A (r) : E + a (r) , (3.53) ⟨ε(x) ⊗ ε(x)⟩ (r) = 1 c r ∂ Ẽ0 ∂L (r) . (3.54)
The effective thermoelastic energy is expressed as :

Ẽ0 = 1 2 E : L : E + τ : E + f , ( 3.55)
with :

                               L = ⟨L ∆ ⟩ + N -1 r=1 c r (L (r) ∆ -L (N ) ∆ ) : (A (r) -I) τ = ⟨τ ∆ ⟩ + N -1 r=1 c r (A (r) -I) T : (τ (r) ∆ -τ (N ) ∆ ) f = ⟨f ∆ ⟩ + N -1 r=1 c r (τ (r) ∆ -τ (N ) ∆ ) : a (r) (3.56)
The heterogeneous thermoelastic problem can be solved using classical methods following linear homogenization schemes or full field computations.

Case of two-phase composites

Following Levin [START_REF] Levin | Thermal Expansion Coefficient of Heterogeneous Materials[END_REF], the localization tensors of a two-phase composite can be expressed as :

A (1) = I + 1 c 1 L (1) ∆ -L (2) ∆ -T : L -⟨L ∆ ⟩ T (3.57) a (1) = 1 c 1 L (1) ∆ -L (2) ∆ -1 : L -⟨L ∆ ⟩ : L (1) ∆ -L (2) ∆ -1 : τ (1)
∆ -τ

∆ and ∆τ = τ

(1) ∆ -τ (2) ∆ : Ẽ0 = ⟨f ⟩+⟨τ ∆ ⟩ : E+ 1 2 E : ⟨L ∆ ⟩ : E+ 1 2 E + (∆L) -1 : ∆τ : L -⟨L ∆ ⟩ : E + (∆L) -1 : ∆τ (3.59)

Discrete representation of fractional viscoelastic behaviours

In link with the EIV method, the relaxation functions of the fractional constituents should be discretized. One could make use of the collocation method, for which its limitations and possible ways of improvements were presented in Chapter 1. This method is surely well-suited for fitting experimental data, or more generally for integral functions in wich the spectrum is not explicitely known. This is not the case for the fractional dashpot, as highlighted in Chapter 1 as well. That is why we rather get interested in a more relevant method based on the discretization of the spectrum. We make use of the discretization procedure developed by Papoulia et al. [START_REF] Papoulia | Rheological representation of fractional order viscoelastic material models[END_REF]. This method consists in rewriting the function as an integral on the whole domain support, and dividing it as three subdomains. The middle one is defined with finite endpoints, so that the function is simply discretized following a composed-midpoint quadrature rule. On the two other subdomains, the function is brought into the spectral domain where the integral can be approximated, then transformed back to the direct domain. The superposition of the three subdomains finally yields the approximate function. However, Papoulia et al. only considered the discretization of the retardation and the relaxation functions of the single fractional dashpot. After presenting the details of their method for this case, we tackle the discretization of the Mittag-Leffler function that is involved in the spectrum of the Maxwell, the Zener and anti-Zener models.

Presentation of the discretization procedure

Creep function of the fractional dashpot

Let us recall the creep function of the Scott-Blair model in one dimension (we consider, say, a purely deviatoric fractional behaviour) :

M(t) = ∞ 0 γ ηΓ(1 + γ)Γ(1 -γ) 1 -e -t/z z 1-γ dz (3.60)
Using the following identities :

Γ(1 + γ) = γΓ(γ), (3.61) Γ(γ)Γ(1 -γ) = π sin(πγ) , (3.62)
the retardance function of the Scott-Blair model (1.80) can be reexpressed as :

M(t) = ∞ 0 sin(πγ) πη 1 z 1-γ 1 -e -t/z dz (3.63)
Let λ and m be the positive endpoints of the interval such that λ < m. Let us recall the notation K for denoting the number of Kelvin units in the discrete model. The discretization procedure is as follows : a change of variable z = e θ is introduced so that the function is defined on (-∞, +∞) and reads :

M(t) = sin(πγ) πη ∞ -∞
e γθ 1 -e -t/e θ dθ.

(3.64)

The midpoint-based approximation is applied on the interval [Λ, Π], where Λ = ln(λ) and Π = ln(m). The precision can be increased by approximating the integral on the truncated portions (-∞, Λ] and [Π, +∞) in the frequency domain, then transforming it back to the direct domain.

On the portion (-∞, Λ]

The Fourier transform of the function (3.63) with respect to the characteristic times z is expressed by :

M * -∞,Λ (ω) = sin(πγ) πη e Λ 0 1 (1 + iωz)z 1-γ dz. ( 3.65)
If the upper bound e Λ is sufficiently low, the range of frequency of interest satisfies ωz ≪ 1.

Making the reasonable assumption ω ≪ e -Λ , the denominator reduces to 1 + iωz ≈ 1.

Note that this hypothesis can be justified a posteriori (we do not provide further details here, but they can be found in [START_REF] Papoulia | Rheological representation of fractional order viscoelastic material models[END_REF]). The integral can be found in closed form :

M * -∞,Λ (ω) ≈ sin(πγ) πηγ e αΛ , ( 3.66)
It follows that the inverse transform is given by :

M -∞,Λ (t) ≈ sin(πγ) πηγ e αΛ , (3.67)

On the portion [Λ, Π]

The change of variable z = e θ leads to the following expression :

M(t) = sin(πγ) πη Π Λ
e γθ 1 -e -t/e θ dθ (3.68)

The interval [Λ, Π] is subdivided in K -1 equal subintervals. The compound trapezoidal method is applied such that the function is estimated on the middle of each subinterval :

M N,Λ,Π (t) = sin(πγ) πη (Π -Λ) N K-1 k=1 e γz k 1 -e -t/e z k (3.69) z k = Λ + (k -1/2) (Π -Λ) (K -1) (3.70)

On the portion [Π, +∞)

The Fourier transform of the integral (3.63) on the interval [Π, +∞) reads :

M * Π,+∞ (ω) = sin(πγ) πη +∞ e Π 1 (1 + iωz)z 1-γ dz (3.71)
Here, the assumption of high values of ω is made [START_REF] Papoulia | Rheological representation of fractional order viscoelastic material models[END_REF], such that ω ≫ e -Π . The approximation 1 + iωz ≈ -iωz can thus be made. The assumption leads to the function :

M * Π,+∞ (ω) ≈ - sin(πγ) πη +∞ e Π 1 iωz 2-γ dz (3.72)
can be calculated in closed form :

M Π,+∞ (t) ≈ sin(πγ) πη e γΠ (1 + iωe Π )(1 -γ) . ( 3.73)

Complete expression

The complete approximation M K (t) is finally given by :

M K (t) = sin(πγ) πη e γΛ γ + e γΠ 1 -γ (1 -e -t/e Π ) + (Π -Λ) N K-1 k=1 e γz k 1 -e -t/e z k (3.74)
The discrete creep function of the fractional Maxwell model can be straightforwardly obtained by the addition of the compliance modulus M eg :

M K,Maxwell (t) = M eg + M K (t) (3.75)

Relaxation function of the fractional dashpot

The same procedure is followed for the relaxation function of the fractional dashpot, of which we also recall the expression in one dimension :

L(t) = ∞ 0 η Γ(1 -γ)Γ(γ) 1 z 1+γ e -t/z dz (3.76)
Once again using the identities (3.61) and (3.62), the relaxation function reads :

L(t) = ∞ 0 ηsin(πγ) πz 1+γ e -t/z dz (3.77)
With the change of variable z = e θ , it can also be expressed as the following :

L(t) = ηsin(πγ) π +∞ -∞
e -t/e θ e -γθ dθ (3.78)

On the portion (-∞, Λ]

The Fourier transform of the function (3.77) reads :

L * -∞,Λ (ω) = ηsin(πγ) π e Λ 0 iω z γ (1 + iωz) dz (3.79)
Assuming ω ≪ e -Λ , the approximation 1 + iωz ≈ 1 + iωe Λ can be made [START_REF] Papoulia | Rheological representation of fractional order viscoelastic material models[END_REF]. The integral can then be calculated in closed form in the frequency domain,

L * -∞,Λ (ω) ≈ ηsin(πγ) π iωe (1-γ)Λ (1 + iωe Λ )(1 -γ) , ( 3.80)
and its expression in the time domain is :

L -∞,Λ (t) ≈ ηsin(πγ) π e -t/e Λ e γΛ (1 -γ) (3.81)

On the portion [Λ, Π]

Using the same discretization as for the creep function, the Eq. (3.78) on the interval

[Λ, Π] is approximated on each point z k = Λ + (k -1/2) (Π -Λ) (N -1) by : L N,Λ,Π (t) = ηsin(πγ) π (Π -Λ) (N -1) K-1 k=1 e -t/e z k e -γz k , ( 3.82)
On the portion [Π, +∞)

In the frequency domain, Eq. (3.77) reads :

L * Π,+∞ (ω) = ηsin(πγ) π +∞ e Π iω z γ (1 + iωz) dz (3.83)
Assuming ω ≫ e -Π , the term 1 + iωz is approximated by iωz [START_REF] Papoulia | Rheological representation of fractional order viscoelastic material models[END_REF] and yields :

L * Π,+∞ (ω) ≈ ηsin(πγ) π e -γΠ γ , (3.84)
which expression is the same in time domain as it doesn't depend on ω.

Complete expression

The complete relaxation function thus reads :

L K (t) = ηsin(πγ) π e -t/e Λ e γΛ (1 -γ) + e -γΠ γ + K-1 k=1 e -t/e z k e -γz k (3.85)
The relaxation function of a fractional Kelvin-Voigt model (KV) is directly obtained by adding the relaxed elastic modulus L er :

L N,Kelvin-Voigt (t) = L er + L N (t) (3.86)

Comparison with the collocation method

For the sake of illustration, let us consider an arbitrary fractional viscoelastic damping element with exponent γ = 0.5, elasticity µ = 1 GPa and relaxation time τ σ = 1 s. The relaxation function of this element is discretized following the collocation method and the improved midpoint-based method for comparison. Let's take the discretization parameters K = 7 and z ∈ [10 -3 , 10 3 ] s, where the characteristic times are distributed with a logarithmic scale. Fig. 3.2 shows that the instantaneous and the long-time responses are accurately reproduced by both approximation methods. However, the midpoint-based method provides a better accuracy in the transient regime, even for a low number of branches. Therefore, the improved midpoint-based method will be kept in the context of this work. 0.0 0.2 0.4

t [s] L(t)
exact midpoint-based collocation

New application of the discretization procedure

The discretized retardance function of the fractional Maxwell constituent and the discretized relaxation function of the Kelvin-Voigt constituent can be expressed from the discretization of the fractional dashpot element. However, as seen in Chapter 1, their dual functions, as well as the Zener's and anti-Zener's, involve the Mittag-Leffler function. Therefore, they can not be directly obtained by a simple combination of Eqs. (3.74) or (3.85) with elastic coefficients. The integral definition of the Mittag-Leffler function displays a rational fraction spectrum with polynoms of non-integer order. In this work, we apply the previous discretization procedure to this particular function. The creep or relaxation functions of the different assemblies can then be completed by adding the asymptotic moduli.

Discretization of the Mittag-Leffler function

As seen in section 1.4.2, the Mittag-Leffler function is defined by the power series :

E γ (z) = +∞ n=0 z n Γ(1 + γn) , ∀z ∈ C, 0 < γ ≤ 1, (3.87)
Alternatively, this function can be represented in its integral form :

E γ (-t γ ) = +∞ 0 H γ (ξ)e -t/ξ dξ, t > 0 with 0 < γ ≤ 1 (3.88)
with

H γ (ξ) = 1 π ξ γ-1 sin(γπ) 1 + 2ξ γ cos(γπ) + ξ 2γ (3.89)
In the following, the denominator will be noted for conciseness :

P γ (ξ) = 1 + 2ξ γ cos(γπ) + ξ 2γ (3.90)
Reexpressing Eq. (3.88) with the argument of interest -(t/τ ) γ and applying the change of variable z = τ ξ, the Mittag-Leffler function reads :

E γ [-(t/τ ) γ ] = sin(γπ) πτ γ +∞ 0 z γ-1 P γ (z/τ ) e -t/z dz (3.91)
As well as the previous section, use will be made of the change of variable z = e θ :

E γ [-(t/τ ) γ ] = sin(γπ) πτ γ +∞ -∞ e θγ P γ (e θ /τ ) e -t/e θ dθ (3.92)
The integral we consider is therefore given by :

I(t) = +∞ -∞ e θγ P γ (e θ /τ ) e -t/e θ dθ (3.93)
For recall, we let Λ = ln(λ) and Π = ln(m) be the endpoints of the interval [Λ, Π] that is subdivided into K -1 equal subintervals. On the portion [Λ, Π], we use the compound midpoint rule. On the ends (-∞, Λ] and [Π, +∞), the integral is evaluated in the frequency domain, then transformed back to the direct domain.

On the portion (-∞, Λ]

The Fourier transform of (3.91) on the interval (-∞, Λ] reads :

I * -∞,Λ (ω) = e Λ 0 z γ P γ (z/τ ) iω iωz + 1 dz (3.94)
Let's assume ω ≪ e -Λ for the range of frequencies of interest. Arguments for this hypothesis will be provided later. Thus, we have iωz + 1 ≈ iωe Λ + 1, so using identities on Fourier transforms, the function is given in time domain by :

I -∞,Λ (t) = e -Λ e -t/e Λ e Λ 0 z γ P γ (z/τ ) dz, (3.95)
where the integral can be computed numerically.

On the portion [Λ, Π]

The compound midpoint-based approximation straightforwardly gives :

I Λ,Π (t) = (Π -Λ) (K -1) N -1 k=1 e γz k P γ (e z k /τ ) e -t/e z k , (3.96)
where

z k = Λ + (k -1/2) (Π -Λ) (K -1
) .

On the portion [Π, +∞)

In the frequency domain, the function reads :

I * Π,+∞ (ω) = +∞ e Π z γ P γ (z/τ ) iω iωz + 1 dz (3.97)
Assuming ω ≫ e -Π , we have iωz + 1 ≈ iωz. Again, this assumption will be motivated later. Moreover, applying the change of variable z γ = X, the function in the time domain reduces to :

I Π,+∞ (t) = +∞ e γΠ 1 P γ (X/τ ) dX, (3.98)
which can be determined analytically :

I Π,+∞ (t) = τ 2γ γ C -B 2 /4 π 2 -arctan e γΠ + B/2 C -B 2 /4 , ( 3.99)
with B = 2τ γ cos(γπ) and C = τ 2γ .

Complete expression

The complete expression reads :

E γ [-(t/τ ) γ ] = sin(γπ) πτ γ e -t/e Λ e Λ e Λ 0 z γ P γ (z/τ ) dz + Π -Λ K -1 K-1 k=1 e γz k P γ (e z k /τ ) e -t/e z k + τ 2γ γ C -B 2 /4 π 2 -arctan e γΠ + B/2 C -B 2 /4 , (3.100)
which is interpreted as a superposition of N Maxwell elements. Its expression in the frequency domain reads :

FT (E γ [-(t/τ ) γ ]) (ω) = sin(γπ) πτ γ e -Λ iω iω + e Λ e Λ 0 z γ P γ (z/τ ) dz+ Π -Λ K -1 K-1 k=1 e γz k P γ (e x k /τ ) iω iω + e z k + τ 2γ γ C -B 2 /4 π 2 -arctan e γΠ + B/2 C -B 2 /4 (3.101)

Justification of the assumptions

On the subintervals (-∞, Λ] and [Π, +∞), we assumed that the pulsation was respectively very low and very high compared to the inverse upper and lower bounds. This hypotheses entail the decreasing value of I * (ω) on these portions as the bounds are large. This section aims at justifying these assumptions a posteriori, by showing that the approximated function is bounded.

Let us decompose I * (ω) = T 1 + T 2 + T 3 , where T 1 , T 2 and T 3 are the integrals over (-∞, Λ], [Λ, Π] and [Π, +∞) respectively. Their integrand is noted q(θ) :

q(θ) = e θγ P γ (e θ /τ ) iωe θ 1 + iωe θ (3.102)
For any complex function q(θ), the highest value of the norm |I * (ω)| is bounded by the expression [START_REF] Papoulia | Rheological representation of fractional order viscoelastic material models[END_REF] :

|I * (ω)|≥ 1 √ 2 (t 1 + t 2 + t 3 ), (3.103)
where

t 1 = Λ -∞ |q(θ)|dθ, t 2 = Π Λ |q(θ)|dθ, t 3 = +∞ Π |q(θ)|dθ (3.104)
The rheological approximation we assessed is given by I * Λ,Π,N = T1 + T2 + T3 , where :

3 EIV approach for fractional viscoelastic constituents

91 T1 = Λ -∞ e θγ P γ (e θ /τ ) iωe θ 1 + iωe Λ dθ T2 = Π -Λ K -1 K-1 k=1 e γz k P γ (e z k /τ ) e -t/e z k T3 = +∞ Π e θγ P γ (e θ /τ ) dθ (3.105)
The approximated function I * Λ,Π,N proves to converge to the exact expression

I * if |I * (ω) - I *
Λ,Π,N | admits a finite limit.

For T 1 |T 1 -T1 | ≤ Λ -∞ e θγ P γ (e θ /τ ) iωe θ 1 + iωe θ - iωe θ 1 + iωe Λ dθ = ω Λ -∞ e θγ P γ (e θ /τ ) 1 1 + iωe Λ (e Λ -e θ )iωe θ 1 + iωe θ dθ ≤ ωe Λ Λ -∞ e θγ P γ (e θ /τ ) iωe θ 1 + iωe θ dθ = ωe Λ t 1
The error is minimized if ωe Λ ≪ 1, i.e for ω ≫ e -Λ .

For T 2

The error of the compound midpoint quadrature approximation over interval I subdivided in an infinite number of intervals of length h → 0 is bounded by :

h 2 I |q ′′ (θ)/24| dθ, (3.106)
where h = Π -Λ K -1 . The second derivative of the integrand is expressed as :

q ′′ (θ) = q(θ) - Υ(θ) 1 + iωe θ - ∂ 2 ∂θ 2 ln P γ (e θ /τ ) + (C γ (θ) -Υ(θ)) 2 , ( 3.107)
with

C γ (θ) = 1 + γ - ∂ ∂θ (ln P γ (e θ /τ )) and Υ(θ) = iωe θ 1 + iωe θ .

92

Discrete representation of fractional viscoelastic behaviours

We find an upper bound of |q ′′ (θ)| by bounding the quantity :

|U | = - Υ(θ) 1 + iωe θ - ∂ 2 ∂θ 2 ln P γ (e θ /τ ) + (C γ (θ) -Υ(θ)) 2 ≤ Υ(θ) 1 + iωe θ + ∂ 2 ∂θ 2 ln P γ (e θ /τ ) + (C γ (θ) -Υ(θ)) 2
Let's evaluate each term. Using the inequality |1 + iωe θ |≥ |iωe θ |, the first term is straightforwardly bounded :

Υ(θ) 1 + iωe θ ≤ 1
For the second term, specifying the derivatives :

∂P γ (e θ /τ ) ∂θ = P ′ γ (e θ /τ )×(e θ /τ
) and

∂ 2 P γ (e θ /τ ) ∂θ 2 = P ′′ γ (e θ /τ )×(e θ /τ ) 2 +P ′ γ (e θ /τ )×(e θ /τ ) ,
it can be reexpressed as :

∂ 2 ∂θ 2 ln P γ (e θ /τ ) = P ′′ γ (e θ /τ ) × (e θ /τ ) 2 + P ′ γ (e θ /τ ) × (e θ /τ ) P γ (e θ /τ ) + (P ′ γ (e θ /τ ) × (e θ /τ )) 2 (P γ (e θ /τ )) 2
Upper bounds must be found for |P ′ γ (e θ /τ )| and |P ′′ γ (e θ /τ )|, meanwhile we must evaluate a lower bound for |P γ (e θ /τ )|.

P ′ γ (e θ /τ ) = γ(e θ /τ ) -1 P γ (e θ /τ ) -1 + (e θ /τ ) 2γ ≤ γ e Π-Λ (2 + e γΠ /τ ) + e 2Πγ-Λ τ 2γ-1
= C (2) , and

P ′′ γ (e θ /τ ) = - γ (e θ /τ ) 2 P γ + γ(e θ /τ ) -1 P ′ γ - γ (e θ /τ ) 2 (e θ /τ ) 2γ -1 + 2γ 2 (e θ /τ ) 2(γ-1) ≤ γ (e Λ /τ ) -2 C (1) + (e Λ /τ ) -1 (C (2) -1 + (e Π /τ ) 2γ ) + 2γ(e Λ /τ ) 2(γ-1)
= C (3) , where C (1) is the upper bound of P γ (e θ /τ ) :

3 EIV approach for fractional viscoelastic constituents 93 (1) .

P γ (e θ /τ ) ≤ 1 + e θ τ γ 2cos(γπ) + e θ τ γ = C
We also have to find a lower bound for the denominator. P γ is defined for 0 < γ ≤ 1. Let's first consider γ < 1, which is the case of interest for this work. The norm can thus be bounded :

|P γ (e θ /τ )|≥ 1 -cos(γπ)
A bound can also be established in the case γ = 1 (although in fractional visoelasticity, γ is strictly lower than 1). To this end, we have to analyse the ratios P ′ γ /P γ and P ′′ γ /P γ . As P ′ γ is a function of P γ , and considering the value of interest θ = ln(τ ), then P ′ γ /P γ = γ. Similarly, the ratio P ′′ γ /P γ can be bounded.

We finally provide the upper bound of the last term :

(C γ (θ) -Υ(θ)) 2 ≤ 2C 2 γ (θ) + 2Υ 2 (θ),
where :

|C γ (θ)|≤ 1 + γ + C (2) (e Π /τ ) 1 -cos(γπ)
So, we finally bounded :

|T 2 -T2 |≤ |U | maj Π -Λ K 2 t 2 , ( 3.108)
|U | maj = 5 + C (3) (e Π /τ ) 2 + C (2) (e Π /τ ) 1 -cos(γπ) + C (2) (e Π /τ ) 2 1 -cos(γπ) 2 + 2 γ + C (2) (e Π /τ ) 2 1 -cos(γπ) (3.109) For T 3 |T 3 -T3 | ≤ +∞ Π e θγ P γ (e θ /τ ) iωe θ 1 + iωe θ -1 dθ = +∞ Π e θγ P γ (e θ /τ ) 1 1 + iωe θ iωe θ iωe θ dθ ≤ ω -1 e -Π +∞ Π e θγ P γ (e θ /τ ) iωe θ 1 + iωe θ dθ = ω -1 e -Π t 3
Thus the error is minimized if ω ≫ e -Π , justifying a posteriori the assumptions that were made for finding the expressions (3.95) and (3.98).

Conclusion

This chapter first presented mean-field homogenization methods for linear viscoelasticity. The fundamental correspondence principle allows to apply the classical homogenization schemes on the heterogeous symbolic elastic problem in the Laplace-Carson domain. Considering a purely imaginary Laplace-Carson variable, the problem is specified for a harmonic loading. Use can therefore be made of variational principles with complex fields, for which case the variational methods take the form of a saddle-point min-max problems. We also mentioned the successful use of these approaches to obtain homogenized estimates for fractional viscoelasticity in the spectral domain. However, the transformation of the averaged fields in the direct domain can generally not be obtained exactly. Besides, the second-order moments of the fields are inaccessible, prohibiting the use of sophisticated estimates. We therefore take interest in the incremental variational approach developed by Lahellec and Suquet [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF]. Their EIV method relies on the derivation of a unique incremental potential, obtained from the dissipation and the free-energy potentials in the framework of the GSM. The difficulty of non uniform viscous fields is tackled with the introduction of a comparison thermoelastic composite. Use is then made of the variational principles of Ponte Castañeda [START_REF] Castañeda | New variational principles in plasticity and their application to composite materials[END_REF] and [4] which take into account the second-moments of the fields and the fields fluctuations. The statistics of the fields are obtained by estimating the error between the real incremental potential and the fictitious one. This method can be implemented by solving an implicit system of non-linear equations that involves the homogenization of a thermoelastic composite. The specification of the type of microstructure suggests the choice of the thermoelastic homogenization scheme.

The objective of this chapter was twofold. On the one hand, we presented the Effective Internal Variational method developed by [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF] in the context of multiple internal variables as generalized by Tressou et al. [START_REF] Tressou | Application of the variational EIV approach to linear viscoelastic phases governed by several internal variables -Examples with the generalized Maxwell law[END_REF]. Then, in the perspective of incorporating fractional laws into the incremental method, we presented the discretization procedure developed

Introduction

In this chapter, we address two different homogenization problems. The first one concerns polymer-based composite materials. Polymeric materials exhibit viscoelastic damping properties, wether they are natural or synthetic. Their mechanical behaviour has been investigated since the early XX th century for the design of structures and vehicles [START_REF] Jones | Handbook of Viscoelastic Vibration Damping[END_REF]. The viscoelasticity of such materials can be accurately described by means of fractional models. For example Rouleau et al. [START_REF] Rouleau | Characterization and Modeling of the Viscoelastic Behavior of a Self-Adhesive Rubber Using Dynamic 4 Bibliography Mechanical Analysis Tests[END_REF] characterized a rubber by means of DMA, and modelled its behaviour as that of a fractional Zener material, while Baudis et al. [START_REF] Baudis | Thorough ultrasonic rheology of soft, visco-elastic materials: Example of crosslinked Polyurethane elastomer[END_REF] measured the properties of a polyurethane at ultrasonic frequencies, and identified a behaviour in agreement with a fractional Kelvin-Voigt law. In polymer-based composites, the interactions between the constituents significantly influences the damping properties, see for example [115] [116] who consider fibrous composites, or [START_REF] Chabert | Filler-filler interactions and viscoelastic behavior of polymer nanocomposites[END_REF] and [START_REF] Baudis | Caractérisation acoustique de matériaux viscoélastiques homogènes et hétérogènes pour la furtivité sous-marine[END_REF] for a viscoelastic matrix reinforced with stainless steel particles. In this chapter, we precisely investigate the latter type of composites. Following Gallican and Brenner [START_REF] Gallican | Homogenization estimates for the effective response of fractional viscoelastic particulate composites[END_REF], we consider a fractional Zener matrix reinforced with spherical rigid particles. As mentioned in Chapter 3, these authors established homogenization estimates by means of the correspondence principle and complex variational principles, based on the Generalized Self-Consistent (GSC) scheme. They obtained accurate results with reference FFT-based computations.

The second case deals with polycrystalline sea ice. The knowledge of ice mechanics is essential for addressing ice engineering problems. For example, the safe design of offshore structures on moving ice sheets, or the use of ice-breakers [START_REF] Sanderson | Ice Mechanics and Risks to Offshore Structures[END_REF]. In earth science, the evolution of sea ice sheets is an essential concern in the context of its interactions with climate change [START_REF] Sun | Evolution of the Global Coupled Climate Response to Arctic Sea Ice Loss During 1990-2090 and its Contribution to Climate Change[END_REF]. The mechanics of sea ice involves a wide range of length and time scales. At large length-scales, significant variability of the mechanical fields are observed within a few months [START_REF] Weiss | Linking scales in sea ice mechanics[END_REF]. This is not the case at the local scale, where the evolution of the fields are caused by the microstructure heterogeneities [START_REF] Jellinek | Viscoelastic Properties of Ice[END_REF]. In [START_REF] Sinha | Rheology of columnar-grained ice: Phenomenological viscoelasticity of columnar-grained ice has been experimentally investigated and analyzed in an effort to explain many apparent peculiarities of the deformation behavior of polycrystalline ice[END_REF], Sinha gave attention to the inital and transient creep regime of sea ice occurring at short times (compared to the few months mentioned earlier). In his article, he carried out uniaxial compressive creep loading experiments on columnar-grained ice. He showed that the material displayed an instantaneous elastic response followed by a delayed viscoelastic strain. He provided a detailed characterization of the microstructure, however he proposed only a macroscopic stress-strain relationship that didn't consider the polycristalline nature of the material. Based on the detailed description of the microstructure, Schapery tackled this polycrystalline ice problem with a micromechanical approach [START_REF] Schapery | Linear Elastic and Viscoelastic Deformation Behaviour of Ice[END_REF]. His model successfully recovers the experimental results obtained by Sinha. Besides, he was able to predict the expected power-law behaviour at long times that is specific to certain polycrystalline materials, mostly metal alloys [START_REF] Da | On The Viscous Flow in Metals, and Allied Phenomena[END_REF]. We take interest in this problem because of the particular assumption made by Schapery, that the local viscoelastic slip system follows a fractional Maxwellian behaviour.

First of all, we provide insights of our numerical implementation, and present an example as a validation case. After that, we apply the Effective Internal Variable method for approaching the two abovementioned problems. For the polymer-based composite, we compare our results to those of the full-field computations carried out in [START_REF] Gallican | Homogenization estimates for the effective response of fractional viscoelastic particulate composites[END_REF]. The matrix follows a fractional Zener behaviour, of which the relaxation function is discretized by means of the improved midpoint-based procedure [START_REF] Papoulia | Rheological representation of fractional order viscoelastic material models[END_REF] presented in the previous chapter. The composite is subject to harmonic loadings. We then address the homogenization of polycrystalline sea ice. The same discretization procedure is applied to the local fractional Maxwell behaviour. For monitoring stress loading control in the EIV, additional iterative instructions are implemented in the algorithm. Our EIV simulations are confronted to the experimental results from [START_REF] Sinha | Rheology of columnar-grained ice: Phenomenological viscoelasticity of columnar-grained ice has been experimentally investigated and analyzed in an effort to explain many apparent peculiarities of the deformation behavior of polycrystalline ice[END_REF].

Numerical validation for a matrix-inclusion microstructure 4.1.1 EIV algorithm and numerical implementation

Preliminary remarks on the EIV method (a) Numerical limitations

The resolution of the non-linear implicit system (3.49) can cause numerical issues. Badulescu et al. [START_REF] Badulescu | Field statistics in linear viscoelastic composites and polycrystals[END_REF] suggest to rewrite the expression of θ (r) as a polynom, in order to avoid prohibited values inside the square root :

θ (r) k -1 2 = ⟨(α k,n -α (r) k,n ) : (α k,n -α (r) k,n )⟩ (r) ⟨(α k -α (r) k,n ) : (α k -α (r) k,n )⟩ (r) (4.1)
They also propose a 2-step resolution in the context of the RVP, found in Appendix B of their article, for a two-phase composite and for a unique internal variable. Because of the similarity of the equations between the RVP and the EIV methods, this can be applied to the latter, as done by Boudet et al. [START_REF] Boudet | An incremental variational formulation for the prediction of the effective work-hardening behavior and field statistics of elasto-(visco)plastic composites[END_REF]. These authors consider only one internal variable (K = 1). Let us omit the subscript . k for the case of a single internal variable. The first step consists in eliminating α (r) n in Eq. (3.49) 2 by substituting the term ⟨α⟩ (r) with its expression provided by Eq. (3.50). We rewrite the latter equation by substituting the expressions of the coefficients with B (r) 1 and B (r) 2 for conciseness :

⟨α(x)⟩ (r) = B (r) 1 : ⟨ε(x)⟩ (r) + B (r) 2 : α (r) n (4.2)
The term ⟨ε(x)⟩ (r) is itself linear to α (r) n through the expression of the auxiliary polarization tensor (3.52) 2 . Let us also rewrite the average strain by replacing the complete expressions of the coefficients with B (r) 3 and C (r) n :

⟨ε⟩ (r) = A (r) : E + a (r) = A (r) : E + B (r) 3 : ∆τ (r) = A (r) : E + B (r) 3 : C (r)
n α (r) n 100

Numerical validation for a matrix-inclusion microstructure

The effective internal variable α (r) n can thus be written in terms of known quantities through the linear relation :

⟨α(x)⟩ (r) = B (r) 1 : E + B (r) 1 : B (r) 3 C (r) n + B (r) 2 : α (r) n (4.3)
After incorporating the last expression in Eq. (4.1), the resolution finally reduces to the minimization of a unique function f (θ (r) ) that depends on θ (r) only. However, the rewriting of the α (r)

n according to the two-step resolution yields a rational fraction with several poles. Still in [START_REF] Badulescu | Field statistics in linear viscoelastic composites and polycrystals[END_REF], it has been shown that f (θ (r) ) can be regularized by multiplying both sides with the denominator. Finally, although the expression (4.1) seems more convenient, its resolution can lead to two solutions. These authors also point out that there is no theoretical criterion for discriminating the physical solution. They empirically suggest that the smallest positive value of θ (r) is generally the best choice, as was already suggested in [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF].

However, note that the formulation of the system as a unique objective function is not straightforwardly applicable for a generalized Maxwell constituent. Indeed, the expression of α k with respect to the polarization field τ Idiart et al. provide an evaluation of the EIV method in the context of linear viscoelasticity in [START_REF] Martín | Model reduction by meanfield homogenization in viscoelastic composites. I. Primal theory[END_REF] and [START_REF] Martín | Model reduction by meanfield homogenization in viscoelastic composites. II. Application to rigidly reinforced solids[END_REF]. In the first paper, they expose the mathematical structure of this method by deriving a reduced-model relying on a Cauchy-Schwarz inequality rather than Legendre transforms. This reformulation allows to highligt that, meanwhile the dissipation potential is always convex, it is not the case for the free-energy. Indeed, whenever the intraphase fluctuations rate of the inelastic strain field are negative, the free-energy is no longer convex. Also, they addressed the issue of multiplicity of the solutions and prescribed to choose the root that minimizes the free-energy density. This has been done by specifying that it is not restrained to a positive value, contrary to the previous suggestions in litterature. The consequences of the loss of convexity are investigated in the second paper [START_REF] Martín | Model reduction by meanfield homogenization in viscoelastic composites. II. Application to rigidly reinforced solids[END_REF] in which the authors investigated the case of loading-unloading histories of deformation. For a history of loading where the applied deformation vanishes, the EIV method displays spurious estimates. For oscillating loadings without vanishing total strain, the EIV method provides estimates in good agreement with reference results. The author stress out that this accuracy is contradictory with the periodic loss of convexity of the free-energy density, and suggest that the conditions on the inelastic and total strain fluctuations may not be accurate.

Initialization

Inputs :

-Time vector -Loading path E(t) -Material properties -Fields at time : -1st and 2nd moments of -1st and 2nd moments of -1st and 2nd moments of

L (r) , c r t = 0 ε n=0 α k,n=0 σ n=0
Class « Linear media of comparison » Input :

Method : Equations of the thermoelastic media Output: For validation of the numerical implementation, we consider the same example as in [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF], which is a matrix-fibers microstructure such that the linear effective properties are exactly given by the lower Hashin-Shtrikman bound. The constituents are isotropic, incompressible, and they are characterized by the shear moduli µ (r) . The matrix is denoted as phase r = 2 and the inclusions as r = 1. The fibers are oriented in the e 3 direction such that the composite is transversely isotropic and is characterized by the unique effective shear modulus μHS , as shear strain loadings will be considered : μHS = µ (1) µ (1) + µ (2)

L (r) , c r , θ (r) k , α (r) k,n L (r) Δ , τ (r) Δ , f ( r)
α (r) k,n ⟨α⟩ (r) , ⟨α ⊗ α⟩ (r) ⟨α⟩ (r) , ⟨α ⊗ α⟩ (r) K θ (r) k , α (r) k,n
c 1 µ (1) + c 2 µ (2) µ (2) + c 2 µ (1) + c 1 µ (2) (4.5)
Now, consider the elastic inclusions embedded in a linear viscoelastic matrix following a Maxwellian behaviour. The latter is characterized by the viscosity and elastic moduli η (2) and µ (2) respectively. The viscoelastic effective shear modulus can be calculated through the correspondence principle. The shear modulus of the matrix reads in the direct and the Laplace-Carson domain :

µ (2) (t) = µ (2) e -tµ (2) /η (2)
μ(2) (p) = pµ (2) p + µ (2) /η (2) (4.6)

The effective shear modulus in the Laplace-Carson domain thus reads :

μHS (p) = µ (1) µ (1) + μ(2) (p) c 1 µ (1) + c 2 μ(2) (p) μ(2) (p) + c 2 µ (1) + c 1 μ(2) (p) (4.7)
We consider the same material data as in [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF] :

2µ (1) = 266.6 GPa 2µ (2) = 46.6 GPa and η (2) = 10 GPa (4.8)

as well as different shear strain loadings E(t) along the path t ∈ [0, T ] such that :

E(t) = E 11 (t)(e 1 ⊗ e 1 -e 2 ⊗ e 2 ) ∀ t ∈ [0, T ] (4.9)
The effective stress response can be analytically found by calculating the inverse Laplace-Carson transform -using Mathematica for example -of the effective constitutive law : Σ(p) = μHS (p) Ê(p).

Let us consider the constitutive law of the generalized Maxwell :

µ (2) (t) = µ ∞ + K k=1 µ (2) k e -tµ (2) k /η (2) k μ(2) (p) = K k=1 pµ (2) k p + µ (2) k /η (2) k (4.10)
As a validation for the extension of the algorithm to K internal variables, let us consider the trivial case of identical branches with the same characteristic times τ

(2)

k = η (2) /µ (2) k such that µ (2)
k = µ (2) /K and η

(2) k = η (2) /K, ∀k = 1, ..., K, and the additional spring is removed, i.e. µ ∞ = 0. With this choice, the constitutive law of the superposition of K Maxwell elements (4.10) reduces to the simple Maxwell element (4.6). In this way, the analytical solutions derived for the simple Maxwellian matrix also stand for the formulation with the corresponding K internal variables.

Validation for monotone loadings

Let us consider the volume fraction of fibers c 1 = 0.45, and fix K = 10 the number of Maxwell branches. We first consider monotonous loadings compared with the analytical results obtained through the correspondence principle. We take the time step ∆t = 5×10 -3 s. Note that it is taken as ∼ 85 times smaller than the value of the relaxation time τ (2) = 0.43s, so that the transient regime can be well captured. Let us consider a relaxation loading. A shear strain loading is instantaneously applied at t = 0 and is then held constant : and the inelastic fluctuations read :

E(t) = 0 if t < t 0 E 11 (e 1 ⊗
C (r) (α k ) = ⟨α k : α k ⟩ (r) -⟨α k ⟩ (r) : ⟨α k ⟩ (r) (4.13)
As mentioned above, all the branches are identical so that the inelastic fluctuations of every viscous strain are all equal. Note also that, as expected by the Hashin-Shtrikman estimate, the fluctuations in the fibers are null.

Constant strain rate loading

Let us consider the constant strain rate Ė(t) such that :

E(t) = bt(e 1 ⊗

Validation for harmonic loading and remarks

The material with a Hashin-Shtrikman microstructure is here subjected to the following harmonic loading : Figs. 4.6 show once again the fields provided by the EIV compared to the exact results. In accordance with the exact predictions, the effective stress displays in the first 2 seconds a severe attenuation, followed by a stabilized oscillating regime. A close-up is represented in Fig. 4.7 to show details. On the transient attenuation regime, the accuracy of the EIV results is less pronounced than the monotonous loadings. However, in the stationary regime the numerical results are in good agreement with exact results. Note that the amplitude of the oscillations in the stationary regime provides the effective shear modulus. Indeed, the macroscopic sinusoidal fields can be considered as complex fields Σ(t) = σ m + σ 0 sin(ωt + δ) and E(t) = ε m + ε 0 sin(ωt), of which the amplitudes are related through the effective constitutive law :

E(t) = (
σ 0 = 2 μ(ω)ε 0 (4.16)
and the effective storage shear modulus is given by : μ(ω) = | μ(iω)| (4.17 It is also worth to look at Figs. 4.8 and 4.9. The first one shows a close-up of the average strain in the elastic fibers around the time they first vanish. The two first chosen instants t = 2.77s and t = 2.91s are represented by the dashed line. Fig. 4.9 represents the strain fluctuations. The dashed lines indicate the two instants abovementioned. These instants coincide with the points where C (m) = C (in) . As reported above, Idiart et al. [START_REF] Martín | Model reduction by meanfield homogenization in viscoelastic composites. II. Application to rigidly reinforced solids[END_REF] discussed the domain where the fluctuations do not decrease. These authors define this condition as the domain of convexity of the free-energy density. According to Fig. 4.9, the free-energy density would appear to periodically lose its convexity. However, in link with their remarks, the EIV method provides accurate results regardless of the violation of the convexity conditions.

C (f) ( ) C (m) ( ) C (in) ( k )

Harmonic response of a fractional viscoelastic matrix reinforced with rigid inclusions

Considering a two-phase composite with a fractional consituent subjected to a harmonic loading of pulsation ω, Gallican and Brenner [START_REF] Gallican | Homogenization estimates for the effective response of fractional viscoelastic particulate composites[END_REF] obtain the effective storage shear modulus for a wide range of pulsations. We refer to their results for our application with the EIV method.

Heterogeneous problem

Let us consider elastic spherical inclusions embedded in a viscoelastic matrix that follows a fractional Zener behaviour. The phases are isotropic and compressible. The bulk modulus of the matrix is assumed to be purely elastic, so that all viscoelastic effects arise from the shear modulus. Denoting the inclusions as phase r = 1 and the matrix as phase r = 2, the fractional Zener constitutive law is specified :

109 σ(t) + (τ (2) σ ) γ D γ σ(t) = µ (2) er : ε(t) + µ (2) eg (τ (2) σ ) γ D γ ε(t) (4.18)
The corresponding relaxation function is expressed as :

µ (2) (t) = µ (2) er + (µ (2) eg -µ (2) er )E γ -(t/τ (2) σ ) γ (4.19)
and in the spectral domain :

μ(2) (ω) = µ (2) er + (µ (2) eg -µ (2) er ) (iωτ σ ) γ 1 + (iωτ σ ) γ (4.20)
The material properties are given by : κ (1) = 40 GPa, and µ (1) = 30 GPa κ (2) = 4 GPa, µ

(2) According to Gusev (2016), the elastic moduli correspond to silica particles and an epoxy resin matrix. For the FFT-based computations, Gallican and Brenner considered several realizations of random distribution of size polydisperse particles. These reference results were confronted to their estimate that involved the resolution of an elastic homogenization problem. They evaluated the Mori-Tanaka, the self-consistent and the generalized selfconsistent (GSC) schemes, and they obtained an excellent agreement with the latter one. We thus focus on the GSC estimate for the EIV method. In elasticity, as noted by Hervé and Zaoui [START_REF] Hervé | Modelling the effective behavior of nonlinear matrixinclusion composites[END_REF], the effective bulk modulus provided by the GSC scheme coincides to the lower Hashin-Shtrikman bound, meanwhile the effective elastic shear modulus is the solution of the second-order equation :

er = 0.01 GPa, µ (2) eg = 1 GPa, τ (2
A μ µ (2) 2 + 2B μ µ (2) + C = 0 (4.22)
where the coefficients A, B and C depend on the other mechanical properties [126] [125].

Finally, the material is subjected to the following harmonic loading of pulsation ω :

E(t) = b 1 + b 2 sin(ωt) (e 1 ⊗ e 1 -e 2 ⊗ e 2 ) (4.23)
In [START_REF] Gallican | Homogenization estimates for the effective response of fractional viscoelastic particulate composites[END_REF], Gallican and Brenner consider the range of pulsation ω ∈ [10 -2 ; 10 4 ] rad.s -1 and the different volume fractions c 1 = 0.1, 0.3 and 0.5. They assess the effective real part of the shear modulus with respect to the frequency. In the next section, we specify how to deal with this problem by means of the EIV method.

Discretization of the fractional Zener's relaxation function

In order to incorporate the fractional viscoelastic Zener constituent into the EIV algorithm, we employ the improved midpoint-based discretization procedure presented in section 3.

K (t) = µ (2) er +(µ (2) eg -µ (2) er ) sin(γπ) πτ γ e -t/e Λ e Λ e Λ 0 z γ P γ (z/τ ) dz+ Π -Λ K -1 K-1 k=1 e γz k P γ (e z k /τ ) e -t/e z k + τ 2γ γ C -B 2 /4 π 2 -arctan e γΠ + B/2 C -B 2 /4 , (4.24)
and in the spectral domain :

μ(2) K (ω) = µ (2) er +(µ (2) eg -µ (2) er ) sin(γπ) πτ γ e -Λ iω iω + e Λ e Λ 0 z γ P γ (z/τ ) dz+ Π -Λ K -1 K-1 k=1 e γz k P γ (e x k /τ ) iω iω + e z k + τ 2γ γ C -B 2 /4 π 2 -arctan e γΠ + B/2 C -B 2 /4 (4.25)
where

z k = Λ + (k -1/2) (Π -Λ) (K -1)
and with B = 2τ γ cos(γπ) and C = τ 2γ .

We empirically choose the finite endpoints of the dicrete subdomain λ and m such that they provide the most accurate function for the lowest number of branches. We determine by trial and error λ = 1.53 × 10 -3 s and m = 3.16 × 10 6 s, and the number of branches K = 11. These parameters provide the two extreme values of shear moduli µ The representation in the spectral domain allows to highlight the specific two asymptotic elastic regimes of a Zener constituent. At short times, the relative error between the real glassy shear modulus µ [START_REF] Tressou | Contribution à l'homogénéisation des milieux viscoélastiques et introduction du couplage avec la température par extensions d'une approche incrémentale directe[END_REF] eg and the approximation is 1.8 × 10 -3 %. At long times, the value of the relative error for µ [START_REF] Tressou | Contribution à l'homogénéisation des milieux viscoélastiques et introduction du couplage avec la température par extensions d'une approche incrémentale directe[END_REF] er is 8 × 10 -3 %.

Homogenization estimates with the EIV method

We consider the sinusoidal strain loading given in Eq. (4.23) with b 1 = 0.5 and b 2 = 5 × 10 -3 . Contrary to homogenization methods involving the correspondence principle, the incremental variational approach provides the mechanical response of the fields in the direct domain. As a result, this method gives the thorough evolution of the fields over time, specifically the short transient response ocurring before the stationary regime. In the latter regime, the material displays a delayed harmonic effective stress response with the same pulsation ω as the loading. The effective modulus is given by the amplitude σ 0 of the oscillations such that : According to the discretization of the Zener's relaxation function, the application of the EIV method consists in the resolution of K = 11 systems of 2 non-linear equations. We take ∆t = 0.005 s on the path t ∈ [0, T ] with T = 20s, and the auxiliary thermoelastic problem is solved with the GSC scheme. Let us illustrate the results for c 1 = 0.1 and ω = 15 rad/s. The macroscopic and local stress responses are shown in Fig. 4.11, where the subscript . (m) refers to the matrix, and . (i) to the spherical inclusions. The graph focuses on the transient and the stationary regimes. In the latter, the amplitude of the effective sinusoidal stress response provides the effective storage shear modulus.

μ(ω) = σ 0 2b 2 (4.26) 4
We also look at the fluctuations of the total and viscous strains in Fig. 4.12. The fluctuations of the strains in the spherical particles are negligible compared to those of the matrix, although they are not null. The 11 plain lines refer to the fluctuations associated to each internal variable. The fluctuations of the viscous strain associated to the smallest characteristic time C (m) (α k=1 ) (z k=0 ∼ 10 -3 s) and the strain fluctuations in the matrix amost coincide. On the contrary, the fluctuations associated to higher values of characteristic times are significantly delayed in time. sampling, leading to high computational costs. This constraint does not appear for a low value of ω, but in this case the stationary regime would occur in a long time, requiring a large time vector. This would again cause high computational costs. For these reasons, we consider the reduced range of pulsations : ω = 5, 10, 15, 20 and 25 rad.s -1 . The EIV method is applied for the volume fraction c 1 = 0.1 and compared to the FFT results in Fig. 4.13.

The effective shear moduli obtained by the EIV simulation are overestimated in comparison with the expected properties. Actually, we observe a constant ratio μEIV /μ FFT = 2.7

113

between the two graphs. Although the effective shear moduli are overestimated, their dependance with respect to ω is well predicted.

We also ran EIV simulations for higher volume fractions c 1 = 0.3 and 0.5. However, the oscillating response of the stress was not regular and we observed spurious numerical jumps. It is not clear why higher volume fractions of inclusions lead to numerical issues.

Effective shear creep function of polycrystalline ice with a local fractional viscoelastic slip system

In this section, we first present the microstructure of polycrystalline ice and its effective behaviour in elasticity given by Schapery [START_REF] Schapery | Linear Elastic and Viscoelastic Deformation Behaviour of Ice[END_REF]. We then confront the EIV approach to the experimental results obtained by Sinha. The auxiliary thermoelastic problem is appropriately solved with the self-consistent estimate, and we consider the local viscoelastic slip system as a fractional Maxwell law.

From the elastic to the viscoelastic deformation analysis of sea ice

Microstructure of S2-ice

We consider S2-ice, which is a common form of sea ice. The single crystal grains have a hexagonal structure seen in Fig. 4.14. Attention will be brought on the two slip systems along the prismatic and the basal planes specified on the illustration. In S2-ice, the grains system grow in a direction, say z, that is perpendicular to the c-axes, as illustrated in the left of

Elastic deformation of S2-ice

The grains are transversely isotropic as well. The moduli are temperature-dependant, and their values are obtained by Gammon [START_REF] Gammon | Elastic Constants of Artificial and Natural Ice Samples by Brillouin Spectroscopy[END_REF] for the fixed temperature T = -20 • C. We consider their local coordinate system e 1 -e 2 -e 3 with the c-axis oriented along the e 3 direction. Note that the elastic slip system thus corresponds to the local modulus L 2323 , or L 44 in Voigt notation. The five independant elastic moduli are given in GPa for temperature T = -20 All moduli are assumed to have the same dependance to temperature, such that they all follow the law :

L ij (T ) = L ij (T r ) 1 -qT 1 -qT r , ( 4.28)
where q = 1.418 × 10 -3 ( • C) -1 , T is in • C and T r is the reference temperature.

We consider a uniaxial compressive stress loading Σ y . The effective Young modulus in the y direction is characterized by an in-plane shear and a compressibility contributions, such that the effective strain is related to the stress through the relation :

E y = 1 4μ ip + 1 4k σ Σ y (4.29)
where the subscript . ip denotes the in-plane shear modulus. Assuming that L 11 ≈ L 44 accordingly to (4.27), the effective in-plane shear modulus can be given by :

μip = k D 2 -1 (4.30)
where

D 2 = (k + L 44 )(k + L ′ 66 ) L 44 L ′ 66 , L ′ 66 = L 11 -L 13 2 (4.31)
and the plane stress bulk modulus is given in [START_REF] Schapery | Linear Elastic and Viscoelastic Deformation Behaviour of Ice[END_REF].

Furthermore, the expression of out-of-plane shear modulus was determined by Benveniste [START_REF] Benveniste | Exact connections between polycrystal and crystal properties in two-dimensional polycrystalline aggregates[END_REF] in this very same case of transversely isotropic polycrystal. The following expression is exact, regardless of the value of L 44 : μ = L 44 L 66 (4.32)

Viscoelastic deformation analysis

According to [START_REF] Weertman | CREEP DEFORMATION OF ICE[END_REF], the creep effects induced by shearing are significantly due to local slip along the basal planes. This effect can be emphasized by the presence of brine pockets on basal planes. That is why Schapery made the assumption that all viscoelastic effects arose from the basal slip system, ie. to L 44 . All other moduli remain elastic. The compliance viscoelastic modulus noted D 44 is assumed to follow the power law :

M 44 (t) = 1 L eg 1 + t τ γ (4.33) system
where L eg is the elastic modulus occuring instantaneously. Its value is given by L 44 found in (4.27). This choice of retardance function actually coincides to the fractional Maxwell law.

The estimate developed by Schapery relies on the correspondence principle, so the effective shear modulus is given in the Laplace domain by the expression (4.32). This author thus obtained an expression of the creep compliance that naturally displays a behaviour ∝ t γ/2 at long times. Furthermore, let us recall that the creep function of S2-ice is expected to follow a power law of order 1/3 [START_REF] Da | On The Viscous Flow in Metals, and Allied Phenomena[END_REF]. It can thus be deducted that the fractional order of the grain should be γ = 2/3, as supported by the accuracy of Schapery's estimate with the experimental results. However, Schapery didn't provide any numerical value for the characteristic time. It is only specified that the choice of τ rigidly shifts the macroscopic creep function along logarithmic times.

Results on the creep response of S2-ice with the EIV method

Stress loading control with the EIV method

The treatment of stress control loading into the EIV method requires an additional step into the algorithm 2. Following [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] in the context of the FFT-based algorithm, this can be done iteratively. Starting with an initial guess fot the effective strain, the accordance with the real effective stress is evaluated iteratively. This additional step implies that the K systems of 2 non-linear equations are solved several times at each time step. The additional numerical bloc is presented individually in Algorithm 3.

Algorithm 3 Additional instructions in the EIV method for stress loading control

Require: E i , C 0 , ϵ while C -1 0 : (Σ -S)

||E i+1 || > ϵ do
Calculate the effective stress from the EIV method : S = L : E i + τ E i+1 = C -1 0 : (Σ -S) + E i end while

Comparison with the analytical case of classical viscoelasticity

As mentionned earlier, S2-ice is isotropically transverse. This symmetry is due to the microstructure of S2-ice, which is characterized by the isotropically random distribution of grains orientations. The self-consistent estimate for the anti-plane effective shear modulus given by Eq. (4.32) is exact for any isotropically transverse microstructure. We thus simplify the problem into a polycrystal with only two orientations in equal proportions, with rotation π/2 rad between each one. The plane isotropy is kept by considering a checkerboard arrangement of the grains, as illustrated in Fig. 4.16. According to [START_REF] Suquet | Numerical Simulation of the Effective Elastic Properties of a Class of Cell Materials[END_REF], this constitutes a special case where the expressions of in-plane and the out-of-plane shear moduli coincide. and the inverse transform yields the exact macroscopic creep function : where J 0 and J 1 are the Bessel functions of the first kind, of order 0 and 1 respectively. Regarding the EIV simulation, we consider the shear stress loading :

F (t
Σ(t) = 1 2
Σ(e x ⊗ e z -e z ⊗ e x ) (4.37)

The effective creep function resulting from the EIV method is compared to the analytical function in Fig. 4.17. The moduli L 44 and L 66 are taken at temperature -20 • C, and we take τ = 2.65 × 10 -2 s. The agreement is excellent and allows us to dive into the case of a fractional viscoelastic slip system.

EIV results and comparison with experimental data

The compressive creep loading experiments were carried out by Sinha with the constant load Σ y = 0.49 MPa, and for different temperatures. He took advantage of the William-Landel-Ferry lawor time-temperature equivalence principlethat allows to shift the system different results with respect to a reference temperature. In doing so and taking T = -10 • C, he could obtain a unique graph representing the effective strain response over almost 300 s. The experimental data E exp y take into account the compressible contribution, following Eq. (4.29). As the compressibility is purely elastic, we will consider a purely deviatoric problem, so the numerical results will be compared to :

F exp (t) = 1 2μ = 2 E exp y Σ y - 1 4k σ (4.38)
For the numerical EIV simulation, we finally consider that the basal slip system follows the fractional relaxation Maxwell law :

L 44 (t) = L eg E γ - t τ γ (4.39)
with the parameter γ = 2/3. As mentioned earlier, the characteristic time is not specified. We thus take τ = 0.05, 0.15, 1 and 10 s. The application of the improved midpoint-based discretization procedure is applied to each case with the same number of branches K = 12. The numerical values for the moduli are taken for temperature T = -10 • C. We show the result of this approximation for τ = 0.05 s in Fig. 4.18.

Because of the large number of iterations for the resolution of the 12 systems of 2 non-linear equations required for the EIV estimate, we consider the creep loading over 4 s only. The graphs 4.19 show the macroscopic creep function for different characteristic times obtained with the EIV method, in lin-lin and semi-log scales. The EIV results are compared with the experimental data of which the elastic compressible contribution is removed, following Eq. (4.38). At t = 0, since the instantaneous response of Maxwell constituents is purely elastic, the polycrystal is supposed to display a purely elastic response as well. The creep function of the EIV estimates have the expected exact value F (t = 0) = 1/2 L eg L 66 ≈ 0.16, meanwhile for the experimental data, the value is ≈ 0.14. This low inaccuracy may be due either to an imprecise modelisation of the local laws, or to the experimental conditions. Nevertheless, we observe a correct agreement for the characteristic time τ = 10s. The semi-log graph allows us to make two main remarks :

(1) The change of the characteristic time results in the rigid shift of the creep function along the logarithmic time vector, as well as Schapery's estimate. (1) The creep function is expected to follow a power law of order 1/3 at long times following Andrade's law. For the material parameter τ = 0.05 s, this behaviour can be observed at t = 2s = 40τ .

Conclusion

In this chapter, we obtained homogenization estimates for two composite materials with fractional viscoelastic constituents, by means of the EIV method. The first case consisted in the evaluation of the effective shear modulus of a reinforced Zener matrix through the results of a harmonic loading. The EIV method does not technically allow to consider a range of pulsations as wide as methods relying on the correspondence principle. Nevertheless, we could obtain estimates on a restrained range that followed the expected tendency with respect to ω. However, they are overestimated compared to reference results. The EIV method can be expected to display pathological results for harmonic loadings. Even though the EIV can give accurate results for a unique internal variable, it may not be the case for several internal variables.

We then use the EIV method with control stress loading. Numerically, this consists in an additional iterative block in the global EIV algorithm, resulting in a higher number of resolutions during a simulation. For the 12 internal variables we consider, this importantly impacts the computational cost. However, we could obtain interesting results by comparison with experimental data. Even though the characteristic time was not explicited in the reference article, we could investigate the expected tendancie of the effective properties. More specifically, the EIV accurately predicts the physical power-law behaviour at long times.

Overall conclusion and perspective

Conclusion

This PhD thesis dealt with the homogenization of composite materials with linear fractional viscoelastic constituents. We first gave an overview of homogenization methods in the context of linear elasticity. In this context, we presented our contribution to numerical homogenization methods, specifically FFT-based methods. More precisely, we addressed the inherent numerical issue due to the Gibbs phenomenon that causes spurious oscillations on the interfaces of the microstructure. We took an interest on the smoothing data filter developed by Morin et al. [START_REF] Morin | Periodic smoothing splines for FFT-based solvers[END_REF] on a two-phase material with a checkerboard microstructure. This case is particularly challenging because of the perioding singular points joining the phases. Indeed, the data filter can not significantly improve the accuracy of the local fields. Therefore, we investigated a combination of the filter by introducing a weighting parameter that balances the contribution of the rigid smooth tensor, and the inverse smooth compliance tensor. Regarding the effects on the effective stiffness tensor, any value of the weighted parameter h < 1 successfully leads to the nullity of off-diagonal terms. Also, the value of the effective shear modulus is impacted; however there is no clear criterion for the choice oh h. On the contrary, the weighting parameter does not significantly affect the local fields, and no improvement was found in comparison to the exact solution.

We then investigated the homogenization of viscoelastic materials by means of the Effective Internal Variable method developed by Lahellec and Suquet [START_REF] Lahellec | Effective behavior of linear viscoelastic composites: A time-integration approach[END_REF]. The incorporation of linear fractional viscoelastic materials into this method relies on the formulation of the constituents constitutive laws in the Generalized Standard Materials framework. The rheological interpretation of the fractional damping unit as a superposition of Maxwell elements allows to establish the existence of free-energy and dissipation potentials. They are indeed convex functions of their variables, with an infinite number of internal variables. The latters are associated to characteristic times and their distribution is given by the explicit expression of the spectrum, as a power-law function. Therefore, the choice of the finite number of internal variables, as well as the characteristic times and moduli, can be appropriately made by the discretization of the spectrum. We made use of the discretization procedure developed by Papoulia et al. [START_REF] Papoulia | Rheological representation of fractional order viscoelastic material models[END_REF]. It is an improved midpoint-based method that provides an accurate approximate function with a low number of parameters. We applied this procedure to the Mittag-Leffler function that is involved in the relaxation function of three (among the four) types of fractional viscoelasicity [START_REF] Caputo | Linear Models of Dissipation whose Q is almost Frequency Independent-II[END_REF]. The spectra are thus rational functions of polynomials of non-integer order. The discretization procedure required assumptions that justified a posteriori.

Finally, we presented the EIV method for several internal variables, as developed by Tressou et al. [START_REF] Tressou | Application of the variational EIV approach to linear viscoelastic phases governed by several internal variables -Examples with the generalized Maxwell law[END_REF]. We gave insights on our numerical implementation of the EIV algorithm. We used this incremental variational approach on two different problems. We first considered spherical rigid inclusions embedded in a fractional Zenerian matrix, subjected to harmonic loadings. The auxiliary thermoelastic problem was solved by means of the Generalized Self-Consistent estimate. We could recover the expected tendency of the effective shear modulus with respect to the pulsation. However, the obtained estimates are not in agreement with full-field reference results [START_REF] Gallican | Homogenization estimates for polymer-based viscoelastic composite materials[END_REF], and we encountered numerical issues depending on the volume fraction of the inclusions. Then, a polycrystal subjected to a shear creep loading was considered. The stress control loading required to add instructions into the EIV method, following [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. Solving the auxiliary thermoelastic problem with the self-consistent estimate, we obtained effective creep functions in accordance with physical considerations, namely the power law behaviour at long times, and the invariance along the logarithmic scale with respect to the characteristic time.

Perspectives

In this work, we developed the appropriate strategy to tackle the homogenization of composites with fractional viscoelastic constituents by means of the Effective Internal Variable approach. However, the results provided by the numerical application were not systematically satisfying, in particular for harmonic loadings. The following points could be investigated in order to both continue the evaluation of the method for such problems, and improve the numerical treatment of the large systems of non-linear equations :

• A short-term perspective consists in comparing EIV estimates for monotonous loadings with full-field methods as reference results. The local constitutive laws can be incorporated following the discretization procedure we used in this PhD thesis. For a start, the hypothesis of purely deviatoric internal variables could be kept.

• The previous point should be deepened following Tressou et al. [START_REF] Tressou | Application of the incremental variational approach (EIV model) to the linear viscoelastic homogenization of different types of microstructures: long fiber-, particle-reinforced and strand-based composites[END_REF], by evaluating the accuracy of EIV method with fractional constituents for different microstructures or in other words, for different homogenization schemes. This would allow to establish the accuracy of the method for the large systems of non-linear equations that are considered. As the discretization of the fractional relaxation functions leads to moduli and characteristic times with high contrasts, this verification work would be of great interest.

• As mentioned in Chapter 3, numerical issues can be encountered in the resolution of the non-linear system of equations, either due to the square root implied in one

B

Computation of the Mittag-Leffler function

The Mittag-Leffler function is defined by the power series : In the context of viscoelasticity, it reduces to 0 < γ ≤ 1 and β = 1 :

E γ (z) = +∞ n=0 z n Γ(1 + γn) , ∀z ∈ C, 0 < γ ≤ 1, (B.2)
Alternatively, this function can be represented in its integral form : The computation of this function requires the use of different numerical techniques in distinct parts of the complex plane. In the case of viscoelasticty, z is a purely real variable, and 0 < γ ≤ 1. The algorithm developed in [START_REF] Gorenflo | Computation of the Mittag-Leffler function[END_REF] reduces to 4 : Obnosov developed explicit solutions for periodic heterogeneous structures submitted to an imposed conductivity field [START_REF] Yurii | Periodic Heterogeneous Structures: New Explicit Solutions and Effective Characteristics of Refraction of an Imposed Field[END_REF]. This solution is extended to elasticity in [START_REF] Bellis | Eigendecomposition-based convergence analysis of the Neumann series for laminated composites and discretization error estimation[END_REF]. with the following definitions :

E γ (-t γ ) =
∆ = µ -1 1 -µ -1 2 µ -1 1 + µ -1 2 (C.3)

. 2 )Figure 1 . 1 :

 211 Figure 1.1 : Creep test. Stress step-function (left) and strain response (right).

. 5 )Figure 1 . 2 :

 512 Figure 1.2 : Relaxation test. Strain step-function (left) and stress response (right).

  ) : σ(p) = L(p) : ε(p) (1.25) providing the expression of the relaxation function L(t) : L(p) = P (p) -1 : Q(p). (1.26) Dually, the relation between the relaxation function and the creep function is provided in the Laplace-Carson domain with : L(p) : M(p) = M(p) : L(p) = I (1.27)

Figure 1 . 4 :

 14 Figure 1.4 : Generalized Maxwell.Figure 1.5 : Generalized Kelvin-Voigt.

Figure 1 . 5 :

 15 Figure 1.4 : Generalized Maxwell.Figure 1.5 : Generalized Kelvin-Voigt.

Figure 1 . 6 :

 16 Figure 1.6 : Representation of the fractional dashpot or spring-pot element

1 . 7 .Figure 1 . 7 :Figure 1 . 8 :

 171718 Figure 1.7 : Creep function of the Maxwell constituent.

Figure 1 . 9 :Figure 1 . 10 :

 19110 Figure 1.9 : Creep function of the Kelvin-Voigt constituent.

Figure 1 . 11 :Figure 1 . 12 :

 111112 Figure 1.11 : Creep function of the Zener constituent.

Figure 1 . 13 :Figure 1 . 14 :

 113114 Figure 1.13 : Creep function of the anti-Zener constituent.

Figure 2 . 1 :

 21 Figure 2.1 : Checkerboard-like microstructure

Figure 2 . 2 :

 22 Figure 2.2 : Shear modulus µ(x) for the checkerboard-like modulus with contrast c = 100.

Fig. 2 .Figure 2 . 3 :Figure 2 . 4 : 1 ⊗ e 3 + e 3 ⊗ e 1 ) 2 (e 2 ⊗ e 3 + e 3 ⊗ e 2

 22324133122332 Fig. 2.3 represents the shear modulus along x at y = L 2 /4, showing that the smooth parameter controls the thickness of the interface. Fig. 2.4 shows the component stress σ 13 resulting from the different computations, where the smoothing effect of s is clearly seen. They are compared to the exact field in the top right corner of the squared inclusion.

  (a) Exact. (b) s = 10 -3 . (c) s = 5.10 -3 .

Figure 2 . 6 :

 26 Figure 2.6 : Shear modulus µ(x) for the checkerboard modulus with contrast c = 100.

Figure 2 . 7 :

 27 Figure 2.7 : Comparison of ε 13 (x) along the dotted line for the checkerboard microstructure and for c = 100. The figure on the right is a close-up on the oscillations.

Figure 2 . 8 :

 28 Figure 2.8 : Comparison of the field ε 13 (x) for c = 2 (left) and close-up on the center of the microstructure (right). The first line shows the exact field, the second line shows the computation with s = 0 and the third line for s = 5.10 -3 .

Figure 2 . 9 :

 29 Figure 2.9 : Comparison of the field ε 13 (x) for c = 100. As the previous figure, the first line shows the exact field, the second line shows the computation with s = 0 and the third line for s = 5.10 -3 ; and the right column is a close-up of the fields in the center of the microstructure.

8 Figure 2 . 10 : 8 Figure 2 .

 821082 Figure 2.10 : Effect of the smoothing on the checkerboard-like microstructure.

( 1 )

 1 No oscillations on the local stress and strain fields, (2) Accuracy of the effective stiffness tensor with the exact solution : on the effective shear modulus and on the off diagonal terms The effect of different values of h on the microstructure is shown on Figs. 2.12a, 2.12b and 2.12c; and the detail on the interface is shown in Fig. 2.13. (a) Exact. (b) s, h = 0.9. (c) s, h = 0.1.

Figure 2 . 12 :Figure 2 . 13 :

 212213 Figure 2.12 : Effect of the parameter h for s = 5 × 10 -3 on the shear modulus µ 2323 (x) with contrast c = 100.

Figure 2 . 14 :

 214 Figure 2.14 : Effective shear modulus depending on h, and s = 5 × 10 -3 fixed.

1 Figure 2 . 15 :

 1215 Figure 2.15 : Exact ε 13 (x) along the dotted line, compared to the FFT-based simulations with and without smoothing, for c = 2. The figure on the right is a close-up on the singular point. The smoothing parameter is fixed : s = 5 × 10 -3 .

Figure 2 . 16 : 3 . 2

 21632 Figure 2.16 : Comparison between the exact and numerical ε 13 (x) along the dotted line, for c = 100. The figure on the right is a close-up on the softer phase. The smoothing parameter is fixed : s = 5 × 10 -3 .

. 3 )

 3 and the effective constitutive law : σ(p) = L(p) : ε(p), where L = L(p) : Â(p)(3.4)

Figure 3 . 1 :

 31 Figure 3.1 : Rheological representation of the fractional dashpot's relaxation function.

α 1 ,

 1 α 2 ,... ⟨E 0 ⟩ is referred to as the condensed incremental potential.The local auxiliary properties L

  following effective energy, only expressed by means of the local properties, where ∆L = L

Figure 3 . 2 :

 32 Figure 3.2 : Comparison between the collocation and the improved midpoint-based approximations

∆

  may be considered together with all the other unknown internal variables α (r) k,n . Indeed, in this case, Eq. (4.3) is written as :⟨α k ⟩ (r) = B (r) k,1 : A (r) E + B , ∀k = 1, ..., K(4.4)(b) Inherent limitations of the method

Δ

  , ⟨ε⟩ (r) , ⟨ε ⊗ ε⟩ (r) Method : systems of the 2 non-linear equations Output : optimized values of ⟨ε⟩ (r) , ⟨ε ⊗ ε⟩ (r)

Figure 4 . 1 : 103 4. 1 . 2

 4110312 Figure 4.1 : Illustration of the structure of the numerical implementation in Python

Figure 4 . 2 :Figure 4 . 3 :

 4243 Figure 4.2 : Constant strain loading, local mean strains and inelastic mean strain.

Figure 4 . 4 :

 44 Figure 4.4 : Strain fluctuations in the matrix and in the inclusions for a relaxation loading.The inelastic fluctuations are also shown.

Fig. 4 .

 4 Fig.4.2 shows the macroscopic strain loading where we fixed E 11 = 1, as well as the local strains in the matrix. The fibers are denoted by the superscript . f and the matrix by . m . This graph also shows the viscous mean strain. As expected, at long times, the viscous means strain and the total strain in the matrix coincide. Note that, as the branches are identical, all the inelastic strains are equal. The macroscopic and local stress responses are shown in Fig.4.3 and are in excellent agreement with the analytical solution. Fig.4.4 shows the fluctuations of the strain fields in the matrix and in the fibers. Let us recall that the fluctuations are expressed by :

e 1 -

 1 e 2 ⊗ e 2 ), with b = 5.75 × 10 -2 (4.14) Here again, Figs. 4.5 show the imposed strain and the local averaged strains. The global and local stresses are satisfactorily compared with the analytical solution on the right figure.

Figure 4 . 5 :

 45 Figure 4.5 : Constant strain rate loading and mean strains (left). Local and global stress response and comparison with the analytical solution (right).

  b 1 + b 2 sin(ωt)) (e 1 ⊗ e 1 -e 2 ⊗ e 2 ), with b 1 = 0.5, b 2 = 5 × 10 -3 and ω = 10.5 rad/s (4.15)

  )

Figure 4 . 6 :

 46 Figure 4.6 : Sinusoidal strain loading and resulting local strains (left). Global and local stress responses with comparison with the exact solution (right).

Figure 4 . 7 :

 47 Figure 4.7 : Close-up on the stress response for a harmonic loading. The vertical dashed lines report the first two instants where all stress components vanish, at t = 2.77s and t = 2.91s.

Figure 4 . 8 :

 48 Figure 4.8 : Close-up of the averaged strain field in the elastic fibers around when it first vanishes. These instants are represented by the dashed lines at t = 2.77s and t = 2.91s

Figure 4 . 9 :

 49 Figure 4.9 : Fluctuations of the strain fields in the matrix and in the inclusions, as well as the inelastic fluctuations. The dashed lines coincide with C (m) = C (in)

  ) σ = 0.0265 s and γ = 0.7(4.21)

3 . 2 .

 32 The discretization of the Mittag-Leffler function given by Eq. (3.100) and (3.101) are integrated in Eq. (4.[START_REF] Eldred | Kelvin-Voigt versus fractional derivative model as constitutive relations for viscoelastic materials[END_REF]) and (4.20) respectively, such that :

  µ

( 2 )

 2 k=1 = 2.3 × 10 -2 GPa, and µ (2) k=11 = 2.6 × 10 -6 GPa. Given the material properties of the matrix (4.21), we present the approximate functions in Figs. 4.10.

Figure 4 . 10 :

 410 Figure 4.10 : Exact relaxation function of the Zener constituent compared to the discretized function obtained from the improved midpoind-based approximation in the direct domain (left) and in the spectral domain (right). We take K = 11 branches, λ = 1.53 × 10 -3 s and m = 3.16 × 10 6 s. The dotted lines represent the moduli µ eg (red)and µ er (green).

Figure 4 . 11 :Figure 4 . 12 :

 411412 Figure 4.11 : Macroscopic and local stress responses to a harmonic loading. The transient regime and the stationary regime can be distinguished

Figure 4 . 13 :

 413 Figure 4.13 : Effective storage shear modulus of a matrix-inclusions viscoelastic composite for c 1 = 0.1. Comparison between the FFT-based computations (GB) from [93] and the EIV simulation.

Fig. 4 .

 4 15. Therefore, the basal planes of the single grains are in vertical planes. The c-axes are randomly oriented in the horizontal x-y plane, as indicated by the arrows in the right of Fig. 4.15. The distribution of the orientations is isotropically random. The polycrystal is thus characterized by a transversely isotropic symmetry.

Figure 4 . 14 :

 414 Figure 4.14 : Hexagonal structure of a grain ice.

Figure 4 . 15 :

 415 Figure 4.15 : Illustration of S2-ice from Schapery [17] (left). Optimal micrograph of the cross section of the polycrystal in the x-y plane obtained by Sinha [123] (right). The arrows indicate the c-axes of the grains.

Figure 4 . 16 :

 416 Figure 4.16 : Illustration of the checkerboard polycrystal. The arrows indicate the c-axes of the grains.

Figure 4 . 17 :

 417 Figure 4.17 : Creep function in semi-log scale of S2-ice with a classical viscoelastic basal slip behaviour. Comparison between the EIV results and the analytical function.

Figure 4 . 18 :

 418 Figure 4.18 : Exact fractional Maxwell relaxation function compared to the improved midpoind-based approximation in the direct domain (left) and in the spectral domain (right), for τ = 0.05 s and γ = 2/3. We take K = 12 branches, λ = 1.53 × 10 -3 s and m = 3.16 × 10 6 s. The dotted line represent the elastic glassy moduli L eg .

Figure 4 . 19 :

 419 Figure 4.19 : Creep function of S2-ice with a fractional viscoelastic basal slip, in linear scale (left) and semi-log scale (right). Comparison between the EIV results for different characteristic times and the experimental data.

E

  γ,β (z) = +∞ n=0 z n Γ(β + γn) , ∀z ∈ C, β ∈ C, γ > 0, (B.1)

+∞ 0 H

 0 γ (ξ)e -t/ξ dξ, t > 0 with 0< γ ≤ 1 (B.3) with H γ (ξ) = 1 π ξ γ-1 sin(γπ) 1 + 2ξ γ cos(γπ) + ξ 2γ (B.4)

142 Algorithm 4

 4 Computation of the Mittag-Leffler functionRequire: 0 < γ ≤ 1, z ∈ R, ρ > 0 the machine precision. if z = 0 then E γ (z) = 1 else if |z|< 1 then k 0 = max 0, ln(ρ(1 -|z|)) ln(|z|) E γ (z) = k 0 k=0 z k Γ(1 + γk) else if |z|> 10 + 5γ then k 0 = -ln(ρ) ln(|z|) if |arg(z)|< γπ 4 max{1, 2|z|, (-ln(πρ/6)) γ } K(γ, χ, z) = -1 γπ e -χ 1/γ zsin(γπ) χ 2 -2χzcos(γπ) + z 2 if |arg(z)|> γπ then E γ (z) = χ 0 0 K(γ, χ, z)dχ else if |arg(z)|< γπ then E γ (z) =

Figure C. 1 : 2 (

 12 Figure C.1 : Checkerboard-like microstructure

Table 1 .

 1 1.

Table 1 .

 1

1 : The four types of classical linear viscoelasticity Constituent Representation M eg M vr L er L vg Maxwell 1 E

Table 1 . 2 :

 12 The four types of fractional viscoelasticity

	Constituent	Representation	M eg	M fr	L er	L fg

2

  Homogenization in elasticity and numerical contribution

51 2.2 Contribution to numerical homogenization : smoothing filter for FFT-based methods 2.2.1 Presentation of the periodic FFT-based solver

Table 2 . 1 :

 21 Comparison of the effective shear modulus for c = 2

		Effective modulus (Pa)	ϵ (%)	Off diagonal term (Pa)
	μexact	1.4142	0	0
	μs=0	1.4142	4.4859 × 10 -5	10 -9
	μs=0.005	1.4173	0.2209	10 -12

Table 2 .

 2 2 : Comparison of the effective shear modulus for c = 10

		Effective modulus (Pa) ϵ (%) Off diagonal term (Pa)
	μexact	3.1623	0	0
	μs=0	3.1624	0.0033	10 -5
	μs=0.005	3.3409	5.6472	10 -3

Table 2 .

 2 3 : Comparison of the effective shear modulus for c = 100

		Effective modulus (Pa)	ϵ (%)	Off diagonal term (Pa)
	μexact	10	0	0
	μs=0	10.0059	0.0588	10 -1
	μs=0.005	18.2571	82.5714	10 -1

  • C : L 11 = 13.81, L 12 = 7.02, L 13 = 5.72, L 33 = 14.89, L 44 = 2.99, L 66 = 3.40 (4.27) 4 EIV estimates of composite materials with fractional viscoelastic constituents 115

by Papoulia et al.[START_REF] Papoulia | Rheological representation of fractional order viscoelastic material models[END_REF] and applied it to the Mittag-Leffler function that is involved in the four types of viscoelasticity. The difficulty relied on the calculation of the integral in the portions (-∞, Λ] and [Π, +∞). For the first subdomain, the function can only be assessed numerically, meanwhile on the other one we could calculate its analytical expression. Besides, in order to express the functions on these subdomains, we made assumptions that we were able to justify a posteriori.

Remerciements

Algorithm and numerical implementation of the EIV method

We provide insights of our numerical implementation, performed in Python. For understanding purpose, the numerical implementation is presented in 2 as a fixed-point algorithm. However note that in practice, the non-linear system of equations requires a more sophisticated solver. Precisely, we do not make use of a fixed-point algorithm and rather choose the solver "root" from the Scipy's optimization package, for which we specify the use of the Levenberg-Marquardt algorithm. Still for clarity, we do not get into the precise root method, but provide a brief illustration of the numerical implementation of the implicit objective function in Fig. 4.1.

Regarding the system of non-linear equations, we regularize the expression of the α (r) k,n by multiplying both sides by the denominator θ (r) k , following Badulescu et al. [START_REF] Badulescu | Field statistics in linear viscoelastic composites and polycrystals[END_REF]. However, we do not follow these authors on the reexpression of θ (r) k as a polynomial (see Eq. 4.1). We rather keep the expression of θ (r) k as (3.49) 1 , i.e. with the square root and the minus sign. Indeed, as mentioned above, the polynomial form (4.1) leads to the existence of multiple roots. As we consider K > 1 systems of equations, we strategically choose to circumvent the consideration of multiple solutions.

Algorithm 2 Fixed-point algorithm of the EIV method for comprehension purpose

Require: α (r)

n and θ (r) (initial guesses), ⟨α n ⟩ (r) and ⟨α n :

n and θ (r) are not steady do Define the auxiliary media L ∆ , τ ∆ and f ∆ Apply the linear elastic homogenization scheme to calculate L, τ , A (r) and a (r) Calculate ⟨ε⟩ (r) and ⟨ε : ε⟩ (r) Calculate ⟨α⟩ (r) and ⟨α : α⟩ (r) Calculate α

Confrontation with reference results

We compare the effective shear modulus obtained from the EIV simulations to the reference FFT results provided by Gallican and Brenner [START_REF] Gallican | Homogenization estimates for the effective response of fractional viscoelastic particulate composites[END_REF]. As mentioned above, the authors considered a wide range of pulsation loadings. As the EIV approach is applied in the direct domain, it is not feasible to consider the whole range. Indeed, large values of ω would require an infinitesimal time step in the EIV implementation for an appropriate

Appendix

A Usual Laplace-Carson transform

The Laplace transform of the function f (t) is defined by the integral :

where p is a complex variable. The Laplace-Carson transform f (p) is defined as :

with μ the effective shear modulus and θ a function, respectively defined by :

Γ(z) being the standard gamma function. Finally, the function χ(z) is defined by :

and dn(z|m) being the complex Jacobian elliptic delta function. For the complex argument z = x + iy, it is expressed as :

C.2 Checkerboard microstructure

We consider a checkerboard field, of which the exact solution was developed by [START_REF] Craster | Four-Phase Checkerboard Composites[END_REF]. We consider the following double periodic rectangular component on the complex plane z as seen in

In the specific case of l = h, then m = m 1 = 1 2 and the solution of the problem reads :

, Γ being the classical gamma function. χ(z) reads :