Effective behavior of long and short fiber-reinforced viscoelastic composites

We study the homogenized properties of linear viscoelastic composite materials in three dimensions. The composites are assumed to be constituted by a non-aging, isotropic viscoelastic matrix reinforced by square or hexagonal arrangements of elastic transversely isotropic long and short fibers, the latter being cylindrical inclusions. The effective properties of these kind of materials are obtained by means of a semi-analytical approach combining the Asymptotic Homogenization Method (AHM) with numerical computations performed by Finite Elements (FE) simulations. We consider the elastic-viscoelastic correspondence principle and we derive the associated local and homogenized problems, and the effective coefficients in the Laplace–Carson domain. The effective coefficients are computed from the microscale local problems, which are equipped with appropriate interface loads arising from the discontinuities of the material properties between the constituents, for different fibers’ orientations in the time domain by inverting the Laplace–Carson transform. We compare our results with those given by the Locally Exact Homogenization Theory (LEHT), and with experimental measurements for long fibers. In doing this, we take into consideration Burger’s and power-law viscoelastic models. Additionally, we present our findings for short fiber reinforced composites which demonstrates the potential of our fully three dimensional approach.


Introduction
Materials characterized by both a viscoelastic response and a composite-like geometrical arrangement are found in several biological contexts driven by natural evolution, see, e.g., ( Atthapreyangkul et al., 2021;Ojanen et al., 2017;Sherman et al., 2017 ). Especially, long and short fiber-reinforced composites are being increasingly exploited in a variety of engineering and manufacturing processes because of their capability of optimising properties such as light weight, stiffness, and strength. On the one hand, high performance composites are typically made of long continuous fibres embedded in a polymer matrix and exhibit viscoelastic properties (see, e.g., Ornaghi Jr. et al., 2020;Wang et al., 2020 ). On the other hand, reinforcement via short fibers (i.e., fiber-shaped inclusions) can provide a valuable alternative in the modelling of failures appearing in composites reinforced by long fibers, see, e.g., Cepero-Mejías et al., 2020;Nonato Da Silva et al., 2020 . Material composites reinforced by short fibers can also provide significant economical and manufacturing advantages over continuous fiberreinforced composites without compromising high performance, as long as the aspect ratio is high enough to support load transfer ( Huang and Huang, 2020;Wang and Smith, 2019;Yu et al., 2014 ).
The multiscale modelling of viscoelastic composites has been increasingly addressed in recent contributions. In this respect, micromechanical models are particularly suitable whenever the aim is to determine the effective response of materials on the basis of individual constituents' properties, such as viscoelastic moduli and fibers properties in terms of geometrical arrangement, volume fraction and orientation. For instance, in Sevostianov et al. (2016) , the effective viscoelastic properties of short fiber reinforced composites are investigated by means of the fraction-exponential operators of Scott Blair-Rabotnov. Moreover, in Kern et al. (2019) a frequency-domain finite element simulations are considered to determine the effective moduli of viscoelastic coated fiberreinforced composites. The investigation of the polymer alignment with the aid of direct numerical simulations of the turbulent channel flow of a viscoelastic FENE-P fluid is conducted in Pereira et al. (2020) . In addition, the effective viscoelastic creep behavior of aligned short fiber composites is obtained in Wang and Smith (2019) via a RVE-based Finite Element algorithm. In Ornaghi Jr. et al. (2020) , the Authors eval-uate the creep, recovery, and viscoelastic properties of unidirectional carbon/epoxy filament wound composite laminates under controlled stress, time, and temperature. Recently, a probabilistic micromechanics damage framework to predict the macroscopic stress-strain response and progressive damage in unidirectional glass-reinforced thermoplastic polymer composites has been proposed in Chen et al. (2021) .
Among the most used techniques addressing the calculation of the effective properties of viscoelastic heterogeneous structures we find the Asymptotic Homogenization Method (AHM). For example, analytical closed form expressions for the effective coefficients of fibrous viscoelastic composites are obtained in Rodríguez-Ramos et al. (2020) , Otero et al. (2020) by means of the two-scale AHM. The theoretical bases of the method are found in the contributions of several authors (see, e.g., Allaire and Briane, 1996;Auriault et al., 2009;Bakhvalov and Panasenko, 1989;Bensoussan et al., 1978;Sanchez-Palencia, 1980 ). In general, the technique takes advantage of the scales separation assumption for decoupling the spatial variables into a microscopic and a macroscopic one. Thus, the solution of the original problem at the macroscale is approximated by considering the solution of the corresponding local problems on a periodic cell and of the homogenized problem. This permits to decrease the computational complexity of the problem at hand, and to encode the information at the microscale into the effective coefficients of the macroscale model.
The main disadvantage of the AHM is that the analytical solution of the local problems can be derived for a few of simple composite structures which can be reduced to one or two dimensions analysis (as those for laminated composites and long fibers). For instance, in studies related with viscoelastic composites for laminated and fiber reinforced composites ( Cruz-González et al., 2020a;Otero et al., 2020 ). For this reason, in order to handle more complex microstructures, numerical approaches based on the Finite Elements (FE) provide a robust alternative to solve the local problems for more general microstructures. Based on these considerations, in the present work, we aim to study the effective properties of viscoelastic composites by means of the combination of the AHM and the FE method, the latter allows us to find the numerical solution of the microscale periodic local problem for different three-dimensional arrangements. This approach provides a new and efficient computational platform for computing the effective properties of viscoelastic composites in three dimensions.
The main aim of this work is to calculate the effective properties of non-aging, linear viscoelastic composites in three dimensions. For our purposes, we employ the modeling approach introduced in our previous work ( Cruz-González et al., 2020b ), wherein a combined framework based on theoretical and computational techniques for computing the effective properties of viscoelastic composites is employed. Specifically, in Cruz-González et al. (2020b) , we analyzed several types of composite structures reinforced by fibers and inclusions. Here, however, we go further in our investigations and consider hexagonal periodic cells and unidirectionally aligned short fibers. Furthermore, in the present framework, we consider the fibers to be transversely isotropic, while in Cruz-González et al. (2020b) , the study was focused on composites with isotropic constituents. This extension requires a suitable generalisation of the computational set up of the problem found, for example, in Penta and Gerisch (2016) for elastic composites. In particular, the interface loads related to the auxiliary local cell problems are obtained for more general orthotropic materials, and then specialised to transverse isotropic constituents in our calculations. Another extension of the present work with respect to Cruz-González et al. (2020b) is that, in the present study, we take into consideration different fiber's orientations.
We further notice that, to the best of our understanding, the major novelty of this work with respect to others existent in the literature and focused on the study of viscoelastic composites (see, e.g. Amiri-Rad et al., 2019;Ornaghi Jr. et al., 2020;Otero et al., 2020;Rodríguez-Ramos et al., 2020;Sevostianov et al., 2016;Tang and Felicelli, 2015;Tran et al., 2011;Wang and Pindera, 2016a;Wang and Smith, 2019;Yancey and Pindera, 1990;Yi et al., 1998 ) is that here,  (i,ii) or hexagonal (iii,iv) arrays of non-overlapping long and short fibers, respectively. (c) -structural level. Microscale: periodic cell for long and short fibers inclusions that do not intersect the boundaries.
using a semi-analytical approach, we are capable to compute the effective properties of three-dimensional, non-aging, viscoelastic composite materials. Particularly, it should be pointed out that although we considered some of the results given in Yancey and Pindera (1990) and Wang and Pindera (2016a,b) for comparison with ours, our methodology provides further developments because of the possibility of considering more complex cell geometries. Indeed, we calculate the effective properties for a viscoelastic composite reinforced with perfectly aligned short fibers and this geometrical configuration is not treatable under the two-dimensional formulation reported in Wang and Pindera (2016a,b) .
The manuscript is organized as follows. In Section 2 , we present the geometrical description of the model and we formulate the linear viscoelastic heterogeneous problem. In Section 3 , we apply the twoscale AHM to obtain the macroscale functional form of the effective viscoelastic coefficients. In Section 4 , we calculate the effective properties by solving appropriate local (cell) problems in three dimensions. The theoretical framework is illustrated in general by considering that both phases are viscoelastic, although the results are then presented by considering purely elastic fibers embedded in a viscoelastic matrix for the sake of comparison against previous analytical results and relevant experiments. In Section 5 , we compare our results against alternative homogenisation techniques in the case of long fibers, and emphasise the potential of our new approach by illustrating the results in the case of composites reinforced by short fibers (i.e., cylindrical inclusions). Finally, in Section 6 , we summarize our findings and highlight the limitations of the current model and possible further developments of the work.

Model description
We identify the heterogeneous, linear viscoelastic material with an open, bounded set Ω ⊂ ℝ 3 (see Fig. 1 (a)). In particular, we consider Ω as a two-constituent composite made of a matrix reinforced by square ( Fig. 1 (b) (i,ii)) or hexagonal ( Fig. 1 (b) (iii, iv)) arrays of unidirectional and periodically distributed long and short fibers in Ω (see Fig. 1 (b)). Furthermore, we consider the existence of two distinct, well-separated length scales and which are related with the characteristic size of the periodic micro-structure and that of the whole composite, respectively (see Fig. 1 ). In this framework, we introduce the dimensionless scaling parameter as follows, and the microscopic spatial variable where is said to be the macroscopic spatial variable.
In particular, we set Ω = Ω 1 ∪ Ω 2 with Ω 1 ∩ Ω 2 = Ω 1 ∩ Ω 2 = ∅, and where Ω 2 denotes the matrix and Ω 1 = ∪ =1 Ω 1 represent the inclusions with ∈ ℕ . Additionally, the interface between Ω 1 and Ω 2 (see Fig. 1 (b)) is denoted by Γ . Referring to Fig. 1 (c), the unitary periodic cell is considered to be constituted by a fiber (long or short) 1 in the matrix 2 so that the periodic cell is given by = 1 ∪ 2 with 1 ∩ 2 = 1 ∩ 2 = ∅. At this scale, the interface between 1 and 2 is denoted by Γ , see also Di Stefano et al. (2020) for an illustration of various periodic cell arrangement in the context of electro-active composites.

Statement of the problem
For the sake of simplicity, we neglect inertia and external volume forces in the model, and impose continuity conditions for displacements and tractions on the interface Γ , i.e., the matrix and the sub-phases are in ideal contact. Therefore, the balance of linear momentum in Ω ⧵ Γ together with the interface conditions read where represents the second-order stress tensor and is the displacement field. Moreover, ( ) denotes the outward unit vector to the interface Γ , and the operator describes the jump of across the interface Γ between Ω 1 and Ω 2 . Notice that the superscript is used to indicate the notation ( , ) = ( , , ) (refer to Di Stefano et al., 2020 for a comprehensive discussion regarding this notation). The system of Eqs. (3a) -(3c) has to be supplemented with boundary conditions on Ω × ℝ and initial conditions in Ω × {0} . However, these conditions do not play a role in the derivation of the effective coefficients, and they are typically to be specified explicitly only when the aim is to obtain a specific solution of the macroscale system of homogenized PDEs, which is not the case here.
In the present framework, the composite behaves as a non-aging viscoelastic material so that ( Christensen, 1982 ) where  is the fourth-order tensor of relaxation moduli, which here is scale-dependent. The discontinuities of the properties between the host medium and the subphases are encoded in the tensor  and, thus, its components are assumed to be smooth functions of in (Ω ⧵ Γ ) × ℝ , but discontinuous on Γ × ℝ . Furthermore, we notice that, in Eq. (4) , denotes the second-order strain tensor for small displacements, namely and we require both minor and major symmetry properties for  , i.e., The integral in Eq. (4) , standing for the stress-strain relationship for non-aging, viscoelastic materials, can be manipulated by means of integral transforms. In particular, the Laplace-Carson transform, which is given by where is the variable in the Laplace-Carson space, reduces (4) to an algebraic equation representing the constitutive relations in classical elasticity theory (see, for instance, Lakes, 2009 ). This methodology originally proposed by Hashin (1965) and known as the elastic-viscoelastic correspondence principle, continues to gain interest in the scientific literature (see Liu et al., 2020;Vilchevskaya et al., 2019;Yang et al., 2019 and references therein). Hence, based on the above considerations, the original system (3a) -(3c) written in the Laplace-Carson domain is given by

Asymptotic homogenization approach
In this section, we summarize the methodology described in Cruz-González et al. (2020b) and write the local and homogenized problems resulting from the application of the AHM to Eqs. (8a) -(8c) .
Before going further, we remark that using the chain rule the following relation for the spatial derivatives holds Moreover, Eq. (5) becomes, where we have introduced the notation The AHM ( Bakhvalov and Panasenko, 1989;Cioranescu and Donato, 1999 ) proposes the solution of the viscoelastic heterogeneous problem (8a) -(8c) as a formal series expansion in powers of . In the Laplace-Carson domain it reads where the coefficients ̂ ( ) ( , , ) are assumed to be periodic in the microscopic variable . Thus, following the standard procedure in asymptotic homogenization (see Bakhvalov and Panasenko, 1989;Cioranescu and Donato, 1999 ), after substitution of the series expansion (12) in the original problem (8a) -(8c) and by equating the result in the same powers of , we obtain that, in the limit → 0 , where ̂ is a smooth vector function of and , and the third order tensor ̂ is the solution of the -local problem given by The uniqueness of the solution of the local problem (14a) -(14d) is guaranteed by enforcing either, the condition ⟨̂ ⟩ = 0 or by fixing the value of ̂ at one point of the reference periodic cell (see Penta and Gerisch, 2016;Penta and Gerisch, 2017 ). In particular, in the subsequent sections we will use the latter because of its advantage when dealing with numerical simulations. Note that, the -local problem has to be supplemented with an initial condition in × {0} .
For completeness in our analysis, we report the homogenized equation at the macroscale in the Laplace-Carson space, which is obtained after equating in the same powers of 0 . Specifically, this can be written as where is the effective relaxation modulus in the Laplace-Carson space. In Eq.
(17) , the notation ⟨ ⟩ denotes the cell average operator and is defined by the expression with | | being the volume of the periodic cell .

Calculation of the effective properties
For simplicity in our calculations, we consider that the relaxation modulus,  , is -constant in each constituent of the periodic cell , i.e.
where the superscript indicate the corresponding constituent, "(1) " for the matrix and "(2) " for the inclusion (see Fig. (1) (c)). Then, the local problem (14a) -(14d) is rewritten as follows, We notice that, in Eq. (20d) , the stress jump conditions lead to the occurrence of interface loads, i.e., which arise as a consequence of the discontinuities in the coefficients of  between the host medium and the sub-phases, and represent the driving force to obtain nontrivial solutions of the six elastic-type local problems ( ( , ) , ≥ ) (see Penta and Gerisch, 2016;Penta and Gerisch, 2017 ). In particular, when the matrix and subphases are orthotropic materials and considering Voigt's notation, the interface loads ( ) read ( ) where { } 3 =1 represents the standard vector basis.

Numerical approach
At this point, it is possible to solve numerically the set of elastictype local problems (20a) -(20e) in the Laplace-Carson space and then, to compute the effective viscoelastic properties by using (17) . For this purpose, we use the finite element software COMSOL Multiphysics®in conjunction with LiveLink TM for Matlab ® scripting (see Cruz-González et al., 2020b;Penta and Gerisch, 2016;Penta and Gerisch, 2017 ). Particularly, once  ( * ) is known, the effective creep compliance  ( * ) in the Laplace-Carson space can be computed through the relationship where denotes the components of the fourth-order identity tensor (see Hashin, 1972 ).
The inversion of the effective coefficients to the original time domain is performed by employing the MATLAB's function INVLAP (see Juraj, 2020;Valsa and Bran ĉik, 1990 and referred here to as Valsa's method. The steps are summarized as follows, where Re indicates the real part of a complex variable. Moreover, to take into account the different orientations that the unidirectional viscoelastic composites may have, we rotate the effective viscoelastic tensors  ( * ) and  ( * ) by an angle and obtain  ( * ) and  ( * ) . In this respect, the following transformations are useful, where and  ( , = 1 , 2 , 3) are the coefficients of the orthogonal rotation tensor (see Ting, 1996 ).
To conclude with this section, it is worth to remark that steps (a)-(d) in the inversion process are equivalents to the stage (IV) in Cruz-González et al. (2020b) . Here, we illustrate through a flowchart (see Fig. 2 ) the methodology described in (I)-(IV) of Cruz-González et al. (2020b) with more details. In the following sections, we refer to as AHMFE the semi-analytical approach proposed in the present work, which combines the Asymptotic Homogenisation Method (AHM) and Finite Elements (FE) simulations.

Instant elastic response
To begin with our analysis, in this section, we compute the effective instant elastic response of a composite with a hexagonal arrangement of transversely, purely elastic isotropic long fibers (see Fig. 1 (b)-(iii)). Although the theoretical framework introduced in the previous section holds for viscoelastic fibers as well, we focus on purely elastic fibers for the sake of comparison of our results with alternative analytic tech-  niques. In this context, the elastic limit case is reached by considering = 0 in Eq. (4) , which implies that  ( * ) becomes the effective stiffness tensor.
Here, we determine the instant elastic effective response of a graphite/epoxy system with hexagonal architecture and different fiber volume fractions. It is worth noticing that the following results differ from the ones obtained in Cruz-González et al. (2020b) since therein it was considered a square array of inclusions in the matrix phase comprising isotropic constituents. The values of the parameters reported in Table 1 are obtained from Wang and Pindera (2016b) . The notation ( ) and ( ) is used for the axial (transverse) Young's and shear moduli, respectively, and represents the axial Poisson's ratio. Before we proceed with the results of the effective coefficients, we gather information on two main features related with the computational approach. The solution's convergence behavior, through three types of mesh discretization, and the execution times required for these calculations. With this purpose, we only provide the results of the computation of the effective transverse Young's modulus ( * ) since those related with the transverse Poisson's ratio ( * ) and the axial shear modulus ( * ) show similar characteristics. Furthermore, with regards to the analysis of the solution's convergence, we use the meshes A , B and C reported in Fig. 3 .
In Fig. 4 , we show the effective transverse Young's modulus ( * ) , normalized by the corresponding matrix Young's modulus . Specifically, in Fig. 4 (a) the effective transverse Young's modulus is computed, for each of the meshes discretizations A , B and C , as a function of the volume fraction . Moreover, in Fig. 4 (b) we compare the results obtained with the present model (AHMFE) with those produced by Wang and Pindera (2016b) (see Fig. 7 in Wang and Pindera, 2016b ) using the finite-volume direct averaging micromechanics (FVDAM) theory. A closer look at the data in Fig. 4 (b) reveals that the relative error between the solutions the present approach (AHMFE) and that in Wang and Pindera (2016b) (FVDAM), namely reaches his maximum value when the coarser mesh (Mesh A ) is considered and is less than 1 . 3% . On the other hand, the relative error computed with the meshes discretizations B and C is less than 0 . 2% . Then, we can conclude that our simulations provide a good agreement with FV-DAM, and that, in this case, our results do not undergo large variations after subsequent mesh discretizations.
To continue with our analysis, in Fig. 4 (c) we provide the total execution times needed to obtain the entire set of effective moduli in relation to the three mesh discretizations. The simulations are set up to take into account a finite number of fiber volume fractions ranging from = 0 . 05 to = 0 . 7 with an increment of 0.05, and they were performed using a machine running Windows 10 Professional 64-bit operating system, with 32.0 GB RAM and Intel(R) Core(TM) i5-8350U CPU 1.70GHz. As illustrated in Fig. 4 (c), the performance using the mesh C could be considered inefficient in terms of time, in part due to the fact that the mesh C is significantly finer compared to the meshes A and B . Particularly, charts (d), (e) and (f) in Fig. 4 offer more precise details on the computing time for each volume fraction and each mesh. The differences in the results are due to the fact that different volumetric fractions modify the geometry of the periodic cell for the corresponding fixed mesh. We also highlight the contrast in relation to the mean computing time. Taking into account both the relative errors and the execution times, we conclude that the mesh B is the best possible choice for our simulations.
Finally, to complete our analysis, in Fig. 5 , we show the numerical values of the effective moduli ( * ) ∕ , ( * ) and ( * ) ∕ for a hexagonal array of unidirectional long fibers and different contrasts in the constituents. In particular, we study the pairs glass/epoxy, boron/aluminum, aluminum/porosity and graphite/epoxy. We refer to Table 1 for the material properties of these constituents. By referring to Table 2 , in which summarize maximum relative errors, we note that our numerical results are in agreement with the results provided in Wang and Pindera (2016b) using the locally exact homogenization theory (LEHT). As observed in Table 2 the maximum relative error is less than 1 . 46% .

Viscoelastic response
In this section, we compute the effective properties of linear viscoelastic composite materials. Particularly, in Section 5.2.1 we compute the effective properties for a composite material made of isotropic constituents where the viscoelastic behaviour of the matrix is described by means of a Burger's model. We consider a long fiber reinforcement with square and hexagonal geometrical arrangements, and consider the case in which Poisson's ratio or the bulk modulus of the viscoelastic matrix are constants. Furthermore, in Section 5.2.2 , we deal with transverse isotropic long fibers with different orientations and consider the power-law model given in Yancey and Pindera (1990) for the characterization of the creep compliance of the viscoelastic matrix. In doing this, we obtained good agreements with both the LETH and experimental results. Finally, in Section 5.2.3 , we address the calculation of the effective properties for a composite made of perfectly aligned short fibers embedded into a viscoelastic matrix with transversely isotropic behaviour. So, we show the potential of our approach in the solution of fully three-dimensional problems involving inclusions, which cannot be addressed by means of analytical methods such as the LEHT. We mention that in the upcoming simulations, we employ Mesh B .

Comparison with locally-exact homogenization theory
In this section, we analyze the influence of the square and hexagonal arrangement of long fibers (see Fig. 1 ) in the calculation of the effective creep compliance in polymeric matrix composites. In addition, as in Wang and Pindera (2016a) , we consider that either Poisson's ratio ( ) or the bulk modulus ( ) of the viscoelastic matrix is constant. This assumption is justified by the necessity of modeling the potentially time-independent response of polymeric matrices under hydrostatic loading Wang and Pindera (2016a) . On the one hand, we take = 0 . 38 as the value for the matrix and therefore the time-dependent bulk modulus is given by ( ) = ( )∕(3(1 − 2 )) , where ( ) stands for the one-dimensional relaxation modulus. On the other hand, if we assume a constant bulk modulus for the matrix, we use the equation = 0 ∕(3(1 − 2 )) to determine its value, and then, the time-dependent Poisson's ratio ( ) arises from the equation ( ) = 1∕2 − ( )∕(6 ) . Here, 0 represents the instantaneous elastic relaxation modulus.
In the following, we investigate the effective properties of unidirectionally-reinforced glass/epoxy composites with linear isotropic constituents, where elastic glass fibers are embedded into a viscoelastic polymeric matrix (see, e.g. Cavalcante and Marques, 2014;Chen et al., 2017;Cruz-González et al., 2020b;Wang and Pindera, 2016a ). Here, in contrast with Cavalcante and Marques, 2014;Chen et al., 2017;Cruz-González et al., 2020b , we consider hexagonal periodic cells in the computations. The mechanical properties of the fibers are given by Young's modulus = 68 . 77 GPa and Poisson's ratio = 0 . 21 . Furthermore, we describe the viscoelastic matrix by assuming the relaxation representation of the four-parameter model or Burger's model, i.e. two Maxwell elements set in parallel (see Mainardi and Spada (2011) for further details). Specifically, the expression of the relaxation modulus is given as follows, where the material properties are taken from Wang and Pindera (2016a) by means of the scalar form of (23) (see Park and Kim, 1999 ) and some algebraic transformations. The data set is reported as follows, 1 = 1 . 12511 GPa , 2 = 2 . 14489 GPa , , 1 = 6999 . 34 h and , 2 = 58 . 2551 h , where ( = 1 , 2 ) represents the elastic modulus of the spring and , is a relaxation time Mainardi and Spada (2011) . It is worth to remark that from eq. (27) , we obtain 0 = 1 + 2 .
In Fig. 6 , we show the curves corresponding to the effective creep compliances  ( * ) 22 ,  ( * ) 44 and  ( * ) 66 for constant Poissons ratio (left charts) and constant bulk modulus (right charts). Specifically, in Fig. 6 , we compare our results with those obtained in Wang and Pindera (2016a) via the LEHT. In these comparisons, hexagonal and square arrays of fibers are studied for a fiber volume fraction equal to 0.6. As it can be noticed, there is a good agreement between the two approaches, which is further evidenced by the maximum relative errors provided in Table 3 .

Power-law model. Comparison with experiments
In this section, we follow the analysis adopted in Wang and Pindera (2016a) , and compare our results with the experimental creep measurements obtained in Yancey and Pindera (1990) for off-axis graphite/epoxy specimens, and with the numerical results given in Table 4 Elastic properties of the transversely isotropic T300 graphite fiber at room ( 22 • ) and elevated ( 121 • ) temperature.  Wang and Pindera (2016a) . In Yancey and Pindera (1990) , the authors observed that at 22 • and 121 • the T300 graphite fiber present an elastic behavior, whereas the creep response of the 934 epoxy matrix is fitted by the power-law where and are experimentally measured parameters and 0 is the instantaneous elastic relaxation modulus ( Yancey and Pindera, 1990 ). Before proceeding, it is worth mentioning that, even though in the present model we do not deal with fractional viscoelasticity (we refer the Reader to Atanackovi ć et al., 2016;Beltempo et al., 2019;Bouras et al., 2018;Mainardi, 2010;Mainardi and Spada, 2011 and the references therein), we find convenient to use some of the results given in these works for the calculation of the relaxation modulus ̂ ( ) which is needed in our simulations.
With this purpose, we introduce the notation = 0 ∕2 , = −1∕ , = 0 Γ(1 + ) and = , so that Eq. (28) can be equivalently rewritten as where Γ denotes the Gamma function, ∈ ]0 , 1] , and ( ) represents the fractional creep compliance ( Mainardi and Spada, 2011 ). Hence, by referring to the results obtained in Mainardi and Spada (2011) , the fractional relaxation modulus is given by the expression where  denotes the Mittag-Leffler function of order , which is defined by the expression ( Gorenflo et al., 2014 ) and for = 1 reduces to  ((− ∕ ) ) = exp (−t∕τ) . We notice that, from Eq. (30) and by virtue of the parameter identifications established for obtaining Eq. (29) , obtain However, Eq. (32) does not provide, in a direct way, an incremental solution for the problem in the time domain. Therefore, we rely on the expression for the relaxation modulus ̂ ( ) , expressed with respect to the Laplace-Carson domain, for the computation of the effective properties. Thus, by employing the results given in Mainardi and Spada (2011) , the Laplace-Carson transform of Eq. (30) is which implies that In this way, the Laplace-Carson transform for the relaxation modulus can be analytically obtained in an explicit form, which highly reduces the numerical complexity of the problem. In the following, we consider a composite with hexagonal arrangement of long fibers where the properties of the constituents, i.e. the elastic fibers and the viscoelastic matrix given in Wang and Pindera (2016a) are summarized in Tables 4 and 5 , respectively.   Similarly to the previous section, we consider the cases of constant Poisson's ratio and constant bulk modulus. In addition, we analyze the influence of two dif-ferent temperatures, i.e. 22 • (room) and 121 • (elevated). The results correspond to composites with fibers that are rotated 10 • and 90 • counterclockwise about the 2 -axis. In this case, the fiber volume fraction of the fiber is fixed to = 0 . 6 . A comparison of our results with those obtained in Wang and Pindera (2016a) using LEHT shows a good agreement between both approaches. Additionally, the qualitative behavior of the curves is very similar to the experimental data. In Table 6 , we provide the maximum relative errors between the results obtained via two methods.

Modeling of short fiber reinforcement
The locally-exact homogenization theory (LEHT) is based on a twodimensional formulation which is only capable of taking into account  long cylindrical fibers (see, e.g. Chen et al., 2017 ). In this section, we show the potential of the AHMFE approach in the modeling of viscoelastic composites for three-dimensional geometrical configurations. Particularly, we consider a viscoelastic composite material with square and hexagonal arrangement of perfectly aligned short fibers (see Fig. 1 (b) (ii,iv)) represented by cylindrical inclusions. In addition, we assume that the constituents behave as those in Section 5.2.2 , that is, we consider a viscoelastic matrix (934 epoxy) with power-law creep compliance as given in (28) , and reinforced by transversely isotropic elastic fibers (T300 graphite). In this context we define the parameters 1 ∶= ℎ 1 ∕ 1 and 2 ∶= ℎ 2 ∕ 2 which relate the length measurement of the fiber and the matrix in the square and hexagonal periodic cell, respectively (see Fig. 8  (a)). In addition, we assume the fibers to be centered in the periodic cells. Therefore, we notice that 0 ≤ 1 , 2 ≤ 1 , where a zero value represents a homogeneous material made only with the matrix and a one value reproduces the case of long fibers as particular case of this approach. Moreover, as observed in Fig. 8 (a), we study counterclockwise uniform rotations of the short fibers about 2 -axis of 0 ≤ < .
In Fig. 8 (b) and (c), we show our findings in the calculation of the effective moduli ( * ) 1 and ( * ) 32 for square and hexagonal periodic cell. In particular, we assume the room temperature ( 22 • ) and the constant bulk modulus approach discussed in the Section 5.2.2 . The fiber volume fraction and the ratio are fixed to = 0 . 1 and = 1 = 2 = 3∕5 , respectively. In the simulations, we consider the rotation angles ∈ {0 , ∕2} .

Discussion and conclusions
In this work, we studied the effective properties of threedimensional, non-aging viscoelastic composites reinforced with square and hexagonal arrangements of rotated elastic short and long fibers. We exploited the correspondence principle by transforming the original viscoelastic problem into a "fictitious " elastic one. In particular, by means of the asymptotic homogenization technique, we studied the two-scale problem as a homogenized one in which the information available at the smallest scale is encoded into the effective coefficients. The latter were calculated through the conception of dedicated numerical algorithms combining finite elements with Valsa's inversion method ( Valsa and Bran ĉik, 1990 ).
In order to validate our results, we considered the case of instant elastic response. We studied the solution's convergence behavior and the execution time by means of three different mesh discretizations. Specifically, in our simulations, we took into account transverse isotropic long fibers in a hexagonal arrangement. In addition, we reached good agreements in the comparisons against FVDAM and LEHT, and, concerning the mesh performance, we concluded that Mesh B was the best choice for our computations, see Fig. 3 .
We also studied the effective properties of viscoelastic composite materials by assuming that either Poisson's ratio or the bulk modulus of the viscoelastic matrix can be constants. In this context, we considered the Burger's model with isotropic constituents and compared our numerical results with those obtained in Wang and Pindera (2016a) when considering the locally-exact homogenization theory. We found good agreements in the comparisons, with a maximum relative error less than 2 . 3% (see Table 3 ). We further compare our technique, using the power-law representation given in Yancey and Pindera (1990) , with the experimental results in Yancey and Pindera (1990) . We obtained the relaxation representation of the power-law model in the time and Laplace-Carson domains by considering some of the results given in the context of fractional viscoelasticity ( Atanackovi ć et al., 2016;Beltempo et al., 2017;Bouras et al., 2018;Mainardi, 2010;Mainardi and Spada, 2011 ). Our numerical results agreed with the experimental measurements and those in Wang and Pindera (2016a) using the LEHT. In particular, we obtain maximum relative errors less than 6 . 8% and 0 . 4% , respectively.
Finally, we addressed the effective properties of viscoelastic composites reinforced by perfectly aligned, transversely isotropic short fibers. This is one of the novelties of this work. In fact, the locally-exact homogenization approach adopted in Wang and Pindera (2016a) and other analytical methods (for instance, using the asymptotic homogenization method addressed in Rodríguez-Ramos et al. (2020) ) cannot consider this case. Specifically, they rely on long fibers (which extend from the top to the bottom of the cell) to obtain reduced problems formulated in two-dimensions, which can therefore be solved analytically. For the simulations, we introduced the parameter to manage the size of the fiber inside the periodic cell and we assumed different fiber orientations. Our findings show that the orientation of the fiber at the microscale has a direct influence in the effective behaviour of the composite at the macroscale. As observed in the parametric study conducted for in Fig. 8 (b) and (c), when the rotation increases, we obtain that the values of the effective moduli ( * ) 1 and ( * ) 32 decrease. These analyses have a positive impact on the manufacturing processes and highlight the advantages of micro-mechanical models.
We mention that one of the main limitations of this work lies in the computation times. The forward-back passage between the time domain and the Laplace-Carson domain compromises the efficiency of the model. In fact, in order to invert the effective coefficients via Valsa's method (see Juraj, 2020;Valsa and Bran ĉik, 1990 ), we needed to calculate + + 1 = 40 points in the Laplace-Carson space to determine one point in the time domain. One possible way to deal with this issue is to perform parallel computations. In addition, we could adopt different inversion techniques such as Zakian's method ( Zakian, 1969 ) which reduces the number of points from 40 to 5 and have been used in several works (see, e.g., Chen et al., 2017;Wang and Pindera, 2016a ).
Further generalizations of the present framework include the consideration of an orientation distribution for the short fibers ( Mishurova et al., 2018;Sevostianov et al., 2016 ). Moreover, we can generalize our approach through the consideration of imperfect contact conditions between the interface of the constituents (see Daridon et al., 2016;Escarpini Filho and Marques, 2016;Rizzoni et al., 2021;Serpilli et al., 2019 ). Besides, another development is related to the modeling of three-dimensional viscoelastic composites materials by means of a three scale asymptotic homogenization approach ( Cruz-González et al., 2020a;Ramírez-Torres et al., 2019 ). In this context, the current methodology can address the modeling of short and long fiber reinforcement in hierarchical composite materials. Another challenging extension consists in generalizing the present framework in order to take into account soft viscoelastic composites characterised by a nonlinear response (by potentially extending the framework recently developed in  and typically found in relevant biological and industrial applications, see, e.g. Destrade and Saccomandi (2007) ; Muliana and Rajagopal (2012) ; Parnell and De Pascalis (2019) . Finally, we would like to mention that, in view of the advantages that the theory of Fractional Calculus ( Atanackovi ć et al., 2014 ) offer. A natural further generalization of the present work is to consider and adapt the ideas put forward in Ramírez-Torres et al. (2021b) and in Ramírez-Torres et al. (2021a) , where problems of fractional diffusion at different spatial scales are studied.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.