<figure><img src="image_1.png" /></figure>

HAL
open science

# Homogenization estimates for the effective response of fractional viscoelastic particulate composites

Valentin Gallican, Renald Brenner

► To cite this version:

Valentin Gallican, Renald Brenner. Homogenization estimates for the effective response of fractional viscoelastic particulate composites. Continuum Mechanics and Thermodynamics, 2019, 31, pp.823-840. 10.1007/s00161-018-0741-8 . hal-02023868

HAL Id: hal-02023868

https://hal.science/hal-02023868

Submitted on 18 Feb 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Homogenization estimates for the effective response of fractional viscoelastic particulate composites

Valentin Gallican · Renald Brenner

the date of receipt and acceptance should be inserted later

**Abstract** 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 frequency. 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.

**Keywords** homogenization · particulate composite · fractional viscoelasticity · polymer material

## 1 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

V. Gallican · R. Brenner (✉)

Sorbonne Université, CNRS, Institut Jean le Rond d'Alembert, F-75005, Paris, France

E-mail: renald.brenner@sorbonne-universite.fr
(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 [52,18,57]. It can be evidenced by making use
of the correspondence principle [38,26] 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 serie [53,34] 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 [61,40,49].

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 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 [45,20,21] 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, [50,12,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 [58,6,19]. 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 [58,6].
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 [44,43] performed on polydisperse microstructures.

## 2 Overall harmonic response of viscoelastic composites

### 2.1 Viscoelastic constitutive relations

#### 2.1.1 Classical viscoelasticity

From the Boltzmann superposition principle, the stress response of a non-ageing linear viscoelastic material to a given derivable strain loading path $\varepsilon(u)$, $u \in [0; t]$ with additional discontinuities $[\varepsilon]_i$ at times $t_i$, and initial conditions $\sigma(t = 0) = 0$, classically reads

$$ \sigma(t) = \int_{0}^{t} \mathbf{L}(t-u) : \dot{\varepsilon}(u) \, du + \sum_{i} \mathbf{L}(t-t_{i}) : [\varepsilon]_{i}^{1} \quad (1) $$

with $\mathbf{L}(t)$ the viscoelastic stiffness tensor (i.e. relaxation function) or, with a concise notation,

$$ \sigma(t) = \frac{d}{dt} (\mathbf{L} * \varepsilon)(\mathbf{x}, t) = (\mathbf{L} \otimes \varepsilon)(\mathbf{x}, t) \quad (2) $$

<footnote n="1">¹ Throughout the text, the tensors are indicated in bold notation. Double and quadruple-dot products between tensors **a** and **b** are respectively denoted by **a** : **b** and **a** :: **b**.</footnote>
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

$$ \mathbf{L}(t) = \mathbf{L}_{e_r} + \int_{0}^{+\infty} \mathbf{G}(\tau)e^{-t/\tau} d\tau. \quad (3) $$

$\mathbf{L}_{e_r}$ denotes the relaxed elastic stiffness tensor, by reference to polymer materials, while $\mathbf{G}$ and $\tau$ are the relaxation spectrum and time.

### 2.1.2 Fractional viscoelasticity

Differential constitutive laws with non-integer derivatives (i.e. fractional calculus [46, 47, 13]) 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 $]-\infty, t]$ (Appendix B.2). The non-integer derivative order $\alpha$ 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 $\sigma$ and the fractional derivative of the strain $\mathcal{D}^{\alpha}\varepsilon$ with $\mathcal{D}^{\alpha} \equiv d^{\alpha}/dt^{\alpha}$ and $0 < \alpha < 1$ (Appendix B.4). Its relaxation function can be described as a sum of exponentials weighted with a power-law relaxation spectrum [60,35]. 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 [55,29].

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 $\mathbf{L}_{\alpha}(t)$ given by [31,32]

$$ \mathbf{L}_{\alpha}(t) = \mathbf{L}_{e_r} + \int_{0}^{+\infty} \mathbf{G}(\tau)E_{\alpha}[-(t/\tau)^{\alpha}] d\tau \quad (4) $$

where $E_{\alpha}(t)$ is the Mittag-Leffler function (Appendix B.3). Obviously, when $\alpha = 1$ the relaxation function $\mathbf{L}_{\alpha}(t)$ corresponds to the classical relaxation function $\mathbf{L}(t)$ (3). In the case of a relaxation spectrum $\mathbf{G}(\tau)$ 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 fulfills 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.

## 2.2 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 $\Omega$ which consists of $N$ homogeneous phases with characteristic function $\chi^{(s)}(\mathbf{x})$ and volume $\Omega^{(s)}$ ($s \in [0; N]$). Moreover, it is assumed that $\Omega^{(s)} \ll \Omega$ and that the phases are perfectly bonded. The fractional viscoelastic relaxation function of phase (s) is denoted $L_{\alpha}^{(s)}(t)$. It follows that the pointwise viscoelastic relaxation tensor $L_{\alpha}(\mathbf{x}, t)$ reads

$$L_{\alpha}(\mathbf{x}, t) = \sum_{s=1}^{N} L_{\alpha}^{(s)}(t) \chi^{(s)}(\mathbf{x}) \quad (5)$$

with $\chi^{(s)}(\mathbf{x}) = 1$ if $\mathbf{x} \in \Omega^{(s)}$ and 0 otherwise. The volume averages over $\Omega$ and $\Omega^{(s)}$ are respectively denoted $\langle \bullet \rangle$ and $\langle \bullet \rangle^{(s)}$. By definition of the characteristic function, the volume fraction of phase (s) is $c_s = \langle \chi^{(s)} \rangle$.

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\omega$ ($i^2 = -1$) [27]. By assuming an overall strain loading $\bar{\epsilon}(t) = \bar{\epsilon}^* e^{i\omega t}$, the local problem corresponding to the steady-state regime at angular frequency $\omega$ reads

$$\left\{
\begin{array}{l}
\sigma^*(\mathbf{x}, i\omega) = L_{\alpha}^*(\mathbf{x}, i\omega) : \varepsilon^*(\mathbf{x}, i\omega), \quad \forall (\mathbf{x}, \omega) \in \Omega \times [0; +\infty[, \\
\text{div } \sigma^* = 0, \quad \text{curl}(\text{curl} \, \varepsilon^*) = 0, \quad \forall (\mathbf{x}, \omega) \in \Omega \times [0; +\infty[, \\
\langle \varepsilon^* \rangle = \bar{\varepsilon}^*,
\end{array}
\right.
\quad (6)$$

with prescribed boundary conditions. The complex fields $(\sigma^*, \varepsilon^*, L_{\alpha}^*)$ are the time LC transforms of the fields $(\sigma, \varepsilon, L_{\alpha})$.
## 2.3 Viscoelastic composites with fractional Zener constituents

### 2.3.1 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

$$ \sigma(t) + \mathbf{L}_f : (\mathbf{L}_{e_g} - \mathbf{L}_{e_r})^{-1} : \mathcal{D}^{\alpha} \sigma(t) = \mathbf{L}_{e_r} : \varepsilon(t) + \mathbf{L}_{e_g} : \mathbf{L}_f : (\mathbf{L}_{e_g} - \mathbf{L}_{e_r})^{-1} : \mathcal{D}^{\alpha} \varepsilon(t), \quad 0 < \alpha < 1, \quad (7) $$

with $\mathbf{I}$ the fourth-order identity tensor. It exhibits asymptotic elastic behaviours with glassy moduli $\mathbf{L}_{e_g}$ at short times ($t \to 0$) and relaxed moduli $\mathbf{L}_{e_r}$ at long times ($t \to +\infty$). $\mathbf{L}_f$ represents the constitutive fractional viscous modulus. The viscoelastic relaxation tensor of phase ($s$), defined by relation (7), reads

$$ \mathbf{L}_{\alpha}^{(s)}(t) = \mathbf{L}_{e_r}^{(s)} + \mathbf{G}^{(s)} E_{\alpha} \left[ -(t/\tau^{(s)})^{\alpha} \right] \quad \text{with} \quad \mathbf{G}^{(s)} = \mathbf{L}_{e_g}^{(s)} - \mathbf{L}_{e_r}^{(s)}. \quad (8) $$

Also, it is worth noting that the eigenvalues of $\mathbf{L}_f^{(s)} : (\mathbf{G}^{(s)})^{-1}$ correspond to the fractional relaxation times $\tau^{\alpha(s)}$ (units: $s^\alpha$) of the constituent. Following [31], it is assumed that, in the general case, the effective relaxation tensor $\tilde{\mathbf{L}}(t)$ can be written in the form

$$ \tilde{\mathbf{L}}_{\alpha}(t) = \tilde{\mathbf{L}}_{e_r} + \int_{0}^{+\infty} \tilde{\mathbf{G}}(\tau) E_{\alpha} \left[ -(t/\tau) \right] d\tau. \quad (9) $$

As a consequence, the complex viscoelastic relaxation tensor $\tilde{\mathbf{L}}_\alpha^*$, which characterizes the steady-state harmonic regime at angular frequency $\omega$, is given by

$$ \tilde{\mathbf{L}}_{\alpha}^{*}(i\omega) = \mathcal{L}\mathcal{C}(\tilde{\mathbf{L}}_{\alpha}(t)) = \tilde{\mathbf{L}}_{e_r} + \int_{0}^{+\infty} \frac{(i\omega\tau)^{\alpha}}{1+(i\omega\tau)^{\alpha}} \tilde{\mathbf{G}}(\tau) \, d\tau. \quad (10) $$

The overall complex constitutive relation is

$$ \langle \sigma^* \rangle(i\omega) = \tilde{\mathbf{L}}_{\alpha}^*(i\omega) : \langle \varepsilon^* \rangle(i\omega), \quad \forall \omega \in [0; +\infty[, \quad (11) $$

and the complex relaxation tensor admits the decomposition

$$ \tilde{\mathbf{L}}_{\alpha}^{*}(i\omega) = \tilde{\mathbf{L}}_{\alpha}'(\omega^{\alpha}) + i\tilde{\mathbf{L}}_{\alpha}''(\omega^{\alpha}) \quad (12) $$
with $(\tilde{\mathbf{L}}'_{\alpha}, \tilde{\mathbf{L}}''_{\alpha})$ the overall storage and loss moduli which are respectively proportional to the stored and dissipated energies [60]. By noting that $i^{\alpha} = e^{i\pi\alpha/2}$, they read

$$ \left\{ \begin{aligned} \tilde{\mathbf{L}}'_{\alpha}(\omega^{\alpha}) &= \tilde{\mathbf{L}}_{e_r} + \int_{0}^{+\infty} \frac{1}{q} \left[ (\omega\tau)^{\alpha} \cos\left(\frac{\pi\alpha}{2}\right) + (\omega\tau)^{2\alpha} \right] \tilde{\mathbf{G}}(\tau) \, d\tau, \\ \tilde{\mathbf{L}}''_{\alpha}(\omega^{\alpha}) &= \int_{0}^{+\infty} \frac{1}{q} \left[ (\omega\tau)^{\alpha} \sin\left(\frac{\pi\alpha}{2}\right) \right] \tilde{\mathbf{G}}(\tau) \, d\tau, \end{aligned} \right. \qquad (13) $$

with $q = 1 + 2(\omega\tau)^{\alpha} \cos(\frac{\pi\alpha}{2}) + (\omega\tau)^{2\alpha}$. The effective loss tensor $\tilde{\eta}(\omega^{\alpha})$ which characterizes the damping is classically defined by

$$ \tilde{\eta}(\omega^{\alpha}) = \tilde{\mathbf{L}}''_{\alpha}(\omega^{\alpha}) : [\tilde{\mathbf{L}}'_{\alpha}(\omega^{\alpha})]^{-1}. \qquad (14) $$

Finally, it is important to note that as $\omega \to +\infty$ or $\omega \to 0$ the local fields are asymptotically solutions of purely elastic heterogeneous problems. The pointwise complex strain field $\varepsilon^*(\mathbf{x}, i\omega)$ thus satisfies

$$ \lim_{\omega \to +\infty} \varepsilon^*(\mathbf{x}, i\omega) = \varepsilon_g(\mathbf{x}) \quad \text{and} \quad \lim_{\omega \to 0} \varepsilon^*(\mathbf{x}, i\omega) = \varepsilon_r(\mathbf{x}) \qquad (15) $$

with $\varepsilon_g(\mathbf{x})$ and $\varepsilon_r(\mathbf{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 $\sigma^*(\mathbf{x}, i\omega)$ with asymptotic fields $\sigma_g(\mathbf{x})$ and $\sigma_r(\mathbf{x})$ respectively.

### 2.3.2 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

$$ \left( \begin{matrix} \boldsymbol{\sigma}' \\ -\boldsymbol{\sigma}'' \end{matrix} \right) = \mathbb{L}_{\mathrm{R}} : \left( \begin{matrix} \boldsymbol{\varepsilon}' \\ \boldsymbol{\varepsilon}'' \end{matrix} \right), \quad \mathbb{L}_{\mathrm{R}} = \left( \begin{matrix} \mathbf{L}'_{\alpha} & -\mathbf{L}''_{\alpha} \\ -\mathbf{L}''_{\alpha} & -\mathbf{L}'_{\alpha} \end{matrix} \right) \quad \text{and} \quad \left( \begin{matrix} \boldsymbol{\sigma}' \\ \boldsymbol{\sigma}'' \end{matrix} \right) = \mathbb{L}_{\mathrm{I}} : \left( \begin{matrix} \boldsymbol{\varepsilon}'' \\ \boldsymbol{\varepsilon}' \end{matrix} \right), \quad \mathbb{L}_{\mathrm{I}} = \left( \begin{matrix} -\mathbf{L}''_{\alpha} & \mathbf{L}'_{\alpha} \\ \mathbf{L}'_{\alpha} & \mathbf{L}''_{\alpha} \end{matrix} \right) \qquad (16) $$

where $(\boldsymbol{\sigma}', \boldsymbol{\varepsilon}', \mathbf{L}'_{\alpha})$ and $(\boldsymbol{\sigma}'', \boldsymbol{\varepsilon}'', \mathbf{L}''_{\alpha})$ are real fields, that is

$$ \boldsymbol{\sigma}^*(i\omega) = \boldsymbol{\sigma}'(\omega^\alpha) + i\boldsymbol{\sigma}''(\omega^\alpha), \quad \varepsilon^*(i\omega) = \varepsilon'(\omega^\alpha) + i\varepsilon''(\omega^\alpha) \quad \text{and} \quad \mathbf{L}_\alpha^*(i\omega) = \mathbf{L}_\alpha'(\omega^\alpha) + i\mathbf{L}_\alpha''(\omega^\alpha) \qquad (17) $$

where the position vector $\mathbf{x}$ has been omitted for conciseness. By considering the frequency-dependent complex “energy” $\phi^*$

$$ \phi^*(\omega^\alpha) = \boldsymbol{\sigma}^* : \boldsymbol{\varepsilon}^* = (\boldsymbol{\sigma}' : \boldsymbol{\varepsilon}' - \boldsymbol{\sigma}'' : \boldsymbol{\varepsilon}'') + i (\boldsymbol{\sigma}'' : \boldsymbol{\varepsilon}' + \boldsymbol{\sigma}' : \boldsymbol{\varepsilon}''), \qquad (18) $$
two saddle-point variational principles on the effective complex “energy” $\tilde{\phi}^*$ can be established [9,19]

$$ \mathrm{Re} \left[ \tilde{\phi}^*(\omega^\alpha, \bar{\epsilon}', \bar{\epsilon}'') \right] = \min_{\epsilon', \langle \epsilon' \rangle = \bar{\epsilon}'} \max_{\epsilon'', \langle \epsilon'' \rangle = \bar{\epsilon}''} \langle \sigma' : \epsilon' - \sigma'' : \epsilon'' \rangle = \min_{\epsilon', \langle \epsilon' \rangle = \bar{\epsilon}'} \max_{\epsilon'', \langle \epsilon'' \rangle = \bar{\epsilon}''} \left\langle \begin{pmatrix} \epsilon' \\ \epsilon'' \end{pmatrix} : \mathbb{L}_R : \begin{pmatrix} \epsilon' \\ \epsilon'' \end{pmatrix} \right\rangle \quad (19) $$

and

$$ \mathrm{Im} \left[ \tilde{\phi}^*(\omega^\alpha, \bar{\epsilon}', \bar{\epsilon}'') \right] = \min_{\epsilon', \langle \epsilon' \rangle = \bar{\epsilon}'} \max_{\epsilon'', \langle \epsilon'' \rangle = \bar{\epsilon}''} \langle \sigma'' : \epsilon' + \sigma' : \epsilon'' \rangle = \min_{\epsilon', \langle \epsilon' \rangle = \bar{\epsilon}'} \max_{\epsilon'', \langle \epsilon'' \rangle = \bar{\epsilon}''} \left\langle \begin{pmatrix} \epsilon'' \\ \epsilon' \end{pmatrix} : \mathbb{L}_I : \begin{pmatrix} \epsilon'' \\ \epsilon' \end{pmatrix} \right\rangle. \quad (20) $$

The solution fields are functions of $\omega^\alpha$. 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 $\tilde{\phi}^d = (\omega/2) \langle \sigma'' : \epsilon' - \sigma' : \epsilon'' \rangle$. However, Re($\tilde{\phi}^*$) and Im($\tilde{\phi}^*$) are stationary with respect to the strain fields $\epsilon'$ and $\epsilon''$. 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

$$ \frac{\partial}{\partial h} \mathrm{Re}(\tilde{\phi}^*) = \left\langle \begin{pmatrix} \epsilon' \\ \epsilon'' \end{pmatrix} : \frac{\partial \mathbb{L}_R}{\partial h} : \begin{pmatrix} \epsilon' \\ \epsilon'' \end{pmatrix} \right\rangle \quad \text{and} \quad \frac{\partial}{\partial h} \mathrm{Im}(\tilde{\phi}^*) = \left\langle \begin{pmatrix} \epsilon'' \\ \epsilon' \end{pmatrix} : \frac{\partial \mathbb{L}_I}{\partial h} : \begin{pmatrix} \epsilon'' \\ \epsilon' \end{pmatrix} \right\rangle. \quad (21) $$

### 2.3.3 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 $\langle \epsilon' \rangle = \bar{\epsilon}$ and $\langle \epsilon'' \rangle = 0$, the stationary principle on the real part of $\tilde{\phi}^*$ (19) leads to

$$ \left\{ \begin{array}{l} \lim_{\omega \to 0} \bar{\epsilon} : \tilde{\mathbf{L}}'_{\alpha}(\omega^{\alpha}) : \bar{\epsilon} = \lim_{\omega \to 0} \bar{\epsilon} : \tilde{\mathbf{L}}_{e_r} : \bar{\epsilon} = \sum_s c_s \mathbf{L}_{e_r}^{(s)} :: \langle \epsilon_r \otimes \epsilon_r \rangle^{(s)}, \\ \\ \lim_{\omega \to +\infty} \bar{\epsilon} : \tilde{\mathbf{L}}'_{\alpha}(\omega) : \bar{\epsilon} = \lim_{\omega \to +\infty} \bar{\epsilon} : \tilde{\mathbf{L}}_{e_g} : \bar{\epsilon} = \sum_s c_s \mathbf{L}_{e_g}^{(s)} :: \langle \epsilon_g \otimes \epsilon_g \rangle^{(s)}. \end{array} \right. \quad (22) $$

On the other hand, the stationary principle on the imaginary part of $\tilde{\phi}^*$ (20) implies that

$$ \lim_{\omega \to 0} \tilde{\mathbf{L}}''_{\alpha}(\omega^{\alpha}) = \lim_{\omega \to +\infty} \tilde{\mathbf{L}}''_{\alpha}(\omega^{\alpha}) = 0. \quad (23) $$
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 $\tilde{\mathbf{L}}_{e_r}$ and $\tilde{\mathbf{L}}_{e_g}$ thus simply requires to solve the heterogeneous elastic problems respectively in the relaxed ($\omega \to 0$) and glassy ($\omega \to +\infty$) states. Besides, by choosing $h = \omega^\alpha$, the derivatives of the stationary principles (21) provide the results

$$ \left\{ \begin{aligned} \lim_{\omega \to 0} \bar{\varepsilon} : \frac{\partial \tilde{\mathbf{L}}''_{\alpha}}{\partial \omega^{\alpha}}(\omega^{\alpha}) : \bar{\varepsilon} &= \tan\left(\frac{\pi\alpha}{2}\right) \lim_{\omega \to 0} \bar{\varepsilon} : \frac{\partial \tilde{\mathbf{L}}'_{\alpha}}{\partial \omega^{\alpha}}(\omega^{\alpha}) : \bar{\varepsilon} = \sin\left(\frac{\pi\alpha}{2}\right) \sum_s c_s \mathbf{L}_f^{(s)} :: \langle \varepsilon_r \otimes \varepsilon_r \rangle^{(s)} \\ \lim_{\omega \to +\infty} (i\omega)^{2\alpha} \bar{\varepsilon} : \frac{\partial \tilde{\mathbf{L}}''_{\alpha}}{\partial \omega^{\alpha}}(\omega^{\alpha}) : \bar{\varepsilon} &= \tan\left(\frac{\pi\alpha}{2}\right) \lim_{\omega \to +\infty} (i\omega)^{2\alpha} \bar{\varepsilon} : \frac{\partial \tilde{\mathbf{L}}'_{\alpha}}{\partial \omega^{\alpha}}(\omega^{\alpha}) : \bar{\varepsilon} = \sin\left(\frac{\pi\alpha}{2}\right) \sum_s c_s (\mathbf{G}^{(s)} : \mathbf{L}_f^{(s)}^{-1} : \mathbf{G}^{(s)}) :: \langle \varepsilon_g \otimes \varepsilon_g \rangle^{(s)}. \end{aligned} \right. \quad (24) $$

It can be remarked that unlike classical viscoelasticity [19] the asymptotic values of the derivative of the effective storage modulus $\tilde{\mathbf{L}}'_\alpha$ do not vanish since $0 < \alpha < 1$. These asymptotic values of the derivative of the real and imaginary parts of the overall complex “energy” $\phi^*$ 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) frequency.

### 2.3.4 Consequences on the overall relaxation spectrum and their physical meaning

With the general expression of the overall complex relaxation function $\tilde{\mathbf{L}}^*_\alpha(i\omega)$ (10), the asymptotic relations (22–24) have implications on time integrals of the overall relaxation spectrum $\tilde{\mathbf{G}}(\tau)$

$$ \left\{ \begin{aligned} & \int_{0}^{+\infty} \tilde{\mathbf{G}}(\tau) \, d\tau = \tilde{\mathbf{L}}_{e_g} - \tilde{\mathbf{L}}_{e_r}, \\ & \bar{\varepsilon} : \int_{0}^{+\infty} \tau^\alpha \tilde{\mathbf{G}}(\tau) \, d\tau : \bar{\varepsilon} = \sum_s c_s \mathbf{L}_f^{(s)} :: \langle \varepsilon_r \otimes \varepsilon_r \rangle^{(s)}, \\ & \bar{\varepsilon} : \int_{0}^{+\infty} \tau^{-\alpha} \tilde{\mathbf{G}}(\tau) \, d\tau : \bar{\varepsilon} = \sum_s c_s (\mathbf{G}^{(s)} : \mathbf{L}_f^{(s)}^{-1} : \mathbf{G}^{(s)}) :: \langle \varepsilon_g \otimes \varepsilon_g \rangle^{(s)}. \end{aligned} \right. \quad (25) $$

The first relation have a simple and well-known physical interpretation in the case of a relaxation loading test ($\bar{\varepsilon}(t) = \bar{\varepsilon}_0$). Indeed, we have

$$ \lim_{t \to 0^+} \bar{\sigma}(t) = \left[ \tilde{\mathbf{L}}_{e_r} + \int_{0}^{+\infty} \tilde{\mathbf{G}}(\tau) \, d\tau \right] : \bar{\varepsilon}_0. \quad (26) $$
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 $\alpha$. For the integrals involving time power functions of order $\pm\alpha$ a distinction has to be made between fractional ($0 < \alpha < 1$) and classical ($\alpha = 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 $\tilde{\mathbf{L}}(t)$ which have no geometric interpretation. On one hand, the fractional derivative of the macroscopic relaxation function reads

$$ \mathcal{D}^{\alpha} \tilde{\mathbf{L}}(t) = - \int_{0}^{+\infty} \tau^{-\alpha} \tilde{\mathbf{G}}(\tau) E_{\alpha} [(-t/\tau)^{\alpha}] d\tau \quad (27) $$

where use has been made of the property (B13) on Mittag-Leffler functions. At short times, we thus obtain

$$ \lim_{t \to 0^+} \mathcal{D}^{\alpha} \tilde{\mathbf{L}}(t) = - \int_{0}^{+\infty} \tau^{-\alpha} \tilde{\mathbf{G}}(\tau) d\tau. \quad (28) $$

In the case of a relaxation loading test with constant strain $\tilde{\varepsilon}_0$, this implies, for a classical viscoelastic behaviour ($\alpha = 1$), that

$$ \lim_{t \to 0^+} \dot{\tilde{\sigma}}(t) = - \int_{0}^{+\infty} \tau^{-1} \tilde{\mathbf{G}}(\tau) d\tau : \tilde{\varepsilon}_0. \quad (29) $$

On the other hand, the fractional integral of the effective relaxation function reads

$$ \mathcal{I}^{\alpha} \tilde{\mathbf{L}}(t) = \tilde{\mathbf{L}}_{e_r} \frac{t^{\alpha}}{\Gamma(\alpha + 1)} + \int_{0}^{+\infty} \tau^{\alpha} \tilde{\mathbf{G}}(\tau) (1 - E_{\alpha} [(-t/\tau)^{\alpha}]) d\tau \quad (30) $$

which gives at long times

$$ \lim_{t \to +\infty} \mathcal{I}^{\alpha} \tilde{\mathbf{L}}(t) = \tilde{\mathbf{L}}_{e_r} \frac{t^{\alpha}}{\Gamma(\alpha + 1)} + \int_{0}^{+\infty} \tau^{\alpha} \tilde{\mathbf{G}}(\tau) d\tau. \quad (31) $$

For classical viscoelasticity, in the case of a loading test with a constant macroscopic strain rate $\dot{\tilde{\varepsilon}}_0$, it follows that

$$ \lim_{t \to +\infty} \bar{\sigma}(t) = \left[ \tilde{\mathbf{L}}_{e_r} t + \int_{0}^{+\infty} \tau \tilde{\mathbf{G}}(\tau) d\tau \right] : \dot{\tilde{\varepsilon}}_0. \quad (32) $$

Relations (29) and (32) provide the physical meaning of the two time integrals of the relaxation spectrum $\int_{0}^{+\infty} \tau^{-1} \tilde{\mathbf{G}}(\tau) d\tau$ and $\int_{0}^{+\infty} \tau \tilde{\mathbf{G}}(\tau) d\tau$ [19].
### 3 Approximate viscoelastic homogenization model

#### 3.1 Approximation of the effective relaxation spectrum $\tilde{G}$

To derive an estimate of the overall viscoelastic relaxation function, it is necessary to approximate the unknown relaxation spectrum $\tilde{G}(\tau)$. The most usual approaches make use of box functions, wedge functions or a sum of Dirac delta functions [60, 14]. The latter case corresponds to a line spectrum

$$ \tilde{G}(\tau) = \sum_{k=1}^{K} \overline{G}_k \delta(\tau - \bar{\tau}_k) \qquad (33) $$

which results in the relaxation function

$$ \tilde{L}_{\alpha}(t) = \tilde{L}_{e_r} + \sum_{k=1}^{K} \overline{G}_k E_{\alpha} [-(t/\bar{\tau}_k)^{\alpha}]. \qquad (34) $$

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 constituted 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

$$ \left\{ \begin{aligned} \tilde{\mathbf{L}}'_{\alpha}(\omega^{\alpha}) &= \tilde{\mathbf{L}}_{e_r} + \sum_{k=1}^{K} \frac{1}{q_k} \left[ (\omega\bar{\tau}_k)^{\alpha} \cos\left(\frac{\pi\alpha}{2}\right) + (\omega\bar{\tau}_k)^{2\alpha} \right] \overline{G}_k, \\ \tilde{\mathbf{L}}''_{\alpha}(\omega^{\alpha}) &= \sum_{k=1}^{K} \frac{1}{q_k} \left[ (\omega\bar{\tau}_k)^{\alpha} \sin\left(\frac{\pi\alpha}{2}\right) \right] \overline{G}_k, \end{aligned} \right. \qquad (35) $$

with $q_k = 1 + 2(\omega\bar{\tau}_k)^\alpha \cos(\frac{\pi\alpha}{2}) + (\omega\bar{\tau}_k)^{2\alpha}$. From the approximation (33) and the conditions on time integrals of the effective relaxation spectrum (25), it follows that

$$ \left\{ \begin{aligned} \sum_{k=1}^{K} \overline{G}_k &= \tilde{\mathbf{L}}_{e_g} - \tilde{\mathbf{L}}_{e_r}, \\ \sum_{k=1}^{K} \overline{\tau}_k^\alpha \overline{\varepsilon} : \overline{G}_k : \overline{\varepsilon} &= \sum_s c^{(s)} \mathbf{L}_f^{(s)} : \langle \varepsilon_r \otimes \varepsilon_r \rangle^{(s)}, \\ \sum_{k=1}^{K} \frac{1}{\overline{\tau}_k^\alpha} \overline{\varepsilon} : \overline{G}_k : \overline{\varepsilon} &= \sum_s c_s \left( \mathbf{G}^{(s)} : \mathbf{L}_f^{(s)-1} : \mathbf{G}^{(s)} \right) : \langle \varepsilon_g \otimes \varepsilon_g \rangle^{(s)}. \end{aligned} \right. \qquad (36) $$

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.
## 3.2 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

$$ \tilde{\mathbf{G}}(\tau) = \sum_{k=1}^{K} 3\bar{\kappa}_k \delta(\tau - \bar{\tau}_k^\kappa)\mathbf{J} + 2\bar{\mu}_k \delta(\tau - \bar{\tau}_k^\mu)\mathbf{K} \quad (37) $$

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

$$ \mathbf{L}_{\alpha}^{(s)}(t) = 3\kappa_{e}^{(s)}\mathbf{J} + 2(\mu_{e_r}^{(s)} + (\mu_{e_g}^{(s)} - \mu_{e_r}^{(s)})E_{\alpha}[-(t/\tau^{(s)})^{\alpha}])\mathbf{K}. \quad (38) $$

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 ($\bar{\kappa}_k = 0$), whence the approximate form of the effective relaxation function

$$ \tilde{\mathbf{L}}_{\alpha}(t) \approx 3\tilde{\kappa}_e \mathbf{J} + 2 \left( \tilde{\mu}_{e_r} + \sum_{k=1}^{K} \bar{\mu}_k E_{\alpha} \left[ -(t/\bar{\tau}_k^{\mu})^{\alpha} \right] \right) \mathbf{K}. \quad (39) $$

Since only shear relaxation times are now involved, we will skip the upper index $\mu$ and use the notation $\bar{\tau}_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

$$ \langle \varepsilon_{eq}^2 \rangle^{(s)} = \frac{1}{c_s} \left[ \frac{1}{3} \frac{\partial \tilde{\kappa}}{\partial \mu^{(s)}} \mathrm{tr}(\tilde{\varepsilon})^2 + \frac{\partial \tilde{\mu}}{\partial \mu^{(s)}} \tilde{\varepsilon}_{eq}^2 \right] \quad \text{and} \quad \langle \mathrm{tr}(\varepsilon)^2 \rangle^{(s)} = \frac{1}{c_s} \left[ \frac{\partial \tilde{\kappa}}{\partial \kappa^{(s)}} \mathrm{tr}(\tilde{\varepsilon})^2 + 3 \frac{\partial \tilde{\mu}}{\partial \kappa^{(s)}} \tilde{\varepsilon}_{eq}^2 \right] \quad (40) $$

with $\varepsilon_{eq}^2 = \frac{2}{3} \mathbf{K} :: (\varepsilon \otimes \varepsilon)$ and $\mathrm{tr}(\varepsilon)^2 = 3 \mathbf{J} :: (\varepsilon \otimes \varepsilon)$. 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) are reduced to

$$ 
\left\{
\begin{aligned}
\sum_{k=1}^{K} \bar{\mu}_k &= \tilde{\mu}_{e_g} - \tilde{\mu}_{e_r}, \\
\sum_{k=1}^{K} \bar{\mu}_k \bar{\tau}_k^{\alpha} &= \sum_{s=1}^{N} (\mu_{e_g}^{(s)} - \mu_{e_r}^{(s)}) \tau^{(s)} \alpha \frac{\partial \tilde{\mu}_{e_r}}{\partial \mu_{e_r}^{(s)}}, \\
\sum_{k=1}^{K} \frac{\bar{\mu}_k}{\bar{\tau}_k^{\alpha}} &= \sum_{s=1}^{N} \frac{(\mu_{e_g}^{(s)} - \mu_{e_r}^{(s)})}{\tau^{(s)} \alpha} \frac{\partial \tilde{\mu}_{e_g}}{\partial \mu_{e_g}^{(s)}}.
\end{aligned}
\right.
\quad (41) $$

The minimal number of terms required in the generalized Prony serie to fulfill the three relations (41) is thus $K = 2$. To determine the four parameters $\bar{\mu}_1, \bar{\mu}_2, \bar{\tau}_1^\alpha$ and $\bar{\tau}_2^\alpha$, 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 [58,6]. This minimal approximation is assessed in the sequel for particulate two-phase composites.

# 4 Applications to particle reinforced two-phase composites

## 4.1 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 equations (41) becomes

$$ \bar{\mu}_1 + \bar{\mu}_2 = A, \quad \bar{\mu}_1 \bar{\tau}_1^\alpha + \bar{\mu}_2 \bar{\tau}_2^\alpha = B \quad \text{and} \quad \frac{\bar{\mu}_1}{\bar{\tau}_1^\alpha} + \frac{\bar{\mu}_2}{\bar{\tau}_2^\alpha} = C \quad (42) $$

with

$$ A = \tilde{\mu}_{e_g} - \tilde{\mu}_{e_r}, \quad B = (\mu_{e_g}^{(1)} - \mu_{e_r}^{(1)}) \tau^{(1)} \alpha \frac{\partial \tilde{\mu}_{e_r}}{\partial \mu_{e_r}^{(1)}} \quad \text{and} \quad C = \frac{(\mu_{e_g}^{(1)} - \mu_{e_r}^{(1)})}{\tau^{(1)} \alpha} \frac{\partial \tilde{\mu}_{e_g}}{\partial \mu_{e_g}^{(1)}}. \quad (43) $$

In order to solve the system (42), we set a fractional relaxation time of the serie equal to the one of the matrix phase: $\bar{\tau}_1^\alpha = \tau^{(1)} \alpha$. The three other parameters are then given by

$$ \bar{\tau}_2^\alpha = \frac{A \bar{\tau}_1^\alpha - B}{C \bar{\tau}_1^\alpha - A}, \quad \bar{\mu}_1 = \frac{B - A \bar{\tau}_1^\alpha}{\bar{\tau}_2^\alpha - \bar{\tau}_1^\alpha} \quad \text{and} \quad \bar{\mu}_2 = A - \bar{\mu}_1. \quad (44) $$
## 4.2 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

$$ \tilde{\mu} = \mu^{(1)} + c_2 \frac{\mu^{(2)} - \mu^{(1)}}{1 + c_1 \frac{\mu^{(2)} - \mu^{(1)}}{\mu^{*} + \mu^{(1)}}} \quad \text{and} \quad \tilde{\kappa} = \kappa^{(1)} + c_2 \frac{\kappa^{(2)} - \kappa^{(1)}}{1 + c_1 \frac{\kappa^{(2)} - \kappa^{(1)}}{\kappa^{*} + \kappa^{(1)}}} \quad (45) $$

with $\mu^*$ and $\kappa^*$ the shear and bulk moduli of the Hill constraint tensor

$$ \mu^* = \frac{\mu^0 (9\kappa^0 + 8\mu^0)}{6(\kappa^0 + 2\mu^0)} \quad \text{and} \quad \kappa^* = \frac{4}{3}\mu^0. \quad (46) $$

The choice ($\mu^0 = \mu^{(1)}$, $\kappa^0 = \kappa^{(1)}$) corresponds to the lower Hashin-Shtrikman bound which coincides with the Mori-Tanaka (MT) model [42] while the choice ($\mu^0 = \tilde{\mu}$, $\kappa^0 = \tilde{\kappa}$) 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].

## 4.3 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 [44, 43, 17].
We consider a symbolic heterogeneous elastic problem (6), with periodic boundary conditions, which corresponds to the stationary harmonic regime at angular frequency $\omega$. By introducing a homogeneous reference material $L_0$, the constitutive law can be rewritten as

$$ \sigma^*(x, i\omega) = L_0 : \varepsilon^*(x, i\omega) + \tau^*(x, i\omega) \quad \text{with} \quad \tau^*(x, i\omega) = (L_0^*(x, i\omega) - L_0) : \varepsilon(x, i\omega) \quad (47) $$

with $\tau^*$ the polarization field. The complex strain solution field thus reads, in real and Fourier spaces,

$$ \left\{ \begin{aligned} \varepsilon^*(x, i\omega) &= \bar{\varepsilon} - [\Gamma_0 * \tau^*](x, i\omega) \\ \hat{\varepsilon}^*(\xi, i\omega) &= -\hat{\Gamma}_0 : \hat{\tau}^*(\xi, i\omega), \quad \forall \xi \neq 0, \quad \hat{\varepsilon}(0, i\omega) = \bar{\varepsilon} \end{aligned} \right. \quad (48) $$

with $\xi$ the wave vector. $\hat{\Gamma}_0$ represents the Fourier transform of the strain Green operator associated with the reference medium $L_0$

$$ \hat{\Gamma}_0 = [\xi \otimes \kappa^{-1} \otimes \xi]_s \quad (49) $$

with the acoustic tensor $\kappa = \xi \cdot L^0 \cdot \xi$. The notation $[.]_s$ indicates minor and major symmetrizations. To solve the Lippman-Schwinger equation for complex strain field $\varepsilon^*$ (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

$$ \mu_0 = \sqrt{\mathrm{Re}(\mu^{(1)}) \mathrm{Re}(\mu^{(2)})} \quad \text{and} \quad \kappa_0 = \sqrt{\mathrm{Re}(\kappa^{(1)}) \mathrm{Re}(\kappa^{(1)})} \quad (50) $$

## 4.4 Overall behaviour of viscoelastic composites

### 4.4.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 ($\omega \to +\infty$), to 3000, in the relaxed state ($\omega \to 0$). Distributions of particles in a cubic unit-cell, with geometric periodic conditions, are obtained by using the Random Sequential
<table><caption>Table 1: Constitutive parameters of reinforced particulate composites with a fractional Zener viscoelastic matrix (38). Bulk and shear elastic moduli (κ, μ) are expressed in GPa.</caption><thead><tr><th></th><th>κ<sub>e</sub></th><th>μ<sub>e<sub>r</sub></sub></th><th>μ<sub>e<sub>g</sub></sub></th><th>τ</th><th>α</th></tr></thead><tbody><tr><td>Viscoelastic matrix (1)</td><td>4</td><td>0.01</td><td>1</td><td>2.65 × 10<sup>-2</sup> s</td><td>0.7</td></tr><tr><td>Elastic inclusions (2)</td><td>40</td><td>30</td><td>30</td><td>—</td><td>—</td></tr></tbody></table>

Adsorption (RSA) algorithm [62,59] (Appendix C). Full-field FFT simulations, for a prescribed angular frequency ω, are performed on discretized unit-cells with a regular grid of 255³ 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 $\tilde{\mu}_{\alpha}^*(i\omega) = \tilde{\mu}_{\alpha}'(\omega) + i\tilde{\mu}_{\alpha}''(\omega)$ 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 $\tilde{\mu}_{\alpha}'(\omega)$ in the range $\omega \in [10^{-2}; 10^{-4}]$ are reported in Figure 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.

### 4.4.2 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 frequency. The approximation of the real and imaginary parts of the effective complex shear moduli $\tilde{\mu}_{\alpha}^*(i\omega)$ reads

$$ \left\{ \begin{aligned} \operatorname{Re}(\tilde{\mu}_{\alpha}^*(i\omega)) &= \tilde{\mu}_{\alpha}'(\omega^{\alpha}) = \tilde{\mu}_{e_r} + \sum_{k=1}^{2} \frac{1}{q_k} \left[ (\omega \bar{\tau}_k)^{\alpha} \cos\left(\frac{\pi \alpha}{2}\right) + (\omega \bar{\tau}_k)^{2\alpha} \right] \bar{\mu}_k, \\ \operatorname{Im}(\tilde{\mu}_{\alpha}^*(i\omega)) &= \tilde{\mu}_{\alpha}''(\omega^{\alpha}) = \sum_{k=1}^{2} \frac{1}{q_k} \left[ (\omega \bar{\tau}_k)^{\alpha} \sin\left(\frac{\pi \alpha}{2}\right) \right] \bar{\mu}_k. \end{aligned} \right. \qquad (51) $$
<figure><img src="image_3.png" /><figcaption>Fig. 1: Mean-field estimates and FFT computations of the effective storage shear moduli $\tilde{\mu}'_\alpha(\omega)$ of particulate composites with different volume fractions $c$ (a): 0.1, (b): 0.3 and (c): 0.5. MT: Mori-Tanaka model and (G)SC: (Generalized) Self-Consistent schemes.</figcaption></figure>

with $q_k = 1 + 2(\omega\bar{\tau}_k)^\alpha \cos\left(\frac{\pi\alpha}{2}\right) + (\omega\bar{\tau}_k)^{2\alpha}$. 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 $\tilde{\mu}_{e_r}$ and $\tilde{\mu}_{e_g}$. The comparison with reference FFT computations shows an overall good agreement of the approximate GSC model (Figure 2). It is worth mentioning that the exact shear relaxation spectrum of the GSC model presents a continuous part [2].

## 5 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
<figure><img src="image_4.png" /><figcaption>Fig. 2: Effective storage moduli $\tilde{\mu}'_{\alpha}$ (a) and loss factor $\tan \delta = \tilde{\mu}''_{\alpha}/\tilde{\mu}'_{\alpha}$ (b) of reinforced viscoelastic composites for particle volume fractions $c = 0.1, 0.3$ and $0.5$. Comparisons between FFT full-field simulations and approximate model (51).</figcaption></figure>

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

$$ (f \odot g)(t) = \frac{d}{dt}(f * g)(t) = \int_{0}^{t} f(t-u)\dot{g}(t) du \qquad (A1) $$
If $g$ is only piecewise continuous and differentiable, its time derivative in (A1) contains Dirac masses at discontinuity points $t_n$ and the Stieljes convolution product must be understood as

$$ (f \otimes g)(t) = \int_{0}^{t} f(t-u) \dot{g}(u) du + \sum_{n} f(t-t_n)[g]_{n} \quad (A2) $$

where $[g]_n$ is the discontinuity of $g$ at time $t_n$ and $\dot{g}(u)$ is the usual derivative of $g$.

The Laplace-Carson (LC) transform of a function $f(t)$ is defined by

$$ \mathcal{L}\mathcal{C}(f(t)) = f^*(p) = p \int_{0}^{+\infty} e^{-pt} f(t) dt \quad (A3) $$

From (A1) and (A3), it follows that

$$ (f \otimes g)^*(p) = f^*(p)g^*(p). \quad (A4) $$

## B Useful relations for fractional calculus

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

### B.1 Gamma function

The Gamma function $\Gamma$ is a continuation of the factorial function to complex numbers. It is defined by

$$ \Gamma(z) = \int_{0}^{+\infty} u^{z-1} e^{-u} du, \quad \forall z \in \mathbb{C} \setminus \{0\} \cup \mathbb{Z}^{-} \quad (B1) $$

Integration by parts leads to the following property

$$ \Gamma(z+1) = z\Gamma(z) \quad (B2) $$

It is also noted that $\Gamma(1) = 1$ and $\Gamma(n) = (n-1)!, \forall n \in \mathbb{N}$.

### B.2 Fractional integral and derivative operators

Cauchy's formula for successive integrations of a causal function $f$, with null initial conditions, reads

$$ \mathcal{I}^n f(t) = \frac{1}{(n-1)!} \int_0^t (t-u)^{n-1} f(u) du, \quad \forall n \in \mathbb{N}, \quad (B3) $$

namely

$$ \mathcal{I}^n f(t) = (Y_n * f)(t) \quad \text{with} \quad Y_n(t) = \frac{t^{n-1}}{(n-1)!}, \quad \forall n \in \mathbb{N}. \quad (B4) $$

By using the function $\Gamma$, the integral operator $\mathcal{I}$ can be extended in a straightforward manner to non-integer order $\alpha > 0$

$$ \mathcal{I}^{\alpha} f(t) = (Y_{\alpha} * f)(t) \quad \text{with} \quad Y_{\alpha}(t) = \frac{t^{\alpha-1}}{\Gamma(\alpha)}, \quad \alpha > 0. \quad (B5) $$
It can be also noted that

$$ \mathcal{I}^{\alpha} f(t) = \mathcal{I}^{\alpha+1} \dot{f}(t) = (Y_{\alpha+1} * \dot{f})(t) = (Y_{\alpha+1} \circledast f)(t) \quad (B6) $$

The non-integer derivative operator $\mathcal{D}^{\alpha}$ is defined by successive derivation (of order 1) and integration of order $1-\alpha$, that is

$$ \mathcal{D}^{\alpha} f(t) = (Y_{1-\alpha} * \dot{f})(t) = (Y_{1-\alpha} \circledast f)(t). \quad (B7) $$

From (A4), the LC transform of the fractional derivative thus simply reads

$$ \mathcal{L}\mathcal{C}(\mathcal{D}^{\alpha} f(t)) = (\mathcal{D}^{\alpha} f)^*(p) = Y_{1-\alpha}^*(p) f^*(p). \quad (B8) $$

By noting that $\mathcal{L}\mathcal{C}(t^{-\alpha}) = \Gamma(1-\alpha) p^{\alpha}$, it gives

$$ (\mathcal{D}^{\alpha} f)^*(p) = p^{\alpha} f^*(p), \quad \alpha > 0. \quad (B9) $$

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_{\alpha}$ is defined by the following power series

$$ E_{\alpha}(z) = \sum_{n=0}^{+\infty} \frac{z^n}{\Gamma(1 + \alpha n)}, \quad \forall z \in \mathbb{C}, \quad 0 < \alpha \le 1. \quad (B10) $$

Alternatively, it admits an integral representation [24,37] which reads

$$ E_{\alpha}(-t^{\alpha}) = \int_{0}^{+\infty} H_{\alpha}(\theta) e^{-t/\theta} d\theta, \quad t > 0, \quad 0 < \alpha < 1. \quad (B11) $$

with

$$ H_{\alpha}(\theta) = \frac{\theta^{\alpha-1} \sin(\alpha\pi)}{\pi(1+2\theta^{\alpha}\cos(\alpha\pi)+\theta^{2\alpha})}. \quad (B12) $$

It can be also noted that $\int_{0}^{+\infty} H_{\alpha}(\theta) d\theta = 1$ since $E_{\alpha}(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

$$ \mathcal{D}^{\alpha} E_{\alpha}(t^{\alpha}) = E_{\alpha}(t^{\alpha}). \quad (B13) $$

The function $E_{\alpha}$ naturally appears in fractional calculus since it is solution of the fractional differential equation

$$ \mathcal{D}^{\alpha} f(t) + \lambda^{\alpha} f(t) = 0, \quad \forall t > 0 \quad \text{with} \quad f(0) = 1. \quad (B14) $$

The solution reads [13] $f(t) = E_{\alpha}(-\lambda t)^{\alpha}$. The asymptotic expansion of the function $E_{\alpha}(-t^{\alpha})$ at short and long times reads

$$ \lim_{t \to 0} E_{\alpha}(-t^{\alpha}) = 1 - \frac{t^{\alpha}}{\Gamma(1+\alpha)} \quad \text{and} \quad \lim_{t \to +\infty} E_{\alpha}(-t^{\alpha}) = \frac{t^{-\alpha}}{\Gamma(1-\alpha)} \quad (\alpha \neq 1). \quad (B15) $$

The Laplace-Carson transform of the Mittag-Leffler function satisfies

$$ \mathcal{L}\mathcal{C}(E_{\alpha} [(-t/\tau)^{\alpha}]) = \frac{(p\tau)^{\alpha}}{1+(p\tau)^{\alpha}}. \quad (B16) $$
## B.4 Fractional dashpot constitutive relations

The constitutive relation of a 1D fractional element reads [7]

$$ \sigma(t) = E\tau^{\alpha}\mathcal{D}^{\alpha}\varepsilon(t) \qquad (B17) $$

with elastic modulus E and fractional relaxation time $\tau^{\alpha}$. From the definition (B7), it follows that

$$ \sigma(t) = \int_{0}^{t} R(t-u)\dot{\varepsilon}(u) du \qquad (B18) $$

where the relaxation function $R(t)$ reads

$$ R(t) = E Y_{1-\alpha} \left(\frac{t}{\tau}\right) = \frac{E}{\Gamma(1-\alpha)} \left(\frac{\tau}{t}\right)^{\alpha}. \qquad (B19) $$

The inversion of the constitutive law (B17) simply yields

$$ \varepsilon(t) = \frac{1}{E\tau^{\alpha}}\mathcal{I}^{\alpha}\sigma(t) = \frac{1}{E\tau^{\alpha}}\mathcal{I}^{\alpha+1}\dot{\sigma}(t). \qquad (B20) $$

From relation (B6), it may be written as

$$ \varepsilon(t) = \int_{0}^{t} F(t-u)\dot{\sigma}(u) du \qquad (B21) $$

where the creep (or retardation) function $F(t)$ reads

$$ F(t) = \frac{1}{E} Y_{1+\alpha} \left(\frac{t}{\tau}\right) = \frac{1}{E \Gamma(1+\alpha)} \left(\frac{t}{\tau}\right)^{\alpha}. \qquad (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 nonoverlapping geometric objects into a fixed volume. The size distribution of the spherical inclusions follows a lognormal density function $\phi(r)$ with r the radius of the particles

$$ \phi(r) = \frac{1}{r \sigma \sqrt{2\pi}} e^{-(\ln(r)-\mu)^2/2\sigma^2}, \quad \sigma > 0, \qquad (C1) $$

where the parameters $\mu$ and $\sigma$ 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 $\mathcal{A}(+\infty) = \int_{0}^{+\infty} \phi(r; \mu, \sigma) dr = 1$. In practice, it is thus necessary to choose a maximal radius $r_{max}$ corresponding to a given value $\mathcal{A}(r_{max})$. We have chosen a maximal radius corresponding to $\mathcal{A}(r_{max}) = 0.99$

$$ r_{max} = e^{\sigma\psi+\mu} \quad \text{with} \quad \psi = \sqrt{2\pi} \operatorname{erf}^{-1} \left( \frac{22}{25} \right) \qquad (C2) $$
<figure><img src="image_5.png" /><figcaption>Fig. C.1: Polydisperse microstructure composed of spherical inclusions with volume fraction $c = 0.3$. Parameters: $\mu = 0$, $\sigma = 1$ and $P = 100$.</figcaption></figure>

By considering evenly spaced radii, the knowledge of the limit radius allows to determine each radius associated to 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 to an inclusion size family $i$ read

$$ \left\{ \begin{array}{l} r_i = \frac{r_{\max}}{2n} (2i - 1), \quad n > 2, \\ P_i = A_i P \quad \text{with} \quad A_i = \int_{r_i^-}^{r_i^+} \phi(r) \, dr, \quad r_i^{\pm} = r_i \pm \frac{r_{\max}}{2n}. \end{array} \right. \qquad (C3) $$

References

1. Bagley, R.L., Torvik, P.J.: A theoretical basis for the application of fractional calculus to viscoelasticity. J. Rheol. 27, 201–210 (1983)

2. Beurthey, S., Zaoui, A.: Structural morphology and relaxation spectra of viscoelastic heterogeneous materials. Eur. J. Mech. A/Solids 19, 1–16 (2000)

3. Blair, G.W.S., Caffyn, J.: The classification of the rheological properties of industrial materials in the light of power-law relations between stress, strain and time. Journal of Scientific Instruments 19(6), 88 (1942)

4. Bobeth, M., Diener, G.: Field fluctuations in multicomponent mixtures. J. Mech. Phys. Solids 34, 1–17 (1986)

5. Brenner, R., Masson, R.: Improved affine estimates for nonlinear viscoelastic composites. Eur. J. Mech. A/Solids 24, 1002–1015 (2005)

6. Brenner, R., Suquet, P.: Overall response of viscoelastic composites and polycristals: exact asymptotic relations and approximate estimates. Int. J. Solids Struct. 50, 1824–1838 (2013)

7. Caputo, M., Mainardi, F.: Linear models of dissipation in anelastic solids. Rivista del nuovo cimento 1, 161–198 (1971)
8. Caputo, M., Mainardi, F.: A new dissipation model based on memory mechanism. *Pure and Applied Geophysics* **91**, 134–147 (1971)

9. Cherkaev, A.V., Gibiansky, L.V.: Variational principles for complex conductivity, viscoelasticity, and similar problems in media with complex moduli. *J. Math. Phys.* **35**, 127–145 (1994)

10. Christensen, R.M.: A critical evaluation for a class of micromechanics models. *J. Mech. Phys. Solids* **38**, 379–404 (1990)

11. Christensen, R.M., Lo, K.H.: Solutions for effective shear properties in three phase sphere and cylinder models. *J. Mech. Phys. Solids* **27**, 315–330 (1979)

12. Dinzart, F., Lipinski, P.: Self-consistent approach of the constitutive law of a two-phase visco-elastic material described by fractional derivative models. *Arch. Mech.* **62**, 135–156 (2010)

13. Dubois, F., Galucio, A.C., Point, N.: Introduction à la dérivation fractionnaire. *Théorie et applications. Techniques de l'Ingenieur* (2010)

14. Eyre, D., Milton, G.W., Lakes, R.S.: Bounds for interpolating complex effective moduli of viscoelastic materials from measured data. *Rheol. Acta* **41**, 461–470 (2002)

15. Eyre, D.J., Milton, G.W.: A fast numerical scheme for computing the response of composites using grid refinement. *Journal of Physique III* **6**, 41–47 (1999)

16. Ferry, J.D.: Viscoelastic properties of polymers. *John Wiley and Sons* (1980)

17. Figliuzzi, B., Jeulin, D., Faessel, M., Willot, F., Koishi, M., Kowatari, N.: Modelling the microstructure and the viscoelastic behaviour of carbon black filled rubber materials from 3d simulations. *Technische Mechanik* **32**(1-2), 22–46 (2016)

18. Francfort, G., Leguillon, D., Suquet, P.: Homogénéisation de milieux viscoélastiques linéaires de Kelvin-Voigt. *C.R. Acad. Sci. Paris* **296**, 287–290 (1983)

19. Gallican, V., Brenner, R., Suquet, P.: Exact asymptotic relations for the effective response of linear viscoelastic heterogeneous media. *C.R. Mecanique* **345**, 742–751 (2017)

20. Gemant, A.: A method of analyzing experimental results obtained from elasto-viscous bodies. *Physics* **7**, 311–317 (1936)

21. Gemant, A.: On fractional differentials. *Phil. Mag.* **25**, 540–549 (1938)

22. Ghossein, E., Lévesque, M.: A fully automated numerical tool for a comprehensive validation of homogenization models and its application to spherical particles reinforced composites. *Int. J. Solids Struct.* **49**, 1387–1398 (2012)

23. Gibiansky, L.V., Milton, G.W.: On the effective viscoelastic moduli of two-phase media. *I. Rigorous bounds on the complex bulk modulus*. *Proc. R. Soc. Lond. A440*, 163–188 (1993)

24. Gorenflo, R., Mainardi, F.: Fractals and fractional calculus in continuum mechanics (CISM courses and lectures, vol. 378), chap. "Fractional calculus: integral and differential equations of fractional order", pp. 223–276. Springer-Verlag Wien GmbH (1997)

25. Gusev, A.A.: Time domain finite element estimates of dynamic stiffness of viscoelastic composites with stiff spherical inclusions. *Int. J. Solids Struct.* **88-89**, 79–87 (2016)

26. Hashin, Z.: Viscoelastic behavior of heterogeneous media. *J. Appl. Mech.* **32**, 630–636 (1965)
27. Hashin, Z.: Complex moduli of viscoelastic composites - I. general theory and application to particulate composites. Int. J. Solids Struct. **6**, 539–552 (1970)

28. Hervé, E., Zaoui, A.: Modelling the effective behavior of nonlinear matrix-inclusion composites. Eur. J. Mech. A/Solids **9**, 505–515 (1990)

29. Heymans, N., Bauwens, J.C.: Fractal rheological models and fractional differential equations for viscoelastic behavior. Rheol. Acta **33**, 210–219 (1994)

30. Hill, R.: A self-consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solids **13**(4), 213 – 222 (1965)

31. Koeller, R.C.: Applications of fractional calculus to the theory of viscoelasticity. J. Appl. Mech. **51**, 299–307 (1984)

32. Koeller, R.C.: Polynomial operators, Stieltjes convolution, and fractional calculus in hereditary mechanics. Acta Mech. **58**, 251–264 (1986)

33. Kreher, W.: Residual stresses and stored elastic energy of composites and polycrystals. Journal of the Mechanics and Physics of Solids **38**, 115 – 128 (1990)

34. Laws, N., McLaughlin, R.: Self-consistent estimates for the viscoelastic creep compliances of composite materials. Proc. R. Soc. Lond. **A359**, 251–273 (1978)

35. Lion, A.: On the thermodynamics of fractional damping elements. Continuum Mech. Thermodyn. **9**, 83–96 (1997)

36. Mainardi, F.: Fractional Calculus and Waves in Linear Viscoelasticity. Imperial College Press (2010)

37. Mainardi, F.: On some properties of the Mittag-Leffler function $E_\alpha(-t^\alpha)$, completely monotone for $t > 0$ with $0 < \alpha < 1$. Discrete. Cont. Dyn.-B **19**, 2267–2278 (2014)

38. Mandel, J.: Cours de mécanique des milieux continus. Gauthier-Villars Editeur, Paris, France (1966)

39. Masson, R., Brenner, R., Castelnau, O.: Incremental homogenization approach for ageing viscoelastic polycrystals. C. R. Mecanique **340**, 378–386 (2012)

40. Masson, R., Zaoui, A.: Self-consistent estimates for the rate-dependent elastoplastic behaviour of polycrystalline materials. J. Mech. Phys. Solids **47**, 1543–1568 (1999)

41. Milton, G.W., Berryman, J.G.: On the effective viscoelastic moduli of two-phase media. II. Rigorous bounds on the complex shear modulus in three dimensions. Proc. R. Soc. Lond. A **453**, 1849–1880 (1997)

42. Mori, T., Tanaka, K.: Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metallurgica **21**(5), 571 – 574 (1973)

43. Moulinec, H., Silva, F.: Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials. Int. J. Num. Meth. Engng **97**, 960–985 (2014)

44. Moulinec, H., Suquet, P.: A numerical method for computing the overall response of nonlinear composites with complex microstructures. Comput. Methods Appl. Mech. Engrg. **157**, 69–94 (1998)

45. Nutting, P.G.: A new general law of deformation. J. Franklin Inst. **191**, 679–685 (1921)

46. Oldham, K.B., Spanier, J.: The fractional calculus. Theory and applications of differentiation and integration to arbitrary order. Academic Press (1974)
47. Podlubny, I.: Fractional differential equations. Academic Press (1999)

48. Ponte Castañeda, P., Suquet, P.: Nonlinear composites. Adv. Appl. Mech. **34**, 171–302 (1998)

49. Rekik, A., Brenner, R.: Optimization of the collocation inversion method for the linear viscoelastic homogenization. Mech. Res. Comm. **38**, 305–308 (2011)

50. Remillat, C.: Damping mechanism of polymers filled with elastic particles. Mech. Mater. **39**, 525–537 (2007)

51. Ricaud, J.M., Masson, R.: Effective properties of linear viscoelastic heterogeneous media: Internal variables formulation and extension to ageing behaviours. Int. J. Solids Struct. **46**, 1599–1606 (2009)

52. Sanchez-Hubert, J., Sanchez-Palencia, E.: Sur certains problèmes physiques d'homogénéisation donnant lieu à des phénomènes de relaxation. C. R. Acad. Sci. Paris Série A **286**, 903–906 (1978)

53. Schapery, R.A.: Viscoelastic behavior and analysis of composite materials. In: G.P. Sendeckyj (ed.) Mechanics of composite materials, vol. 2, pp. 85–167. New York: Academic Press (1974)

54. Schapery, R.A.: Linear elastic and viscoelastic deformation behaviour of ice. J. Cold Reg. Eng. **11**, 271–290 (1997)

55. Schiessel, H., Blumen, A.: Hierarchical analogues to fractional relaxation equations. J. Phys. A: Math. Gen. **26**, 5057–5069 (1993)

56. Sevostianov, I., Levin, V., Radi, E.: Effective properties of linear viscoelastic microcracked materials: Application of Maxwell homogenization scheme. Mech. Mater. **84**, 28–43 (2015)

57. Suquet, P.: Homogenization Techniques for Composite Media (Lecture notes in Physics, vol. 272), chap. “Elements of homogenization for inelastic solid mechanics”, pp. 194–278. Springer-Verlag (1987)

58. Suquet, P.: Four exact relations for the effective relaxation function of linear viscoelastic composites. Comptes Rendus Mecanique **340**(4), 387 – 399 (2012). Recent Advances in Micromechanics of Materials

59. Torquato, S.: Random heterogeneous materials. Springer (2002)

60. Tschoegl, N.: The phenomenological theory of linear viscoelastic behavior – An introduction. Springer-Verlag (1989)

61. Turner, P.A., Tomé, C.N.: Self-consistent modeling of visco-elastic polycrystals : application to irradiation creep and growth. J. Mech. Phys. Solids **41**(7), 1191–1211 (1993)

62. Widom, B.: Random sequential addition of hard spheres to a volume. J. Chem. Phys. **44**, 3888–3894 (1966)

63. Willis, J.R.: Variational and related methods for the overall properties of composites. Adv. Appl. Mech. **21**, 2–78 (1981)

64. Yin, D., Wu, H., Cheng, C., Chen, Y.Q.: Fractional order constitutive model of geomaterials under the condition of triaxial test. Int. J. Numer. Anal. Met. **37**, 961–972 (2013)