Does a Fractal Microstructure Require a Fractional Viscoelastic Model?

Abstract

The question addressed by this paper is tackled through a continuum micromechanics model of a 2D random checkerboard, in which one phase is linear elastic and another linear viscoelastic of integer-order. The spatial homogeneity and ergodicity of the material statistics justify homogenization in the vein of the Hill–Mandel condition for viscoelastic media. Thus, uniform kinematic- or traction-controlled boundary conditions, applied to sufficiently large domains, provide macroscopic (RVE level) responses. With computational mechanics, this strategy is applied over the entire range of the relative content of both phases. Setting the volume fraction of either the elastic phase or the viscoelastic phase at the critical value (≃0.59) results in fractal patterns of site-percolation. Extensive simulations of boundary value problems show that, for a viscoelastic composite having such a fractal structure, the integer (not fractional) calculus model is adequate. In other words, the spatial randomness of the composite material—even in the fractal regime—is not necessarily the cause of the fractional order viscoelasticity.

1. Introduction

Dating back to when [1,2] was written, fractional calculus (i.e., non-integer order derivatives and integrals) models have been known to provide good fits to constitutive responses of many viscoelastic materials, both natural and man-made (see also [3,4,5]) On the other hand, many natural materials display fractal spatial patterns [6], which has motivated some researchers to hypothesize that fractal structures dictate the fractional response. Thus, a question arises: Are fractional viscoelastic models dictated by fractal (micro)structures of materials? Put differently, do fractal materials made of non-fractional phases require integer-order or fractional-order viscoelastic models? This issue has been studied for fractal arrangements of rheological elements [7,8,9,10], albeit without a 2D or 3D field theory—such as continuum mechanics, homogenization, and associated boundary value problems—truly taken into account.

The question posed above (and alluded to in the title of this paper) is approached in the present study through a micromechanics model of a random viscoelastic microstructure, homogenized according to a Hill–Mandel condition [11,12]. More specifically, we take the microstructure as a two-phase planar random chessboard, in which one phase is elastic and another viscoelastic. By setting the volume fraction of either the elastic phase or the viscoelastic phase at the critical value (≃0.59), we have a Bernoulli percolation, which is known to result in a fractal pattern with fractal dimension D≃1.9 D≃1.9 [13]. At the same time, the spatial homogeneity and ergodicity of the material statistics allow for the application of uniform kinematic- or traction-controlled boundary conditions, which, for sufficiently large domains, provide nearly macroscopic (i.e., almost scale-independent) responses. In other words, this involves the upscaling of a statistical volume element (SVE) (i.e., from a mesoscale level) to a representative volume element (RVE) (i.e., to a macroscale level). Overall, with this methodology and the tools at hand, we study the constitutive response not only at the Bernoulli percolation point but also for the entire range of volume fractions of the viscoelastic phase.

We begin the paper by setting up the random viscoelastic microstructure model, followed by a Hill–Mandel condition for materials and an application of the computational methodology for determination of mesoscale responses over the entire range of volume fractions. Next, classical versus fractional calculus models of viscoelasticity are discussed. Finally, we explore the applicability of a fractional calculus model for our random microstructure at critical points of percolation of the elastic phase and the viscoelastic phase.

2. Viscoelasticity of Random Microstructure

The onset of percolation depends on a particular geometry and distribution of component phases. For instance, in conductive polymer composites, the polymer acts as the insulating phase while the carbonaceous fillers are conductive. While the elastic properties of two-phase percolating systems of springs or beams are known (e.g., Chapter 4 in [11]), an analogous transition in mechanical properties of elastic-viscoelastic composites—i.e., from an elastic to a viscoelastic type as a function of volume fraction—has not yet been investigated. In order to clearly identify an elastic-viscoelastic transition, a significantly large contrast will be mixing a linear viscoelastic phase with a comparatively stiff linear elastic phase. This setup is consistent with the fact that carbon can be much harder than most polymers. Thus, with the volume fraction _{vf vf} signifying the content of the viscoelastic phase, _vf=0% _vf=0% means “all elastic”, while _vf=100% _vf=100%means “all viscoelastic”.

In the sequel, we focus on one simple planar microstructure: a random checkerboard (also called random chessboard). This microstructure is rather easy to simulate with square-shaped finite elements, although large domains needed for homogenization at large scales are computationally challenging. Technically speaking, the random checkerboard is generated from a Bernoulli process on a Cartesian lattice, with each square cell occupied, independently of realizations at all other cells, with probability p and 1−p 1−pby viscoelastic and elastic phases, respectively. Letting L denote the domain size and d the edge length of each square cell, we introduce a non-dimensional mesoscale parameter:

δ=L/d(d=sizeofonechessboardsite).

Hence, for a square δ.δ δ×δ lattice, the number of different realizations (or the size of the entire sample space Ω Ω ) is Ω=^2δ.δ Ω=^2δ×δ . Each elementary event ω ω gives rise to one realization _Bδω _Bδω of a random medium _Bδ={_Bδω;ω∈Ω} _Bδ={_Bδω;ω∈Ω} . _Bδω _Bδω occurs with probability 1/^2δ.δ 1/^2δ×δ . Effectively, according to the probability p, the volume fraction _{vf vf}of the viscoelastic phase is determined.

Following [11,14,15], a computational mechanics method is employed to determine the response of any specific _Bδω _Bδω . Clearly, the random checkerboard microstructure is spatially statistically homogeneous and ergodic, such that a scale-dependent homogenization according to the Hill–Mandel condition is valid. However, since the finite-size scaling is not of primary interest here, we simulate mesoscale domains at mesoscales that are as large as computationally possible: δ=64 δ=64 for weaker fluctuations, and δ=128 δ=128and 256 for stronger fluctuations. Note that viscoelastic boundary value problems are much more demanding in terms of computing time than the elastic ones.

Assuming the viscoelastic response on microscale (i.e., of each phase) to be isotropic, the time-dependent stiffness tensor C(t) C(t) is characterized by the relaxation shear modulus μt μt and the relaxation bulk modulus k(t) k(t)

C(t)=2μtK+3k(t)J

where J J and K K are the spherical and deviatoric parts of the unit fourth-order tensor. The detailed time-dependence of μ(t) μ(t) and k(t) k(t)can be further modeled numerically via the Prony series

μ(t)=_μ01−∑n=1N_gn¯(1−^e−t_τn),k(t)=_k01−∑m=1M_km¯(1−^e−t_τm)

where _{μ0 μ0} and _{k0 k0} are instantaneous moduli, with _gnk¯ _gnk¯ , _km¯ _km¯ and _{τi τi} being the fitted parameters that describe the material’s behavior at different time scales. Generally speaking, the time dependence of shear and bulk is not necessarily the same. Different sets of ( _gnk¯ _gnk¯ , _km¯ _km¯ , _{τi τi}) result in different frequency-dependent viscoelastic responses.

The elastic phase is set to be one thousand times harder in both shear and bulk compared to the instantaneous modulus of the viscoelastic phase, that is,

_{μE μ0}=_{kE k0}=1,000.

For simplicity, when modeling the viscoelastic phase, only one term of the Prony series is considered. Effects of time constants in both the time and frequency domains are shown in Figure 1. To grasp a rapid and sharp reduction in stiffness under a relaxation test, we choose the case _μ1k¯=_k1¯=0.9,_τ1=0.25 _μ1k¯=_k1¯=0.9,_τ1=0.25, which corresponds to a more rapid and sharp reduction in relaxation tests and a frequency-sensitive viscoelastic phase.

Figure 1. Viscoelastic responses for different time dependent parameters: (a) Normalized relaxation modulus; (b) Normalized complex modulus magnitude and phase angle.

First, consider two boundary value problems for the body _Bδ(ω) _Bδ(ω):

1. kinematic uniform boundary conditions (KUBC):

_ui(x,t)=_ϵij0(t)_xj,∀x∈∂_Bδ

* static uniform boundary conditions (SUBC):

_ti(x,t)=_σij0(t)_nj∀x∈∂_Bδ,

where _σij0(t) _σij0(t) and _ϵij0(t) _ϵij0(t)are the imposed constant tensors. The KUBC and SUBC problems are used, respectively, to conduct relaxation and creep tests in shear.

Applying the loading as a Heaviside (step) function, we determine the apparent properties by taking a volume average of the composite’s response:

1. under KUBC, letting _ui(x,t)=_ϵij0H(t)_xj(x) _ui(x,t)=_ϵij0H(t)_xj(x) on ∂_Bδ(ω) ∂_Bδ(ω),

_μδd(ω,t)=_σ12(ω,t)¯_ϵ120

* under SUBC, letting _ti(x,t)=_σij0H(t)_xj(x) _ti(x,t)=_σij0H(t)_xj(x) on ∂_Bδ(ω) ∂_Bδ(ω),

_Jδt(ω,t)=_ϵ12(ω,t)¯_σ120.

The relaxation stress and creep strain, after stabilization ( t=2 t=2 s), are plotted in Figure 2 for a wide range of volume fractions. When the volume fraction of the viscoelastic phase exceeds ∼60%, the elastic-viscoelastic composite relaxes significantly with much less stress remaining in the body after 2 s under KUBC. Similarly, when the elastic phase occupies more than ∼60%, the composite behaves much more stiffly, and there is a drastic accumulation of specimen strain under SUBC. Both situations show the dominance of the actual percolating phase.

Figure 2. Random checkerboards ( 64×64 ) at different volume fractions: (a) relaxation stress and creep strain after stabilization of response; (b) microstructures ( 128×128 ) showing the (black) viscoelastic phase embedded in the (white) elastic phase.

Since the absolute values of the composite’s moduli vary with different volume fractions, it is convenient to normalize the relaxation modulus by the initial value. In this way, all the normalized relaxation moduli scale between 0 and 1, and indicate the relative reduction in the overall stiffness under relaxation. Since a pure elastic material should exhibit no time dependence and sustain a constant modulus under relaxation, the normalized relaxation modulus can be utilized as an indicator of transition from elasticity to viscoelasticity.

The normalized relaxation moduli for all volume fractions are plotted in Figure 3. It is seen that, when the volume fraction of the elastic phase is less than 40% ( _vf>0.59 ), the composite shows significant softening after relaxation, which is similar to the viscoelastic matrix. However, once the content of elastic suspension is greater than 60% ( _vf<0.41), a reduction in stiffness becomes negligible and elasticity dominates the composite’s response. In this case, the damping properties associated with viscoelasticity are almost locked away and do not take effect.

Figure 3. Normalized relaxation modulus of random checkerboard at different volume fractions.

The critical point for Bernoulli site percolation is at p≃0.59 [16], so the volume fraction from ≃0.41 to ≃0.59 is the range where a strong transition from an elastic-dominated to a viscoelastic-dominated behavior takes place and a competition between both phases is strongest.

To further quantify this elastic-viscoelastic transition, we propose the parameter relative viscoelasticity, which is calculated by dividing the reduction in normalized modulus in the composite by the maximum possible reduction, and that is the reduction in a normalized modulus if only a pure viscoelastic phase exists. Hence, the relative viscoelasticity for the pure elastic phase is 0, while the relative viscoelasticity for the pure viscoelastic phase is 100%.

The relative viscoelasticity is evaluated for the random checkerboard across different volume fraction through Monte Carlo sampling under KUBC. From Figure 4, it can be seen that, for the large contrast of 1000, the mesoscale δ=64 is still not enough for the RVE since the statistical variation exists. At δ=128 , a much sharper change in terms of viscoelasticity occurs between the volume fraction 0.41 and 0.59 . With the viscoelastic phase exceeding 60% , more than 90% of the viscoelasticity of the component phase has been incorporated in the composites. This observation is consistent with the percolation theory, which states that the site percolation on the square lattice occurs at the critical probability _pc≃0.59.

Figure 4. Reduction in normalized modulus after relaxation of random checkerboards at the full range of volume fractions and at two mesoscales: δ=64 and 128.

Note that a sharp change in mechanical response cannot be observed because the elasticity and viscoelasticity are inherently coupled in our model. In the generalized Maxwell model, at least one spring is added, which forces the relaxation stress to be non-vanishing at infinite time. Hence, the constitutive relation of the viscoelastic phase will always be categorized as “solids” as opposed to “fluids”. However, as the relaxation stress after stabilization decreases, the transition slope in Figure 4 tends to become steeper and finally approaches a step shape, which represents the percolation problem of elastic solids and viscous fluids.

3. Fractional Calculus Models in Viscoelasticity 3.1. Background

Considering that the speculations about generalizing the notion of a derivative to non-integer order started with Gottfried W. Leibnitz and Leonhard Euler in the 17th and 18th centuries, the fractional calculus appears to be almost as old as the conventional (Newtonian) calculus. In recent years, significant interest in fractional calculus has been stimulated by applications in physics, engineering, biology, and economics. The mathematical development of fractional calculus is reviewed in many places, e.g., [17,18]. The earliest proponent of fractional calculus in viscoelasticity appears to be Nutting [1], who observed that the stress relaxation phenomenon could be modeled by fractional powers of time. Later, Gemant [2] also found that stiffness and damping properties of viscoelastic materials appeared to be proportional to fractional powers of frequency, which led him to suggest the use of time differentials of fractional order. Bagley and Torvik [3,4] showed that the fractional calculus models of viscoelastic materials are in harmony with the molecular theories that describe the macroscopic behavior of viscoelastic media.

Generalized fractional calculus models for isotropic and anisotropic viscoelasticity are developed and reviewed in [17,18,19]. Compared to classical viscoelasticity, the fractional calculus models of viscoelastic constitutive relations employ derivatives (and integrals) that are of some non-integer order. Constitutive relationships employing derivatives of fractional order relating stress and strain lead to well-posed problems and have evolved as an empirical method of describing properties of many viscoelastic materials, including plastics, polymers, bituminous binder and cartilage [20,21,22,23,24]. Analogous and extensive studies have been conducted in diffusion-wave phenomena [25] and thermoelasticity, typically with an eye to improving the Fourier heat flow model [26]. Overall, a connection of the fractional response model to the underlying material microstructure—be it random or possibly percolating and fractal—is still missing.

3.2. Classical Rheological Models

Classical viscoelastic models require the constitutive relations to be Green’s functions associated with the integer-order derivative in the integral. Starting with

σ(t)∼^Dn[ϵ(t)]

we have two special cases: (i) when n=0 , the model is a spring element, whose stress–strain relation is σ(t)=E ϵ(t) ; (ii) when n=1 , the model is a linear viscous dash-pot element: σ(t)=η dϵ(t)/dt.

Classical linear viscoelastic material models are represented by a system of springs and dash-pots. Given the support by many experimental observations, one very popular model is the generalized Maxwell model, where a number of Maxwell models are connected in parallel with one isolated spring element introduced for the non-vanishing asymptotic relaxation modulus at infinite time (Figure 5a). Solution of classical linear viscoelastic models leads to exponential type responses. For example, the relaxation modulus for an N-series generalized Maxwell model is expressed as

E(t)=_E0+∑i=1n_Eiexp(−t_τi)

which is also called the Prony series. In a generalized Maxwell model with N+1 elements, there is a total of 2n+1 parameters: _E0,_E1,_τ1,…,_En,_τn . The advantage of classical rheological models is that exponential functions are easy to handle under mathematical operations such as Laplace and Fourier transforms. However, exponential functions predict a strong frequency dependence of damping properties, whereas measurements of viscoelastic materials reveal a very small change in their dissipative behavior with varying frequency [27]. Thus, classical viscoelastic models have to involve many higher-derivative terms to better agree with measured data, which generally leads to a (much) greater number of material parameters that need to be determined, e.g., [28].

Figure 5. (a) Generalized Maxwell model; (b) fractional calculus model (Smyth and Green, 2015).

3.3. Fractional Calculus Models

The basic building block of a fractional calculus model in viscoelasticity is a so-called spring-pot described by the equation

σ(t)=η^Dα[ϵ(t)]

where α∈(0,1) and the (conventionally adopted) fractional derivative of order α of a function x(t) with respect to time t is given by the Riemann–Liouville definition [29]:

^Dα[x(t)]=1Γ(1−α)ddt_∫0tx(τ)^(t−τ)αdτ, 0<α<1

where Γ(x)is the Gamma function. The Fourier transformation of this derivative leads to

ϝ[^Dαf(t)]=^(iω)αf˜

where f˜is the Fourier transform of function f. In the fractional calculus model, the spring-pot replaces the classical dash-pot of classical viscoelasticity.

Starting with the generalized Maxwell model, there is an additional free parameter of the derivative order _αi on each Maxwell model and the total number of parameters increases to 3n+1 . The relaxation modulus of the fractional calculus model in Figure 5b is given by [22]

E(t)=_E0+∑i=1n_{Ei εα}(−_{Ei ηi}^tα)

where

_εα(x)=∑k=0∞^xkΓ(αk+1).

Given that the time and frequency dependencies can be reduced with a decreasing derivative order α, the viscoelastic behavior of many polymers is generally modeled with fewer parameters using fractional calculus models than in the integer-order models.

4. Fractional Derivative in Elastic-Viscoelastic Composites

One important feature of a fractional calculus model is that it has the flexibility of a continuous passage from a constitutive relation of an ideal solid state to an ideal fluid state. By changing the fractional derivative order αfrom 0 to 1, the fractional calculus model smoothly transitions from an elastic spring to a viscoelastic dash-pot. This property can be utilized to calibrate the response of viscoelastic random composites.

Now, consider an elastic-viscoelastic composite whose viscoelastic phase is represented by a Zener model ( α=1 ); note here that the Zener model has the impact response. Decreasing α from 1 to 0 models the increasing content of the elastic phase. For any intermediate volume fraction of the viscoelastic versus elastic phase, the effective behavior of the composite should be well explained with a four-parameter fractional calculus model of Figure 5b. The corresponding constitutive equation is

σ(t)+η_E1×^Dα[σ(t)]=_E0ϵ(t)+_E0+_E1E1η×^Dα[ϵ(t)].

For the elastic-viscoelastic random composite with an arbitrary volume fraction, the behavior of the composite is fully captured by this equation with model parameters α∈(0,1] , _E0 , _E1 , and η . Upon taking the Fourier transform of Equation (16), it is seen that this model predicts a frequency-dependent complex modulus of the following form:

^E∗(ω)=_E0+_E0+_E1E1η^(iω)α1+η_E1^(iω)α.

Using the relation ^(iω)α=^ωαcos(πα/2)+isin(πα/2) , the parameters can be fitted by evaluating the loss factor (the tangent of the phase angle ϕ ) at multiple frequencies [4]:

η=(_E1−_E0b)sinπα2^ωα_E0+(_E1+_E0b)cosπα2^ωα+_E1b^ω2α.

For a quantitative study, we continue to use the random checkerboard, with responses of the composite at various volume fractions of the viscoelastic phase determined by the computational mechanics model of Section 2. In order to capture the broad range of frequency dependence and obtain a better fit of model parameters, a steady state analysis is performed from 0.01to 10 Hz for each composition.

Figure 6 displays the phase angle’s dependencies on frequency, plotted on semi-log and log-log graphs as best fits obtained through a least-squares method. In Figure 6a, we show the case of four volume fractions. In Figure 6b, we display the cases of critical volume fractions (i.e., percolation points) of elastic ( _vf=0.41 ) and viscoelastic ( _vf=0.59 ) phases; note that _vf pertains to the viscoelastic phase. Simulations have been conducted on 256×256 checkerboard domains, which was the largest size we could handle with our computational resources and allowed us to grasp the RVE-level responses of fractal patterns. At the same time, we considered high contrasts ( 1,000 ) to grasp the physics of critical phenomena. Overall, the fitted fractional order α is still in the range [0.99,1.0] for low contrast and [0.98,1.0]for high contrast.

Figure 6. Best fits of fractional calculus models in terms of the dependence of phase angle on the frequency: (a) at four different volume fractions of the viscoelastic phase; (b) at two critical volume fractions (0.41 and 0.59) corresponding, respectively, to the fractal percolations of elastic and viscoelastic phases.

Collecting all the results, it was found that, independent of the volume fraction of the viscoelastic phase, the fractional order does not decrease and that the best fit is always practically attained for integer α(=1). The frequency dependence of the elastic-viscoelastic composite remains the same as that of the viscoelastic phase, which is described through classical rheological models. The response of the composite is still in the form of the exponential function

E(t)=_E0+_E1 ^e−t/b

and the fitted parameters for different volume fractions are listed in Table 1. It can be seen that adding the elastic phase to the viscoelastic phase only increases the instantaneous modulus _E0and the time constant b, but it does not change the type of the viscoelastic model.

[ Table omitted. See PDF. ]

Overall, this study shows that, for one particular viscoelastic composite having a fractal structure but made of classical (non-fractional) viscoelastic continuum phases, the fractional calculus model is not necessary. In other words, the randomness of a composite material—even in the fractal regime (!)—is not the cause of the fractional order viscoelasticity. Rather, it is inferred that such a fractional response must be due to the fractional nature of the constituent phase(s).

5. Conclusions

In this paper, the passage between elasticity and viscoelasticity is investigated in elastic-viscoelastic random composites with strong contrast. For the random checkerboard, a sharp change in response occurs at the volume fraction ∼0.41, with the response switching from a time-independent to a time-dependent type. This critical volume fraction directly corresponds to the probability threshold of site percolation on ^Z2.

On a separate track, fractional calculus models are applied to describe the change in behavior of the elastic-viscoelastic random checkerboard. Upon examining the whole range of volume fractions, especially at the critical volume fraction, it is found that randomness and fractality of the checkerboard is not the cause of the fractional order viscoelasticity. Therefore, there is at least one random composite material whose response, even at the percolation point having a fractal structure, is adequately described by an integer-order viscoelastic model, which suggests that fractional-order viscoelasticity is likely due to the nature of the constituent phase(s) at the molecular level.

While fractal systems of rheological elements have been studied, they have not been studied from the standpoint of the mechanics of composite materials. Our study shows that a boundary value problem and homogenization in space-time is necessary to obtain an effective constitutive relation. At the same time, models of the aforementioned references are more appropriate for 1D systems. Of course, given that the present study involved one fractal system in 2D, a wealth of other composite media arrangements in 2D and 3D still needs to be analyzed to arrive at general conclusions to the question posed in the title.

Acknowledgments

Expert comments of Harry H. Hilton on fractional viscoelasticity have been helpful in the course of this research. This material is based upon work partially supported by the NSF under grant CMMI-1462749.

Author Contributions

Martin Ostoja-Starzewski conceived the problem; Jun Zhang performed, in consultation with Martin Ostoja-Starzewski, all the simulations and calculations; Jun Zhang and Martin Ostoja-Starzewski wrote the paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Nutting, P.G. A new general law of deformation. J. Frankl. Inst. 1921, 191, 679–685.

2. Gemant, A. A method of analyzing experimental results obtained from elasto-viscous bodies. Physics 1936, 7, 311–317.

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

4. Bagley, R.L.; Torvik, P.J. On the fractional calculus model of viscoelastic behavior. J. Rheol. 1986, 30, 133–155.

5. Bagley, R.L. Power law and fractional calculus model of viscoelasticity. AIAA J. 1989, 27, 1412–1417.

6. Mandelbrot, B.B. The Fractal Geometry of Nature; W. H. Freeman and Company: New York, NY, USA, 1982.

7. Heymans, N.; Bauwens, J.-C. Fractal rheological models and fractional differential equations. Rheol. Acta 1994, 33, 210–219.

8. Schiessel, H.; Blumen, A. Mesoscopic picture of the sol-gel transition: Ladder models and fractal networks. Macromolecules 1995, 28, 4013–4019.

9. Ionescu, C.M.; Kosiński, W.; Di Paola, M. Viscoelasticity and fractal structure in a model of human lungs. Arch. Mech. 2010, 62, 21–48.

10. Butera, S.; Di Paola, M. A physically based connection between fractional calculus and fractal geometry. Ann. Phys. 2014, 350, 145–158.

11. Ostoja-Starzewski, M. Microstructural Randomness and Scaling in Mechanics of Materials; Chapman & Hall/CRC Press: Boca Raton, FL, USA, 2008.

12. Ostoja-Starzewski, M.; Kale, S.; Karimi, P.; Malyarenko, A.; Raghavan, B.; Ranganathan, S.I.; Zhang, J. Scaling to RVE in random media. Adv. Appl. Mech. 2016, 49, 111–211.

13. Takayasu, H. Fractals in the Physical Sciences; Manchester University Press: Manchester, UK, 1990.

14. Zhang, J.; Ostoja-Starzewski, M. Mesoscale bounds in viscoelasticity of random composites. Mech. Res. Commun. 2015, 68, 98–104.

15. Zhang, J.; Ostoja-Starzewski, M. Frequency-dependent scaling from mesoscale to macroscale in viscoelastic random composites. Proc. R. Soc. 2016, 472, 20150801.

16. Ziman, J.M. Models of Disorder: The Theoretical Physics of Homogeneously Disordered Systems; Cambridge University Press: Cambridge, UK, 1979.

17. McBride, A.C. Fractional Calculus and Integral Transforms of Generalized functions; Pitman: London, UK, 1979.

18. Kiryakova, V.S. Generalized Fractional Calculus and Applications; CRC Press: Boca Raton, FL, USA, 1993.

19. Hilton, H.H. Generalized fractional derivative anisotropic viscoelastic characterization. Materials 2012, 5, 169–191.

20. Hilton, H.A.; Hsu, J.; Kirby, J.S. Linear viscoelastic analysis with random material properties. Probab. Eng. Mech. 1991, 6, 57–69.

21. Merrett, C.G.; Hilton, H.H. Fractional order derivative aero-servo-viscoelasticity. Int. J. Dyn. Control 2016, 5, 239–251.

22. Caputo, M.; Mainardi, F. Linear models of dissipation in anelastic solids. La Rivista del Nuovo Cimento 1971, 1, 161–198.

23. Slonimsky, G.L. Laws of mechanical relaxation processes in polymers. J. Polym. Sci. Polym. Symp. 1967, 16, 1667–1672.

24. Mainardi, F. Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models; World Scientific: Singapore, 2010.

25. Povstenko, Y. Linear Fractional Diffusion-Wave Equation for Scientists and Engineers; Birkhäuser: Basel, Switzerland, 2015.

26. Povstenko, Y. Fractional Thermoelasticity; Springer: Berlin/Heidelberg, Germany, 2015.

27. Jones, D.I.G. Handbook of Viscoelastic Vibration Damping; John Wiley & Sons: Hoboken, NJ, USA, 2001.

28. Smyth, P.A.I.; Green, I. Fractional calculus model of articular cartilage based on experimental stress-relaxation. Mech. Time Depend. Mater. 2015, 19, 209–228.

29. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of their Solution and Some of their Applications; Academic Press: Cambridge, MA, USA, 1998.