Effective properties of fractional viscoelastic composites via two‐scale asymptotic homogenization

Driven by the growing interest in fractional constitutive modeling, we adopt a two‐scale asymptotic homogenization approach to study the effective properties of fractional viscoelastic composites. We focus on a purely mechanical setting and derive the cell and homogenized problems corresponding to the balance of linear momentum equation in the absence of body forces and inertial terms. In doing this, we reformulate the original framework in the Laplace–Carson domain and discuss how to obtain the effective coefficients in the time domain. We particularize the general setting of our work by considering memory functions that describe special types of fractional linear viscoelastic behaviors, and after presenting the limit cases of our selections, we framed the homogenization results to account for benchmark problems with different combinations of constitutive models. Specifically, these latter involve elastic, fractional Kelvin–Voigt, fractional Zener and fractional Maxwell constituents. Our results permit us to reinforce the interpretation of the theoretical findings and to elucidate the role of the fractional constitutive models on the effective properties of the composites under investigation.


INTRODUCTION
Recent advances in the investigation of materials possessing a complex micro-structure have shown that their time-dependent constitutive behavior can be described in terms of power laws (see, e.g., [1][2][3]).Within the context of linear viscoelasticity, spring-dashpot networks [3] have been used to mimic power-law models.However, the definitions provided by the theory of fractional calculus [4,5] have permitted to capture power-law behaviors by introducing a number of fractional parameters.
Fractional viscoelasticity has arisen as a promising tool to describe and analyze experimental data related to the rheological behavior of various heterogeneous systems found in biology [6,7], and materials science [8,9], among others.Together with successful experimental data fitting, there have been significant efforts to answer theoretical questions on the foundations of linear fractional viscoelastic models.For instance, Bologna et al. [10] provide thermodynamic restrictions for the orders of the power laws used to characterize the fractional-order hereditary behavior of hydroxyapatite-based composites.On the other hand, the thermodynamical restrictions for a linear constitutive equation containing fractional derivatives of stress and strain of different orders are discussed in Atanackovi ć et al. [11] by using the entropy inequality for isothermal processes.
In spite of the examples reported above, and the several studies concentrating on the effective properties of linear viscoelastic composites (see, for instance, [12][13][14][15][16][17]) using different homogenization approaches, to the best of our knowledge, fractional viscoelasticity is considerably underused in the investigation of the overall behavior of composite materials.For example, Gallican and Brenner [18] propose estimates for the effective response of viscoelastic composites with fractional Zener constituents via micro-mechanical modeling.On the other hand, in Ostoja-Starzewski and Zhang [19], the Hill-Mandel condition is used to determine the overall response of a random viscoelastic micro-structure.In the context of asymptotic homogenization, we mention, for instance, the works [20] and [21] where closed formulas for the effective coefficients of fiber-reinforced viscoelastic and visco-piezoelastic composites are obtained based on the Rabotnov's fractional exponential kernel.Finally, we mention the paper by Wang and Pindera [15], in which, even though the notion of fractional viscoelasticity is not considered, the authors employ the locally exact homogenization theory to investigate the effective properties of a unidirectional composite with the viscoelastic matrix exhibiting a power-law constitutive behavior.
We further remark that it is still lacking a well-established theoretical study connecting the role of the micro-structure in the origins of fractional behavior.In this respect, we mention the paper by Brenner [22] where the author studies the emergence of a fractional viscoelastic behavior resulting from micro-structural information considering Hashin-Shtrikman and self-consistent estimates.We also mention that, in recent investigations, there have been attempts to propose variable-order fractional parameters depending on specific fields or micro-structural information.For instance, to study the effect of temperature on the mechanical behavior of a linear fractional viscoelastic material, Colinas-Armijo et al. [23] propose fractional orders explicitly depending on the temperature.On the other hand, in Ramirez and Coimbra [24], the time-dependent order of the fractional derivative is assumed to be a measure of the rate of change of disorder within a composite material.In particular, the results derived from the statistical mechanical model proposed in Ramirez and Coimbra [24] where in agreement with experimental studies conducted in an epoxy resin and a carbon/epoxy composite.
Here, motivated by previous works (see, e.g., [17,20,21,[25][26][27]), we use the asymptotic homogenization method [28][29][30] to compute the effective properties of a composite with viscoelastic constitutive response.In particular, we conduct our calculations in the general setting of linear viscoelastic materials and specialize them for the case in which the local constitutive response of each constituent is described by a fractional viscoelastic model.Although the approach has common points with that followed in Rodríguez-Ramos et al. [20] and Otero et al. [21], we set the problem differently such that it is possible to explicitly split the elastic and the memory contributions of the composite.This is one of the novelties of this work.Furthermore, we analyze the consequences of limit cases by using a unified memory law which includes the fractional Kelvin-Voigt, Maxwell, and Zener models and that can be generalized to include Prony series and Rabotnov's kernels.In addition, we find analytical expressions for the effective coefficients for special cases of composite media.Within this context, we find that for a layered composite with elastic and fractional Kelvin-Voigt constituents the effective coefficients can be written in terms of a two-parameter Mitta-Leffler function.This constitutes another novelty of this work.
The present work is organized as follows.In Section 2, we introduce the main equations governing our model.In Section 3, we reformulate the original problem in a two-scale fashion and introduce the topology of the composite micro-structure.Afterwards, in Section 4, we derive the cell and homogenized problems and the general form of the effective coefficients for a two-phase viscoelastic composite via the asymptotic homogenization technique.In doing this, we transform the two-scale problem to the Laplace domain using the elastic-viscoelastic correspondence principle and take advantage of the Laplace-Carson transform in the inversion of the effective coefficients to the time domain.Furthermore, we use these theoretical results to derive the analytical expressions for the effective elasticity and relaxation tensors in the case of a composite medium with a three-dimensional layered structure.In Section 5, we present specific constitutive models for linear viscoelastic and fractional viscoelastic materials and discuss some limit cases of the latter.Moreover, we specialize the general theory by presenting different benchmark problems and discussing the results produced by analytical and numerical simulations.Finally, in Section 6, we highlight the main results and outline some future developments.

PROBLEM FORMULATION
Let us denote by B an open and bounded set of the three-dimensional Euclidean space, S, representing a continuum body with periodic micro-structure and constituted by two different solid sub-phases denoted by B 1 and B 2 .The two open subsets B 1 and The bar notation indicates that we consider the topological closure of the subset.Furthermore, we denote by I the interface between B 1 and B 2 and by N I the unit vector normal to I, where the positive orientation is taken from B 1 to B 2 .Additionally, for each point X ∈ B, we associate a triple of real numbers represented by (X 1 , X 2 , X 3 ) so that they represent, with abuse of notation, the coordinates of X in a Cartesian reference system.
In the absence of body forces, the point-wise balance of linear momentum for the composite B reads div where S  (X, t), with  = 1, 2, denotes the Cauchy stress tensor related to the constituent B  and (X, t) ∈ B  × [0, t f [, with t f representing a positive time.In addition, we prescribe continuity of displacements and tractions on the constituents' interface, namely, for all t ∈ [0, t f [, where we have used the notation with X − ∈ I − (X + ∈ I + ) and I − (I + ) being the set of points of I interpreted as the points of I which are in contact with the sub-phase B 1 (B 2 ).The problem defined by Equations ( 1), (2a), and (2b) needs to complemented with appropriate boundary and initial conditions.However, since here we are focused on computing the effective material properties, they do not need to be specified.We also notice that the notation used so far is not general to all possible topological configurations and, thus, may need to be adapted, for example, to include the case of disconnected subsets [31].By adhering to the infinitesimal theory of viscoelasticity [32,33], the stress tensors, S  (X, t), for non-aging viscoelastic media in the absence of residual stresses can be written as where E  (X, t) ∶= sym(grad U  (X, t)) represents the infinitesimal strain tensor, U  is the displacement field, grad U  = [U  ] a ∕X b i a ⊗i b denotes the displacement gradient, and {i a } 3 a=1 is the orthonormal vector basis in Cartesian coordinates.The fourth-order tensors C  (X) and L  (X, t) are referred to as the elasticity tensor and the stress-relaxation tensor [33,34], respectively, which we consider to enjoy of both left and right minor symmetries, and major symmetry, that is, with a, b, c, d = 1, 2, 3.In particular, under the hypothesis of causal histories, that is, by requiring the infinitesimal strain to be zero for t < 0 and non-zero in the interval [0, t[, we can write Ē (X, t) = H(t) Ē (X, t), where H(t) denotes the Heaviside step function and Ē (X, t) is the infinitesimal strain which can have non-zero values for t < 0. Therefore, Equation (4) can be rewritten as For simplicity of notation, in the remainder of this work, we will write E  instead of Ē .In particular, the constitutive expression in Equation ( 6) can be equivalently rewritten in the form where G  (X, t) ∶= C  (X) + L  (X, t).We remark that we follow the approach outlined by Eringen [33] for the presentation of the general form of the viscoelastic constitutive law.Therefore, we adopt the notation and terminology used in Eringen [33].In doing this, instead of referring to G  as the relaxation tensor (as it is usually done), we prefer to call L  the relaxation tensor, since it is the one that encodes the memory properties of the constituents.

MULTI-SCALE FORMULATION
Here, we reformulate the equations specified in the previous section in a multi-scale fashion and describe the topology of the micro-structure for the problem under consideration.

Separation of scales
We start by denoting with L c and  the characteristic length scales of the composite medium and of its internal structure so that where  is referred to as the smallness parameter.Moreover, following Ramírez-Torres et al. [26] and Di Stefano [31], we formally rewrite a given physical quantity Φ ∶ B × [0, t f [→ R in a two-scale fashion, so that the dependence on the characteristic length scales is explicitly taken into account.Specifically, the multi-scale version of a quantity Φ(X, t) is written as [26,31] Φ(X, t) = (x, , t), (9) where the dimensionless variables x ∶= X∕L c and  ∶= X∕ are referred to as the macroscopic variable and the microscopic variable, respectively.We notice that, in this framework, the partial derivatives of Φ with respect to the spatial coordinates X i , i = 1, 2, 3, of X can be expressed as For the problem under study, we introduce the notations for the displacement and stress fields, and for the elasticity and the relaxation tensors.Furthermore, we set G  (X, t) =   (x, , t).

Topology of the micro-structure
As usually done in the context of the asymptotic homogenization technique [28,29], we assume that all the fields of interest are periodic with respect to the micro-scale variable  and introduce an elementary cell, which we denote Y in a non-dimensional formalism.Specifically, we consider that For the sake of simplicity, we adopt the assumption of macroscopic uniformity (see, for instance, Burridge and Keller [35] and Penta et al. [36]), which permits to choose the elementary cell independently of the position in the domain B. This consideration implies that the elementary cell, Y, and, thus, the sub-phases Y 1 and Y 2 as well as the interface I are representative of the micro-structure of B. Consequently, the normal vector n Y results to be independent of the macroscopic variable x.

Two-scale governing equations
In view of the above discussions, we can reformulate the original problem given by Equations ( 1)-(2b) as for the generic field (x, , t), with being the set of points of I Y interpreted as the points of I Y which are in contact with the sub-phase Y 1 (Y 2 ).Furthermore, for  = 1, 2, we have that Before going further, we transform the problem specified in (13a)-(15c) from the time domain to the Laplace transform domain by considering the elastic-viscoelastic correspondence principle (see, e.g., [15,17,37] in the context of asymptotic homogenization).Specifically, by introducing the Laplace transform of the generic field (x, , t), namely, the system of Equations (13a)-(13c) can be rewritten as where, by considering the expressions in formulae (15a)-(15c), the Laplace transform of the two-scale, infinitesimal stress tensor can be written as where ĝ denotes the Laplace-Carson transform of , which is defined as [38]  for the generic field .In particular, since the elasticity tensor is time independent, ĉ (x, , s) =   (x, ), and thus, we can write ĝ (x, , s) = ĉ (x, , s)

THE ASYMPTOTIC HOMOGENIZATION TECHNIQUE
Following the standard procedure in the asymptotic homogenization technique [28,29], we prescribe a formal two-scale expansion for the displacement ũ in power series of the smallness parameter  as follows: where u (i)   represent smooth, Y-periodic fields.The substitution of the expansion ( 21) into (17a) leads, after multiplication by  2 and equating in powers of the smallness parameter , to a set of differential equations holding in the periodic cell, Y, and parametrized by the macro-scale variable x.These equations, in particular, are supplemented with interface conditions resulting from the substitution of the formal expansion ( 21) into (17b)-(17c).In the following, we report the set of problems associated with the asymptotic homogenization procedure.Although the description of the homogenization technique presented in this section is classical, we report its main steps since the presentation of the problem at hand differs from those typically encountered in the existing literature.

The first cell problem
The first cell problem reads where, for i = 0, 1, 2, … , we have introduced the notation The system (22a)-(22c) represents a linear elastic-type boundary problem with zero source term, for which the only periodic solutions are of the form [28,29]

The second cell problem
The second cell problem can be written as where e (0)  x (x, s) ∶= L −1 c sym(grad x u (0) (x, s)).The system (25a)-(25c), which can be also regarded as a linear elastic-type boundary problem in the unknown field u (1)   , admits a unique solution up to an additive -constant field [30].Consequently, we can express u (1)   through the ansatz where   denotes a -periodic third-order tensor field and w  is a vector field.We notice that using the expression introduced in Equation ( 26), we can write e (1)   (x, , s) ∶= sym(grad  u (1)   (x, , s)) =   (x, , s) ∶ e (0)  x (x, s), (27) where the components of the fourth-order tensor field   are From the above expression, we can verify that   enjoys of the left-minor symmetry, that is, [  ] abcd = [  ] bacd .Now, by substituting the ansatz (26) in Equations (25a)-(25c) and taking into account the expression in (27), the second cell problem, usually referred to as the local problem, reads We remark that, to solve the cell problem (29a)-(29c), we take advantage of the symmetry of the tensor fields in these equations to reduce the number of unknowns.Looking at Equation ( 28), we need to find the 81 components of the fourth-order tensor   .However, by definition, the tensor   enjoys the left-minor symmetry.On the other hand, the symmetry of e (0)  x implies that  is symmetric in its last two indices and, consequently,   also enjoys right minor symmetry.Therefore, the number of unknown components of the tensor   reduces to 36.We refer the reader to Di Stefano et al. [31] and Penta and Gerisch [39] for further details.

The homogenized problem
By considering the third cell problem, we can write x (x, s) + ĝ (x, , s) ∶ e (1)   (x, , s) ] (1)  x (x, , s) + ĝ (x, , s) ∶ e (2)   (x, , s) where, for i = 0, 1, 2, … , Before going further in our analysis, let us introduce, for a generic field   (x, , s), with  = 1, 2, the operator where |Y | represents the measure of Y  .In particular, we notice that the average of  over the cell Y is where  ∶=  1  Y 1 +  2  Y 2 and  Y  denotes the indicator function of Y  .Consequently, by applying the operator introduced in (32) to (30a) and considering the assumption of macroscopic uniformity, which permits to commute the integral operator over the periodic cell and the divergence with respect to x, we can deduce that x (x, s) + ⟨ ĝ (x, , s) ∶ e (1)   (x, , s) Thus, summing up over  = 1, 2 the second addend in (34) and considering the periodicity properties of u (1)   and u (2)   , we have that (1)  x + ĝ ∶ e (2)    ]⟩ ĝ ∶ e (1)  x + ĝ ∶ e (2)

𝜂𝑦
] n  d, (35) where we have also invoked the interface condition (30c) and used the notation n 1 ∶= n Y and n 2 ∶= −n Y .We remark that the result (35) holds true for a certain class of periodic cells with specific topological characteristics; see, for example, Di Stefano et al. [31] for further discussion on this issue.Therefore, by summing Equation ( 34) over  = 1, 2, and taking into account (35) and (27), we obtain that Equation ( 34) can be equivalently rewritten as where we have introduced the auxiliary notation We notice that Equation (36) constitutes the homogenized equation for the unknown u (0) (x, s) and holds true in the spatial domain B h , namely, in the homogenized version of B.
Remark 1.We remark that the fourth-order tensor  eff can be written as where  eff and  eff are defined through the expressions and represent, respectively, the effective elasticity and the effective relaxation tensors.It is worth noticing that, in general, the fields u (i) may not represent the result of the application of the Laplace transform to a certain vector field and, thus, the transformation of the effective coefficient to the time domain should be carefully done by selecting an appropriate operator.This, however, is not an easy task as we would want to select a transformation operator which permits us to obtain the equivalent homogenized properties corresponding to the original problem.Here, to avoid this uncertainty, we assume that u (i) represents the result of a Laplace transformation, such that u (i) (x, , s) ≡ ũ(i) (x, , s).Therefore, working with the second addends in (39a) and (39b), we can deduce that where the time convolution is understood between the components of the fourth-order tensor.Then, we can conclude that both  eff and  eff represent the result of Laplace-Carson transforms and we can write Consequently, to obtain the value of  eff in the time domain, we need to consider the inverse Laplace-Carson transform.In the remainder of this work, we will consider the assumptions leading to Equation (41).

Effective coefficients for viscoelastic layered composites
Our main scope is to show how fractional constitutive laws influence the effective properties of a composite.Thus, although we could consider geometric settings more complex than the one factually employed in the sequel (see, e.g., Cruz-González et al. [17]), here, we specialize our theory to the case of a layered composite medium in order to give prominence to the role of the use of a fractional constitutive law with respect to geometric complexity.We notice that the effective coefficients that we obtain for our problem by using the Laplace-Carson transform for the case of strain and stress fields uniform in each layer were previously determined in Backus [40].
The right minor symmetry of the tensor ĝ allows to rewrite the cell problem (29a)-(29c) in the form ) where [41].In particular, the symmetry of the third-order tensor   in its last two indices reduces the number of unknowns of the cell problem from 27 to 18.Then, for a three-dimensional layered composite medium with layers orthogonal to the i 3 direction and with each layer assumed to be homogeneous in  1 and  2 , and, thus, in each plane orthogonal to i 3 , the material properties of the composite change solely along the i 3 direction, and thus, the first equation in the cell problem (42a)-(42c) reduces to where, for simplicity of notation, we have replaced  3 with .Within this framework, the original problem, which is three-dimensional, can be reconceived as one-dimensional and, accordingly, the unit cell can be associated with the interval Y =]0, 1[ and the homogenized domain with B h =]0, L∕L c [, with L > 0. Due to the simplified structure of Equation ( 43), a direct integration with respect to  yields that the term inside the braces is independent of , that is, where   is the "integration constant."This, in turn, leads to where the auxiliary notation h −1  is used to represent the inverse of the second-order tensor with components [h  (x, , s)] ab = [ ĝ (x, , s) ] a3b3 .We notice that, to solve the cell problem, one should determine   .However, since we are interested only in the determination of the effective coefficients, this is not necessary.Indeed, looking at Equations (39a) and (39b), one can just determine   , which reduces to the partial derivative of the third-order tensor   , given by Equation (45).Consequently, from Equation (45), we only need to find the "integration constants" supplied by the components [  ] p3cd .In particular, the substitution of (45) into the interface condition (42c) leads to Thus, by applying the integral operator defined in (32) to (45), summing up over  = 1, 2, and taking into account the local periodicity of   , as well as its continuity at the interface, we deduce that with We notice that the tensorial coefficients of   represent partition coefficients, so that Equation ( 47) defines [] r3cd as the weighted average of [ ĝ ] n3cd .Now, by substituting ( 47) into (45), the effective elastic and relaxation tensors for the three-dimensional layered composite medium are given by the expressions and In particular, by considering Equations ( 49) and ( 50), we can obtain an expression for ĝeff by means of Equation ( 41).

BENCHMARK PROBLEMS
Let us assume to be in the presence of a layered composite with isotropic viscoelastic laminates such that, for  = 1, 2,   (x, ) = 3 e  (x, ) + 2 e  (x, ), (51a) where  e  and  e  represent, respectively, the elastic bulk modulus and the second Lamé's parameter.Furthermore,  v  and  v  denote memory functions, which need to be expressed constitutively.The fourth-order tensors  ∶= , extract, respectively, the spherical and the deviatoric part of a symmetric second-order tensor.Specifically, where A is a generic second-order symmetric tensor.Moreover, the component expressions corresponding to the tensor products I⊗I and I⊗I are (see, for instance, Federico [42]) where  ab denotes the Kronecker delta symbol.

Linear viscoelastic and fractional viscoelatic constitutive laws
Here, following Mainardi [4], we specify the constitutive expressions for some types of viscoelastic and fractional viscoelastic materials.
where (t) is the Dirac delta distribution (with units of inverse of time) and    and    are referred to as relaxation times.
We remark that both a   and a   have physical units where Σ denotes a characteristic stress unit and T is a characteristic time unit.That is, they represent dynamic viscosities.On the other hand, b   and b   have physical units = Σ, that is, they denote elastic parameters.The constitutive representations in (54a) and (54b) admit the possibility that for some materials the viscoelastic spherical and deviatoric contributions can be represented by different viscoelastic models [43].In Table 1, we summarize the types of viscoelastic models that can be specified from (54a) and (54b).We remark that the expressions specified above can be generalized to include other linear viscoleastic constitutive laws (see, for instance, Mainardi [4]).The symbols    and    have been used to represent the dynamic viscosities, and  M  ,  Z  ,  M  , and  Z  are elastic parameters associated with the Maxwell and Zener models.
Remark 2. We notice that the elastic parameters  M  and  M  and  Z  and  Z  are, in general, different from  e  and  e  .Specifically, in the case of the Zener model, we need to comply with the condition at t = 0, [4] 0 <  e  <  Z  and 0 <  e  <  Z  .
(55) Furthermore, we remark that using the coefficients specified in Table 1 associated with the Maxwell model, the tensor   can be written in the more standard form That is, to obtain the expression in (56) from (54a) and (54b), we included the elastic term with elastic parameters  e

Linear fractional viscoelastic models
Let us assume that the elasticity tensor,   , and the relaxation tensor,   , have the same form as in (51a) and (51b).However, in this case, we consider that the constitutive expression for the memory functions  v  and  v  are of fractional type.Taking inspiration from [4], we write where    ,    ∈]0, 1[, Γ(•) is the Gamma function and     and     denote one-parameter Mittag-Leffler functions of t∕   and t∕   .We notice that in Equations ( 58a) and (58b), the dependence of a   , a   , b   , and b   on    and    is understood but omitted for reducing the notational burden and the dot "•" indicates that these functions must be evaluated in (x, ).Considering the discussions in the previous section, some types of fractional viscoelastic models can be obtained by appropriately fixing the values of the coefficients in (58a) and (58b).In Table 2, we summarize such models.
Before going further, we remark that the fractional parameters, in general, can also depend on x and .However, in the writing of the memory functions  v  and  v  , we have preferred not to specify this dependence for the sake of simplicity.We also mention that we are tacitly assuming that the fractional parameters are time independent so that we can compute the Laplace transform of the memory functions by using classical results.
Remark 3 (Physical units).By referring to Table 2, we notice that both a   and a   have physical units Consequently, in the limit, lim Thus, in the first case, the coefficients a   and a   become dynamic viscosities, while in the second case, they reduce to elastic parameters.We remark that it is possible to design experiments that measure anomalous dynamic viscosities (see, e.g., Bonfanti et al. [3]) represented by the coefficients a   and a   .Here, however, to overcome the use of an anomalous dynamic viscosity, whose physical meaning may be unclear, we consider the expressions introduced in Table 2.That is, we write, for the fractional Kelvin-Voigt model where the characteristic times    and    are defined as in Table 2, namely,    =    ∕ e  and    =    ∕ e  .Finally, we remark that both b   and b   represent elastic parameters.
Remark 4 (Limit cases).In the limit, in which the fractional parameters    and    tend to one from below, the fractional viscoelastic models accounted for in this work and specified in Table 2 reduce to the standard linear viscoelastic models examined in Table 1 (see [3,4,11,44] for further details).On the other hand, if we take the limit {   ,    } → 0 + , both the fractional Maxwell and Zener models reduce to constitutive expressions of linear elastic type.Specifically, we have that, in the limit {   ,    } → 0 + , the Laplace transform of the infinitesimal stress tensors corresponding to fractional constitutive laws of Maxwell and Zener types take on the form (see Equations ( 78) and (80) below) respectively, where On the other hand considering the coefficients in Table 2 corresponding to a fractional Kelvin-Voigt material, the limit {   ,    } → 0 + of the Laplace transform of the infinitesimal stress tensor yields with The expression for the limit of σ specified in (64) suggests the existence of a second-order partial time derivative in the constitutive expression of the infinitesimal stress tensor and, thus, a damping effect (see, e.g., Atanackovic et al. [45] and Vaz and de Oliveira [46]).

Numerical results
For the sake of convenience, let us consider that the elastic parameters  e  and  e  , as well as the memory functions  v  and  v  are uniform in the macro-scale variable and piecewise constant in the micro-scale variable, that is, and and and (67)

Benchmark problem I: Fractional viscoelastic-elastic composite
In this first benchmark problem, let us assume that one constituent is constitutively determined by a fractional viscoelastic law, for example, B 1 , while the other constituent, B 2 , is linear elastic.Specifically, we assume that B 1 is characterized by a fractional Kelvin-Voigt viscoelastic law.Examples of the use of the fractional Kelvin-Voigt model in practical cases are summarized in Bonfanti et al. [3] and include, but are not limited to, the modeling of human prostate [47] and of breast cells [48].
Using the coefficients defined in Table 2 for the fractional Kelvin-Voigt model, the Laplace transforms of the memory functions are which, by considering Equation (51b), leads to the following expression for the Laplace transform of the relaxation tensor where we notice that because of (66) and (67), the Laplace transform of the relaxation tensor,  1 , is uniform in the macro-scale variable x.
In particular, we notice that, by using formula (37) and the representations (51a) and (51b), the effective coefficient [ ĝeff ] 3333 can be written in the form where Thus, by substituting (71a) and (71b) in (70), the effective coefficient [ ĝeff ] 3333 can be written in the simplified form In the specific case of the present benchmark problem, since Y 1 represents a fractional Kelvin-Voigt material and Y 2 is of elastic type, we have that Consequently, by considering that  ∶=   1 =   1 and applying the inverse Laplace-Carson transform, we can find an analytical expression for the effective coefficient [ eff ] 3333 , which reads where  2−,1 (•) ≡  2− (•) denotes the two-parameter Mittag-Leffler function in the case in which the second parameter is 1, and [] 3333 constitutes the effective coefficient of an equivalent layered composite made of two different elastic constituents, that is, In Figure 1, we plot the effective coefficient [ eff ] 3333 for different values of the fractional parameter, namely,  ∶=   1 =   1 ∈ {0.01, 0.2, 0.4, 0.6, 0.8, 0.99}, in the time interval [0s, 10 s], where, to highlight the role of the fractional constitutive properties, we assumed that both constituents have the same volumetric fraction, namely, |Y  | = 0.5,  = 1, 2. For our analysis, we used the parameter values reported in Table 3, which have been taken from Alotta et al. [43].In particular, for  near 1, we recover the standard Kelvin-Voigt exponential decay, so that the effective coefficient converges exponentially towards the value [] 3333 , which is the equivalent coefficient that the homogenized medium would have if both of its constituents were elastic materials with the same elastic properties as the materials employed in our model.As  decreases towards zero, the decay is characterized by the Mittag-Leffler function of Equation (74), which implies that the decay is accompanied by an oscillatory behavior.In fact, for  → 0 + , the oscillatory trend tends to increase, so that the convergence towards [] 3333 slows down considerably.Indeed, for  = 0, no decay is observed, and [ eff ] 3333 oscillates indefinitely around [] 3333 , so that the viscous effect is lost.We recall that, as discussed in Remark 4, the oscillatory trend is related to the fact that in the limit  → 0 + , the material property lim →0 + ĝ1 features a damping effect due to the presence of the factor s 2 in (64), which is equivalent to a second-order time derivative.Although it is not shown here, an oscillatory decay for the effective coefficient  eff is predicted by our results also if a fractional Kelvin-Voigt constituent is combined with a fractional Maxwell or fractional Zener one.Furthermore, from Equation (74), we notice that, at time t = 0s, the value  of the effective coefficient of our model is independent of  and coincides with the elastic coefficient of the material B 2 , divided by |Y 2 |.This behavior is due to the fact that, in our model, we have assumed that the composite is made of an elastic medium, that is, B 2 , and a viscoelastic one, that is, B 1 .Indeed, if we had considered the case of two viscoelastic materials, Equation (74) would have been different, and consequently, the effective coefficient at t = 0 s would have been a combination of the elastic coefficients of both materials.We finally notice that the oscillatory decay seems to be a property of the geometrical setting since, in a previous work [17], we did not observe this trend when considering a fiber-reinforced composite featuring power-law constitutive properties.We continue our analysis by plotting the effective coefficient [ eff ] 3333 for different values of the volumetric fraction of the fractional Kelvin-Voigt constituent.In the case reported in Figure 2, we fix  = 0.1 and notice that, for decreasing values of the fractional Kelvin-Voigt constituent's volumetric fraction, the oscillatory behavior of the effective coefficient decreases and tends, as expected, to the effective coefficient of the equivalent elastic composite, that is, a composite made of two elastic constituents with the same material parameters reported in Table 3.
Before going further, we notice that, because of the simplicity of the setting considered so far, i.e. a composite made of a fractional Kelvin-Voigt constituent with a single fractional parameter  and an elastic one, it was possible to find an analytical solution for the effective coefficient [ eff ] 3333 .Nevertheless, for more complex cases, finding such a solution is not a simple task.For this reason, in what follows, we adopt the methodology described in previous works (see, for instance, Cruz-González et al. [17,25]), in which the MATLAB function INVLAP [49,50] was adapted to numerically find the inverse Laplace-Carson ĝeff .We emphasize that we choose this function because of its power laws of the type s  where  is a real number.order to test the accuracy of this numerical tool, we compare, in Figure 3, the numerically obtained coefficient [ eff ] 3333 with the analytical expression given in As it can be observed, both the numerical and the analytical solutions agree for each value of the fractional parameter.

Benchmark problem II: Fractional viscoelastic composite
In this section, we consider that the constituents of the composite medium are characterized by different fractional viscoelastic models.In particular, we examine the case in which B 1 is assumed to be of fractional Maxwell type, whereas B 2 is characterized by a fractional Zener law.Examples of the use of these models in biologically relevant scenarios include the modeling of atrial tissue [6] and brain tissue [51].Here, we numerically compute the effective properties of the composite described above by considering the expressions obtained in (49) and (50) for different situations.In particular, the memory functions for the constituent B 1 (see Equations (58a) and (58b), and Table 2) are which, by means of the Laplace transform, can be written as Therefore, the Laplace transform of the relaxation tensor is curves in Figure 4 show a monotonically decreasing trend of the reported effective properties, which is characteristic of composite materials composed of an elastic constituent and a Maxwell/Zener one (see, for instance, Cruz-González et al. [17]).Furthermore, in Figure 5, we show the influence of different volumetric fractions of the fractional Maxwell constituent on the effective coefficient [ eff ] 3333 .For this purpose, we set   1 =   1 = 0.99 and   2 =   2 = 0.01 (panel on the left) and   1 =   1 = 0.01 and   2 =   2 = 0.99 (panel on the right).This means that the results on the left panel of Figure 5 correspond to a composite constituted by a quasi-standard viscoelastic Maxwell material and a quasi-elastic one, while on the right panel of Figure 5, the same effective coefficient is evaluated for the case of a quasi-standard elastic-viscoelastic Zener composite.Particularly, we notice that the instant of time at which the effective properties of the composite become stationary depends not only on the fractional parameter but also on the constituents' volumetric fractions.

CONCLUSIONS
In this work, we employ the asymptotic homogenization to compute the effective properties of a hypothetical linear viscoelastic composite medium.Specifically, after laying out the main theoretical results, we specialize them to the case of a laminated composite medium, each constituent of which is isotropic and constitutively described through a linear viscoelastic model.This particular geometrical setting permits to reduce the computational complexity of the problem, since it allows to obtain the analytical expressions of the effective elasticity and relaxation tensors (refer to Equations ( 49) and ( 50)).
Motivated by the interest that fractional constitutive modeling has drawn in the investigation of heterogeneous, complex materials, we considered memory functions describing special types of fractional linear viscoelastic behaviors.Following the discussions of the limit cases in Remarks 3 and 4, we accounted for benchmark problems with different combinations of linear viscoelastic laws.In particular, in a first benchmark problem, we assumed a composite medium made of an elastic and a fractional Kelvin-Voigt viscoelastic medium, while in the second benchmark problem, we studied the case of a composite with fractional Maxwell and fractional Zener viscoelastic constituents.We point out that the number of benchmark cases investigated in this work can be extended to include other combinations.
We emphasize that, in our study, we gave prominence to the influence that the fractional parameters    and    , and the constituents' volumetric fractions |Y  | exert on the time-dependent behavior of the effective properties of the homogenized medium.In particular, with reference to the first benchmark problem, for    =    , we found an analytical expression for the effective coefficient [ eff ] 3333 in terms of a two-parameter Mittag-Leffler function that entails, as observed in Figure 1, a transition from an exponential decay trend to an oscillatory one, as the limit  → 0 + is approached.Furthermore, we analyzed the influence of the constituents' volumetric fraction on the effective properties.For this purpose, we fixed the fractional parameter to be near zero, so that the fractional Kelvin-Voigt material included damping effects (refer to Remark 4).In particular, in Figure 2, we showed that the increase of the elastic constituent's volumetric fraction diminishes the damping effect created by the existence of the fractional Kelvin-Voigt material.We also showed that for decreasing values of the fractional Kelvin-Voigt constituent's volumetric fraction, the oscillatory behavior of the effective coefficient decreases and tends, as expected, to the effective coefficient of an equivalent elastic composite.In the second benchmark problem, we examined the case of a composite made of a fractional Maxwell constituent and a fractional Zener one.In this case, following previous works of some of the authors, we numerically computed the effective coefficients and illustrated the influence of different combinations of    and    .As it was for the first benchmark problem, we noticed that the instant of time at which the effective properties remain constant depends not only on the values of the fractional parameters but also on the constituents' volumetric fractions.
As we mentioned in Section 1, there is still lack of theoretical and experimental studies connecting the role of the micro-structure in the origins and description of the fractional behavior of some kinds of materials.Thus, part of our current investigations focuses on modeling the underlying mechanisms that give rise to fractional behavior starting from the micro-structure.On the other hand, we will explore the consideration of micro-structural information in the fractional parameters and study specific scenarios where this information could be valuable.
and Y 1 ∩ Y 2 = ∅, where Y 1 and Y 2 are open subsets.Furthermore, we denote by I Y the interface between Y 1 and Y 2 and with n Y the normal vector to I Y where the positive orientation is taken from Y 1 to Y 2 , that is, n Y is directed outwards the interior of Y 1 .

TABLE 1
Three types of viscoelastic models.

TABLE 2
Three types of fractional viscoelastic models.

TABLE 3
Material properties.
Comparison of the numerical and analytical values of [ eff ] 3333 for different values of .[Colour figure can be at wileyonlinelibrary.com].