Does a Fractal Microstructure Require a Fractional Viscoelastic Model ?

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


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

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 v f signifying the content of the viscoelastic phase, v f = 0% means "all elastic", while v f = 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 by 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: Hence, for a square δ × δ lattice, the number of different realizations (or the size of the entire sample space Ω) is |Ω| = 2 δ×δ .Each elementary event ω gives rise to one realization B δ (ω) of a random medium B δ = {B δ (ω) ; ω ∈ Ω}.B δ (ω) occurs with probability 1/2 δ×δ .Effectively, according to the probability p, the volume fraction v f of the viscoelastic phase is determined.
Following [11,14,15], a computational mechanics method is employed to determine the response of any specific 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 for weaker fluctuations, and δ = 128 and 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) is characterized by the relaxation shear modulus µ (t) and the relaxation bulk modulus k(t) where J and K are the spherical and deviatoric parts of the unit fourth-order tensor.The detailed time-dependence of µ(t) and k(t) can be further modeled numerically via the Prony series where µ 0 and k 0 are instantaneous moduli, with g n , k m and τ 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 (g n , k m , τ 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, 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 First, consider two boundary value problems for the body B δ (ω): • kinematic uniform boundary conditions (KUBC): • static uniform boundary conditions (SUBC): where σ 0 ij (t) and 0 ij (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: • under SUBC, letting The relaxation stress and creep strain, after stabilization (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.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% (v f > 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% (v f < 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.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 p c 0.59.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.

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.

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 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 which is also called the Prony series.In a generalized Maxwell model with N + 1 elements, there is a total of 2n + 1 parameters: E 0 , E 1 , τ 1 , ..., E n , τ 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].

Fractional Calculus Models
The basic building block of a fractional calculus model in viscoelasticity is a so-called spring-pot described by the equation 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]: where Γ(x) is the Gamma function.The Fourier transformation of this derivative leads to 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] where 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.

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 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], E 0 , E 1 , 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: 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]: 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.01 to 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 (v f = 0.41) and viscoelastic (v f = 0.59) phases; note that v f 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.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 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 E 0 and the time constant b, but it does not change the type of the viscoelastic model.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).

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

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

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

Figure 4 .
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.

Figure 6 .
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.

Table 1 .
Parameters of fractional calculus model.