Homogenization estimates for the effective response of fractional viscoelastic particulate composites

This article is devoted to the micromechanical modelling of the time harmonic response of viscoelastic composites made of fractional Zener constituents. By extending previous results in classical viscoelasticity, new exact relations on time integrals of the effective relaxation spectrum are obtained. They are related to the intraphase second moments of the strain field in the asymptotic elastic regimes at low and high frequencies. Based on these relations, the effective relaxation spectrum is approximated by a sum of two weighted Dirac delta functions. An attractive feature of this viscoelastic homogenization model is that it only involves the resolution of two elastic homogenization problems. This model is applied to estimate the response of particle-reinforced two-phase composites. Its relevance is assessed by performing comparisons with FFT full-field simulations for distributions of polydisperse spherical particles.


Introduction
By contrast with elastic or viscous heterogeneous media, the coupling between conservative and dissipative deformation mechanisms in viscoelastic composite materials leads generally to emerging effective properties (i.e. not present at the scale of the individual constituents). In particular, the mixture of elementary viscoelastic constituents (i.e. with a single relaxation time) presents in general at the macroscopic scale a long-memory effect which manifests itself in the effective constitutive relation through an integral kernel [18,52,57]. It can be evidenced by making use of the correspondence principle [26,38] which allows to transform a viscoelastic problem into symbolic elastic ones in the Laplace-Carson domain ("Appendix A"). It thus gives access to the complex overall properties characterizing the response of the composite media subjected to harmonic loadings [27]. The overall integral kernel can be derived in closed form only in specific cases [2,5,39]. In general, it can be approximated by a Prony series [34,53] which turns out to be exact only if the effective spectrum is a sum of weighted Dirac delta functions [2,5]. The well-known collocation method relies on this series development of the effective spectrum [40,49,61].
The viscoelastic behaviour of many monophase viscoelastic materials (single constituents) does not reduce to one of the four elementary classical viscoelastic models [7]. A possible, and widely used, approach is to resort to a generalized Maxwell (Kelvin-Voigt) model with a finite number of relaxation (retardation) times to represent the viscoelastic moduli of single viscoelastic constituents [60]. It does precisely correspond to a Prony series development for the local integral kernel. An alternative approach consists in making use of Communicated by Andreas Öchsner.
V. Gallican · R. Brenner (B) Sorbonne Université, CNRS, Institut Jean le Rond d'Alembert, 75005 Paris, France E-mail: renald.brenner@sorbonne-universite.fr fractional viscoelasticity [36,46]. In this framework, the differential constitutive equation involves non-integer derivatives of the stress and strain. Fractional calculus models have firstly been used as an empirical method to describe the properties of linear viscoelastic materials due to its ability to model long-memory effects. Experimental observations covering a wide range of materials [20,21,45] motivated Scott-Blair [3] to propose fractional derivative models to improve the description of time-dependent responses of materials. It can be also noted that Bagley and Torvik [1] proposed a physical interpretation of fractional viscoelasticity for polymer materials by establishing a link with the Rouse model. As compared with classical viscoelasticity, the study of effective properties of fractional viscoelastic composites, based on the correspondence principle, has received less attention (see, for instance, [12,50,56]).
Exact relations have been recently obtained on the asymptotic behaviour (in the time or frequency domains) of viscoelastic composites made of classical elementary constituents [6,19,58]. They imply conditions which have to be fulfilled by certain time integrals of the effective relaxation (retardation) spectrum. Besides, they have been used to propose "minimal" approximate homogenization viscoelastic models, for mixtures of Maxwell constituents, based on asymptotic decoupled elastic or viscous homogenization problems [6,58].
Following this approach, the present study aims at proposing a micromechanical modelling of the time harmonic response of viscoelastic polymer composite materials whose local constitutive response is described by a fractional Zener model. New asymptotic exact relations on the storage and loss viscoelastic moduli are derived in the context of fractional viscoelasticity by making use of stationary principles for complex viscoelasticity [23]. They lead to exact conditions on time integrals of the overall relaxation spectrum which are related to non-integer derivative and integral of the effective complex viscoelastic moduli. These relations hold for any microstructure. Approximate homogenization estimates over the whole frequency range are then built in a way similar to the case of classical viscoelasticity. By contrast with the common use of the correspondence principle, the proposed approximate viscoelastic model only requires the resolution of the relaxed and glassy elastic heterogeneous problems. An application of this model is presented for the particular case of an isotropic viscoelastic polymer matrix, with negligible bulk relaxation, containing randomly distributed spherical elastic particles. The relevance of the model is assessed by comparison with fast Fourier transform (FFT) full-field computations [43,44] performed on polydisperse microstructures.

Classical viscoelasticity
From the Boltzmann superposition principle, the stress response of a non-ageing linear viscoelastic material to a given derivable strain loading path ε(u), u ∈ [0; t] with additional discontinuities [ε] i at times t i , and initial conditions σ (t = 0) = 0, classically reads with L(t) the viscoelastic stiffness tensor (i.e. relaxation function) or, with a concise notation, where * and denote, respectively, the usual and Stieltjes convolution products. In the sequel, attention will be restricted to viscoelastic materials which present asymptotic elastic regimes at short and long times. The general expression of the relaxation function thus reads L e r denotes the relaxed elastic stiffness tensor, by reference to polymer materials, while G and τ are the relaxation spectrum and time.

Fractional viscoelasticity
Differential constitutive laws with non-integer derivatives (i.e. fractional calculus [13,46,47]) have been proposed to describe the viscoelastic response of various materials such as polymers [1], polycrystalline ice [54] and rocks [64]. By contrast with the usual derivative, the non-integer derivative of a function f at a certain time t depends on the function history on the range ] − ∞, t] ("Appendix B.2"). The non-integer derivative order α is thus sometimes called "memory parameter". Such laws are required to describe creep compliances with a nonlinear time dependence. They are also useful to describe viscoelastic transient regimes, between two elastic asymptotic states, with few parameters as compared to generalized Maxwell (resp. Kelvin-Voigt) models.
The behaviour law of a fractional constituent (fractional dashpot or "spring-pot" [31]) is a linear relation between the stress σ and the fractional derivative of the strain D α ε with D α ≡ d α / dt α and 0 < α < 1 ("Appendix B.4"). Its relaxation function can be described as a sum of exponentials weighted with a powerlaw relaxation spectrum [35,60]. In terms of rheological elements, a fractional dashpot can thus be described by a continuous generalized Maxwell model. Conversely, it can be noted that hierarchical assemblages of simple rheological elements (springs and dashpots) exhibit a fractional constitutive behaviour [29,55].
By extending generalized Maxwell and Kelvin-Voigt models to the fractional case, it can be shown that the constitutive relation is still of the form (2) with a fractional relaxation function L α (t) given by Koeller [31,32] where E α (t) is the Mittag-Leffler function ("Appendix B.3"). Obviously, when α = 1 the relaxation function L α (t) corresponds to the classical relaxation function L(t) (3). In the case of a relaxation spectrum G(τ ) consisting of a single Dirac delta function, relation (4) is the relaxation function of a fractional Zener model (i.e. fractional standard solid) [8]. Lion [35] has proved that this model fulfils Clausius-Duhem inequality (i.e. non-negativity of dissipation). This led him to the conjecture that any fractional constitutive model with relaxation function (4) is thermodynamically admissible.

Local problem for harmonic loadings
From the integral representation of the fractional viscoelastic constitutive relation, the correspondence principle allows to transform the fractional viscoelastic problem into a symbolic elastic one [26,38]. For the particular case of harmonic loadings, the problem is brought from time to spectral domain and enables the determination of storage and loss moduli characterizing the dynamic response of the material. Hereafter, we consider an heterogeneous medium occupying volume Ω which consists of N homogeneous phases with characteristic function χ (s) (x) and volume Ω (s) (s ∈ [0; N ]). Moreover, it is assumed that Ω (s) Ω and that the phases are perfectly bonded. The fractional viscoelastic relaxation function of phase (s) is denoted L (s) α (t). It follows that the pointwise viscoelastic relaxation tensor L α (x, t) reads with χ (s) (x) = 1 if x ∈ Ω (s) and 0 otherwise. The volume averages over Ω and Ω (s) are, respectively, denoted • and • (s) . By definition of the characteristic function, the volume fraction of phase (s) is c s = χ (s) . The response of a fractional viscoelastic heterogeneous media to a sinusoidal loading is classically studied in the spectral domain by considering the Laplace-Carson transform ("Appendix A") of the constitutive equation for a purely imaginary transform variable p = iω (i 2 = −1) [27]. By assuming an overall strain loading ε(t) = ε * e iωt , the local problem corresponding to the steady-state regime at angular frequency ω reads with prescribed boundary conditions. The complex fields (σ * , ε * , L * α ) are the time LC transforms of the fields (σ , ε, L α ).

Local and effective viscoelastic properties
The fractional standard linear solid (fractional Zener model) [8] is defined, in a general three-dimensional context, by the following homogeneous fractional differential equation (7) with I the fourth-order identity tensor. It exhibits asymptotic elastic behaviours with glassy moduli L e g at short times (t → 0) and relaxed moduli L e r at long times (t → +∞). L f represents the constitutive fractional viscous modulus. The viscoelastic relaxation tensor of phase (s), defined by relation (7), reads Also, it is worth noting that the eigenvalues of L (s) f : (G (s) ) −1 correspond to the fractional relaxation times τ α (s) (units: s α ) of the constituent. Following [31], it is assumed that, in the general case, the effective relaxation tensor L(t) can be written in the form As a consequence, the complex viscoelastic relaxation tensor L * α , which characterizes the steady-state harmonic regime at angular frequency ω, is given by The overall complex constitutive relation is and the complex relaxation tensor admits the decomposition with L α , L α the overall storage and loss moduli which are, respectively, proportional to the stored and dissipated energies [60]. By noting that i α = e iπα/2 , they read with q = 1 + 2(ωτ ) α cos πα 2 + (ωτ ) 2α . The effective loss tensor η(ω α ) which characterizes the damping is classically defined by Finally, it is important to note that as ω → +∞ or ω → 0 the local fields are asymptotically solutions of purely elastic heterogeneous problems. The pointwise complex strain field ε * (x, iω) thus satisfies with ε g (x) and ε r (x) the real strain fields which are solutions of the heterogeneous glassy and relaxed elastic problems. The same asymptotic properties hold for the stress field σ * (x, iω) with asymptotic fields σ g (x) and σ r (x), respectively.

Saddle-point variational principles in complex viscoelasticity
Let us introduce the following forms of the complex constitutive relation (6) 1 rewritten as systems of real equations where (σ , ε , L α ) and (σ , ε , L α ) are real fields, that is where the position vector x has been omitted for conciseness. By considering the frequency-dependent complex two saddle-point variational principles on the effective complex "energy" φ * can be established [9,19] Re and The solution fields are functions of ω α . It is stressed out that the functionals (19) and (20) do not have a physical meaning by contrast with the minimal variational principle [9,23,41] which considers the average dissipated energy over a period of oscillation, that is φ d = (ω/2) σ : ε − σ : ε . However, Re( φ * ) and Im( φ * ) are stationary with respect to the strain fields ε and ε . Consequently, use can be made of a lemma on the derivative of the stationary value of a function with respect to a parameter (see "Appendix B". in [48]). Derivatives of (19) and (20) with respect to parameter h thus read

Exact relations on the overall complex storage and loss moduli
We can take advantage of the saddle-point variational principles (19) and (20) and their derivatives with respect to a parameter (21) to obtain asymptotic properties on the overall storage and loss moduli at low and high frequencies. By considering an overall strain such that ε = ε and ε = 0, the stationary principle on the real part of φ * (19) leads to On the other hand, the stationary principle on the imaginary part of φ * (20) implies that Relations (22) and (23) correspond to a well-known result for mixtures of classical Zener constituents stating that the asymptotic overall behaviours at low and high frequencies are purely elastic. This property still holds in the fractional case since the fractional feature of the constitutive behaviour only affects the transient regime.
The determination of the asymptotic properties L e r and L e g thus simply requires to solve the heterogeneous elastic problems, respectively, in the relaxed (ω → 0) and glassy (ω → +∞) states. Besides, by choosing h = ω α , the derivatives of the stationary principles (21) provide the results (24) It can be remarked that unlike classical viscoelasticity [19] the asymptotic values of the derivative of the effective storage modulus L α do not vanish since 0 < α < 1. These asymptotic values of the derivative of the real and imaginary parts of the overall complex "energy" φ * are thus exactly determined with the local properties of the constituents and the intraphase second moment of the strain solution fields of the purely elastic heterogeneous problems at low (relaxed state) and high (glassy state) frequencies.

Consequences on the overall relaxation spectrum and their physical meaning
With the general expression of the overall complex relaxation function L * α (iω) (10), the asymptotic relations The first relation has a simple and well-known physical interpretation in the case of a relaxation loading test (ε(t) = ε 0 ). Indeed, we have The integral of the relaxation spectrum is thus directly related to the overall stress gap between short and long times (see, for instance, [16]). Obviously, it does not depend on the fractional order α. For the integrals involving time power functions of order ±α a distinction has to be made between fractional (0 < α < 1) and classical (α = 1) viscoelasticity. In the fractional case, these time integrals are related to the asymptotic values of the fractional derivative or integral of the overall relaxation function L(t) which have no geometric interpretation. On the one hand, the fractional derivative of the macroscopic relaxation function reads where use has been made of the property (B13) on Mittag-Leffler functions. At short times, we thus obtain In the case of a relaxation loading test with constant strain ε 0 , this implies, for a classical viscoelastic behaviour (α = 1), that lim t→0 +σ On the other hand, the fractional integral of the effective relaxation function reads which gives at long times For classical viscoelasticity, in the case of a loading test with a constant macroscopic strain rateε 0 , it follows that Relations (29) and (32) provide the physical meaning of the two time integrals of the relaxation spectrum +∞ 0 [19].

Approximation of the effective relaxation spectrum G
To derive an estimate of the overall viscoelastic relaxation function, it is necessary to approximate the unknown relaxation spectrum G(τ ). The most usual approaches make use of box functions, wedge functions or a sum of Dirac delta functions [14,60]. The latter case corresponds to a line spectrum which results in the relaxation function It corresponds to a generalized Maxwell fractional calculus model [31]. With the property of the Mittag-Leffler functions E 1 (t) = e t ("Appendix B.3"), the form (34) reduces to a usual Prony series in the case of classical viscoelasticity. As a side note, it must be mentioned that in particular cases the exact effective relaxation spectrum obtained by homogenization is indeed composed of discrete lines; see [2,51] for mixtures of Maxwell constituents. Similar results hold if we consider instead Zener constituents. The storage and loss moduli tensors corresponding to the fractional Prony approximation (34) are thus given by with q k = 1 + 2(ωτ k ) α cos πα 2 + (ωτ k ) 2α . From the approximation (33) and the conditions on time integrals of the effective relaxation spectrum (25), it follows that This set of three relations can be used to construct a minimal approximation of the overall relaxation spectrum (33) for mixtures of fractional Zener constituents.

Approximate model for isotropic composites
The previous tensorial conditions are now specified for composites with overall and local isotropy. In this particular case, the effective relaxation spectrum can be written with J and K the isotropic projectors on hydrostatic and deviatoric symmetric second-order tensors. Moreover, attention is further restricted to composite materials whose local bulk relaxation is negligible. It is thus assumed that the viscoelastic constituents are purely elastic in dilatation and viscoelastic in shear. The relaxation function of a phase (s) is thus given by As it is well known, the overall response is in general viscoelastic both in dilatation and shear [26]. The overall bulk viscoelastic response arises because of the contrast between the elastic bulk moduli of the phases which implies a deviation of the local fields from a purely hydrostatic state. (For porous plastic materials, the same argument explains why plasticity occurs under a macroscopic pressure loading.) However, the bulk viscoelasticity of the composite remains rather small and will be neglected in the sequel (κ k = 0), whence the approximate form of the effective relaxation function Since only shear relaxation times are now involved, we will skip the upper index μ and use the notation τ k hereafter.
In elastic composites, the second moment per phase of the strain field is classically obtained from the partial derivatives of the overall energy with respect to the phase elastic moduli [4,33,48]. In the particular case of local and overall isotropy, its deviatoric and hydrostatic parts read with ε 2 eq = 2 3 K::(ε ⊗ ε) and tr(ε) 2 = 3 J::(ε ⊗ ε). Given the assumptions made on the local and overall viscoelastic moduli, only the relation on the phase average of the square of the equivalent strain (40) 1 is useful for the present study. With (37) and (40), the tensorial relations (36) The minimal number of terms required in the generalized Prony series to fulfil the three relations (41) is thus K = 2. To determine the four parameters μ 1 , μ 2 , τ α 1 and τ α 2 , an additional relation is necessary. Note that this is not the case for mixtures of Maxwell constituents for which the number of terms in the serie can be chosen so that the number of unknown parameters match the number of available equations [6,58]. This minimal approximation is assessed in the sequel for particulate two-phase composites.

Approximate viscoelastic model parameters
The composites considered are made of a viscoelastic matrix (phase 1) reinforced by spherical elastic inclusions randomly distributed. The system of Eq. (41) becomes In order to solve the system (42), we set a fractional relaxation time of the serie equal to the one of the matrix phase: τ α 1 = τ (1) α . The three other parameters are then given by

Mean-field homogenization estimates
The approximate viscoelastic model only requires to estimate the asymptotic effective elastic moduli of the composite, which correspond to the relaxed and glassy states, and their derivatives with respect to the phase moduli. In the following, this is achieved by resorting to a mean-field homogenization scheme whose relevance is assessed with reference numerical estimates on representative microstructures. The classical estimates (bounds) of the Hashin-Shtrikman type for materials with an isotropic distribution of the phases [63] read (1) μ + μ (1) and κ = κ (1) + c 2 κ (2) − κ (1) 1 + c 1 κ (2) − κ (1) κ + κ (1) with μ and κ the shear and bulk moduli of the Hill constraint tensor μ = μ 0 9κ 0 + 8μ 0 The choice (μ 0 = μ (1) , κ 0 = κ (1) ) corresponds to the lower Hashin-Shtrikman bound which coincides with the Mori-Tanaka (MT) model [42] while the choice (μ 0 = μ, κ 0 = κ) defines the self-consistent (SC) estimate [30]. Another widely used model to estimate the effective response of particulate composites is the generalized self-consistent (GSC) model [11,28] which considers coated spherical inclusions. The effective bulk modulus is the one of the Hashin composite sphere assemblage which attains the lower Hashin-Shtrikman bound, and the effective shear modulus is solution of a quadratic equation. It has been shown that this model correctly describes experimental results on the effective (shear) viscosity of polydisperse suspensions with rigid particles [10]. It is also noted that comparisons have been reported between the MT, SC, GSC models and unit-cell computations, in the elastic case, for composites with monodisperse spherical inclusions [22].

Fourier transform-based numerical homogenization
To assess viscoelastic mean-field estimates for particulate composites, full-field computations are performed on unit cells containing a random distribution of polydisperse spherical inclusions. Thanks to the correspondence principle [38], the overall complex viscoelastic properties can be obtained by making use of the fast Fourier transform (FFT) numerical scheme [17,43,44]. We consider a symbolic heterogeneous elastic problem (6), with periodic boundary conditions, which corresponds to the stationary harmonic regime at angular frequency ω. By introducing a homogeneous reference material L 0 , the constitutive law can be rewritten as with τ * the polarization field. The complex strain solution field thus reads, in real and Fourier spaces, with ξ the wave vector. 0 represents the Fourier transform of the strain Green operator associated with the reference medium L 0 with the acoustic tensor κ = ξ .L 0 .ξ . The notation [.] s indicates minor and major symmetrizations. To solve the Lippman-Schwinger equation for complex strain field ε * (48) 1 , we resort to the iterative numerical scheme of Eyre and Milton which is well-suited for materials with a high mechanical contrast [15,43]. It must be noted that, as pointed out by Figliuzzi et al. [17], the reference elastic stiffness L 0 must be real to ensure the symmetry of the physical strain and stress fields (i.e. real parts of the corresponding complex fields). For the studied two-phase composite, the optimal choice for the shear and bulk moduli of an isotropic reference medium L 0 is μ 0 = Re μ (1) Re μ (2) and κ 0 = Re κ (1) Re κ (1) .

Assessment of mean-field models with respect to FFT reference results
We consider two-phase particulate composites with a size polydisperse distribution of elastic spherical inclusions embedded in a fractional Zener viscoelastic matrix. The material parameters are given in Table 1. The elastic moduli, used in [25], correspond to silica particles and a typical epoxy resin matrix. It is noted that the elastic shear moduli contrast vary from 30, in the glassy state (ω → +∞), to 3000, in the relaxed state (ω → 0). Distributions of particles in a cubic unit cell, with geometric periodic conditions, are obtained by using the Random Sequential Adsorption (RSA) algorithm [59,62] ("Appendix C"). Full-field FFT simulations, for a prescribed angular frequency ω, are performed on discretized unit cells with a regular grid of 255 3 voxels. The overall response is averaged over 10 different microstructural realizations for each volume fraction (c = 0.1, 0.3 and 0.5). Besides, since the bulk viscoelasticity is neglected, only shear loadings are considered. An overall isotropic complex shear moduli μ * α (iω) = μ α (ω)+i μ α (ω) is defined by averaging the response of the composite subjected to three independent shear loadings. Comparisons between MT, SC, GSC estimates and FFT computations of the overall storage shear moduli μ α (ω) in the range ω ∈ [10 −2 ; 10 −4 ] are reported in Fig. 1.
At low particles volume fraction, the three mean-field models almost coincide and agree with the numerical results. With increasing volume fraction, only the GSC scheme provides accurate estimates. These results are consistent with the ones obtained in [25] by using finite-element computations with a time-integration approach. It can be also noted that the SC scheme, which presents a percolation threshold for an infinite mechanical contrast, leads to an unrealistic quasi-elastic effective response at volume fraction c = 0.5.

Approximate viscoelastic GSC estimates
Following these assessments, we choose the GSC scheme to build an approximate viscoelastic homogenization model which makes only use of the linear homogenization elastic problems at low and high frequencies. The with q k = 1 + 2(ωτ k ) α cos πα 2 + (ωτ k ) 2α . The unknown parameters of the serie expansions are obtained with relations (43) and (44) where the closed form GSC estimate [11,28] is used for μ e r and μ e g . The comparison with reference FFT computations shows an overall good agreement of the approximate GSC model (Fig. 2). It is worth mentioning that the exact shear relaxation spectrum of the GSC model presents a continuous part [2].

Conclusion
This study is a contribution to the description of the overall properties of composites materials with fractional viscoelastic constituents. Previous results for classical viscoelasticity [19] have been extended to fractional behaviour in the particular case of mixtures of Zener constituents. Exact relations have been obtained on time integrals of the effective relaxation spectrum which are related to non-integer derivative and integral of the complex viscoelastic moduli. However, by contrast with classical viscoelasticity, no simple physical interpretation can be given. By making use of these relations, an approximate model has been proposed for the effective complex moduli by describing the relaxation spectrum of the composite with only two Dirac delta functions. The parameters of this approximation (relaxation times and corresponding intensities) depend on the local properties of the constituents and the intraphase second moments of the strain fields which are solution of the purely elastic heterogeneous problems at low and high frequencies. For particulate composites with polydisperse spherical elastic inclusions, the comparisons of this approximate model with reference FFT full-field simulations show that it delivers accurate estimates in the whole frequency range.

A Stieltjes convolution and Laplace-Carson transform
The Stieltjes convolution product of two functions f and g is the derivative of their usual convolution product. If g is a differentiable function of time, it reads If g is only piecewise continuous and differentiable, its time derivative in (A1) contains Dirac masses at discontinuity points t n and the Stieltjes convolution product must be understood as where [g] n is the discontinuity of g at time t n andġ(u) is the usual derivative of g. The Laplace-Carson (LC) transform of a function f (t) is defined by From (A1) and (A3), it follows that

B Useful relations for fractional calculus
For a comprehensive review on fractional calculus and its applications in viscoelasticity, the reader is referred to [13,36,46,47].

B.1 Gamma function
The Gamma function Γ is a continuation of the factorial function to complex numbers. It is defined by Integration by parts leads to the following property It is also noted that Γ (1) = 1 and Γ (n) = (n − 1)!, ∀n ∈ N.

B.2 Fractional integral and derivative operators
Cauchy's formula for successive integrations of a causal function f , with null initial conditions, reads namely By using the function Γ , the integral operator I can be extended in a straightforward manner to non-integer order α > 0 It can be also noted that The non-integer derivative operator D α is defined by successive derivation (of order 1) and integration of order 1 − α, that is From (A4), the LC transform of the fractional derivative thus simply reads which extends a classical result to the case of derivatives with non-integer order.

B.3 Mittag-Leffler function
The Mittag-Leffler function (or exponential) E α is defined by the following power series Alternatively, it admits an integral representation [24,37] which reads It can be also noted that +∞ 0 H α (θ ) dθ = 1 since E α (0 + ) = 1. The Mittag-Leffler function corresponds to a generalization of the exponential function since E 1 (z) = e z . Besides, its fractional derivative satisfies The function E α naturally appears in fractional calculus since it is solution of the fractional differential equation The solution reads [13] f (t) = E α (−(λt) α ). The asymptotic expansion of the function E α (−t α ) at short and long times reads The constitutive relation of a 1D fractional element reads [7] σ (t) = Eτ α D α ε(t) (B17) with elastic modulus E and fractional relaxation time τ α . From the definition (B7), it follows that where the relaxation function R(t) reads (B19) The inversion of the constitutive law (B17) simply yields (B20) From relation (B6), it may be written as where the creep (or retardation) function F(t) reads (B22)

C Random distribution of size polydisperse particles
Random polydisperse microstructures have been built using the RSA algorithm which consists in placing randomly, irreversibly and sequentially non-overlapping geometric objects into a fixed volume. The size distribution of the spherical inclusions follows a lognormal density function φ(r ) with r the radius of the particles φ(r ) = 1 r σ √ 2π e −(ln(r )−μ) 2 /2σ 2 , σ > 0, (C1) where the parameters μ and σ represent the mean and the standard deviation of the variable's natural logarithm, respectively. The generation of a microstructure composed of P spherical inclusions requires to define inclusion size families characterized by a fixed radius and a number of inclusions. The determination of the inclusion size families and the associated radii, supposed to be evenly spaced hereafter, is assessed from the integral of the lognormal density function A(+∞) = +∞ 0 φ(r ; μ, σ ) dr = 1. In practice, it is thus necessary to choose a maximal radius r max corresponding to a given value A(r max ). We have chosen a maximal radius corresponding to A(r max ) = 0.99 r max = e σ ψ+μ with ψ = √ 2π erf −1 22 25 (C2) By considering evenly spaced radii, the knowledge of the limit radius allows to determine each radius associated with an inclusion size family. The number of inclusions per inclusion size family is then deduced by the calculation of the area sections. For n inclusion size families, the radius r i and the number of particles P i associated with an inclusion size family i read ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ r i = r max 2n (2i − 1) , n > 2, An example of polydisperse microstructure is reported in Fig. 3.