Effective behavior of long and short fiber-reinforced
viscoelastic composites

Oscar Luis Cruz González, A. Ramírez-Torres, R. Rodríguez-Ramos, J.A.

Otero, R. Penta, F. Lebon

► To cite this version:

Oscar Luis Cruz González, A. Ramírez-Torres, R. Rodríguez-Ramos, J.A. Otero, R. Penta, et al. Effective behavior of long and short fiber-reinforced viscoelastic composites. Applications in Engineering Science, 2021, 6, pp.100037. 10.1016/j.apples.2021.100037. hal-03183369

HAL Id: hal-03183369

https://hal.science/hal-03183369

Submitted on 27 Mar 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Effective behavior of long and short fiber-reinforced viscoelastic composites

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

Check for
updates

O.L. Cruz-González<sup>a</sup>, A. Ramírez-Torres<sup>b</sup>, R. Rodríguez-Ramos<sup>c</sup>, J.A. Otero<sup>d</sup>, R. Penta<sup>b,*</sup>,
F. Lebon<sup>a</sup>

## ARTICLE INFO

**Keywords:**

Homogenization

Viscoelasticity

Fiber reinforced composites

Power-law model

Transverse isotropy

Finite elements

## ABSTRACT

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.

## 1. 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 fiber-reinforced 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 fiber-reinforced 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-

<footnote><sup>a</sup>Aix Marseille Univ, CNRS, Centrale Marseille, LMA UMR 7031, Marseille, France</footnote>

<footnote><sup>b</sup>School of Mathematics and Statistics, Mathematics and Statistics Building, University of Glasgow, University Place, Glasgow G12 8QQ, UK</footnote>

<footnote><sup>c</sup>Facultad de Matemática y Computación, Universidad de La Habana, San Lázaro y L, Vedado, La Habana CP 10400, Cuba</footnote>

<footnote><sup>d</sup>Escuela de Ingeniería y Ciencias, Tecnológico de Monterrey, Campus Estado de México, Atizapán de Zaragoza EM CP 52926, Mexico</footnote>

<footnote>* Corresponding author.
E-mail address: Raimondo.Penta@glasgow.ac.uk (R. Penta).</footnote>
uate the creep, recovery, and viscoelastic properties of unidirectional carbon/epoxy filament wound composite laminates under controlled stress, time, and temperature. Recently, a probabilistic micromechan- ics 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 viscoelas- tic 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 assump- tion for decoupling the spatial variables into a microscopic and a macro- scopic 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 per- mits to decrease the computational complexity of the problem at hand, and to encode the information at the microscale into the effective coef- ficients 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 struc- tures which can be reduced to one or two dimensions analysis (as those for laminated composites and long fibers). For instance, in studies re- lated 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 ap- proaches 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 effec- tive properties of viscoelastic composites by means of the combination of the AHM and the FE method, the latter allows us to find the nu- merical solution of the microscale periodic local problem for different three-dimensional arrangements. This approach provides a new and ef- ficient 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 compos- ite 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 generali- sation 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 trans- verse 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 ma- jor novelty of this work with respect to others existent in the liter- ature 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 Fe- licelli, 2015; Tran et al., 2011; Wang and Pindera, 2016a; Wang and Smith, 2019; Yancey and Pindera, 1990; Yi et al., 1998) is that here,

<figure><img src="image_4.png" /><figcaption>Fig. 1. (a) Macroscale: viscoelastic heterogeneous material with (b) square (i,ii) or hexagonal (iii,iv) arrays of non-overlapping long and short fibers, respec- tively. (c) ε-structural level. Microscale: periodic cell for long and short fibers inclusions that do not intersect the boundaries.</figcaption></figure>

using a semi-analytical approach, we are capable to compute the effec- tive properties of three-dimensional, non-aging, viscoelastic composite materials. Particularly, it should be pointed out that although we con- sidered some of the results given in Yancey and Pindera (1990) and Wang and Pindera (2016a,b) for comparison with ours, our methodol- ogy provides further developments because of the possibility of consid- ering 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 two- scale AHM to obtain the macroscale functional form of the effective viscoelastic coefficients. In Section 4, we calculate the effective prop- erties 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). Fi- nally, in Section 6, we summarize our findings and highlight the limi- tations of the current model and possible further developments of the work.

## 2. Model description

We identify the heterogeneous, linear viscoelastic material with an open, bounded set Ω ⊂ R 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 L 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,

$$ \epsilon = \frac{\ell}{L} \ll 1, \qquad (1) $$

and the microscopic spatial variable

$$ y = \frac{x}{\epsilon}, \qquad (2) $$

where x is said to be the macroscopic spatial variable.
In particular, we set $\overline{\Omega} = \overline{\Omega}_1^\epsilon \cup \overline{\Omega}_2^\epsilon$ with $\overline{\Omega}_1^\epsilon \cap \Omega_2^\epsilon = \Omega_1^\epsilon \cap \overline{\Omega}_2^\epsilon = \emptyset$, and where $\Omega_2^\epsilon$ denotes the matrix and $\Omega_1^\epsilon = \cup_{i=1}^N \Omega_1^\epsilon$ represent the inclu- sions with $N \in \mathbb{N}$. Additionally, the interface between $\Omega_1^\epsilon$ and $\Omega_2^\epsilon$ (see Fig. 1 (b)) is denoted by $\Gamma^\epsilon$. Referring to Fig. 1 (c), the unitary peri- odic cell $Y$ is considered to be constituted by a fiber (long or short) $Y_1$ in the matrix $Y_2$ so that the periodic cell is given by $\overline{Y} = \overline{Y}_1 \cup \overline{Y}_2$ with $\overline{Y}_1 \cap Y_2 = Y_1 \cap \overline{Y}_2 = \emptyset$. At this scale, the interface between $Y_1$ and $Y_2$ is denoted by $\Gamma_Y$, see also Di Stefano et al. (2020) for an illustration of various periodic cell arrangement in the context of electro-active com- posites.

## 2.1. 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 $\Gamma^\epsilon$, i.e., the matrix and the sub-phases are in ideal contact. Therefore, the balance of linear momentum in $\Omega \setminus \Gamma^\epsilon$ together with the interface conditions read

$$ \nabla \cdot \sigma^{\epsilon}(x, t) = \mathbf{0} \qquad \text{in } (\Omega \setminus \Gamma^{\epsilon}) \times \mathbb{R}, \qquad (3a) $$

$$ [\mathbf{u}^{\epsilon}(x, t)] = \mathbf{0} \qquad \text{on } \Gamma^{\epsilon} \times \mathbb{R}, \qquad (3b) $$

$$ [\sigma^{\epsilon}(x, t)\mathbf{n}^{(y)}] = \mathbf{0} \qquad \text{on } \Gamma^{\epsilon} \times \mathbb{R}, \qquad (3c) $$

where $\sigma^\epsilon$ represents the second-order stress tensor and $\mathbf{u}^\epsilon$ is the dis- placement field. Moreover, $\mathbf{n}^{(y)}$ denotes the outward unit vector to the interface $\Gamma^\epsilon$, and the operator $[\phi^\epsilon]$ describes the jump of $\phi^\epsilon$ across the interface $\Gamma^\epsilon$ between $\Omega_1^\epsilon$ and $\Omega_2^\epsilon$. Notice that the superscript $\epsilon$ is used to indicate the notation $\phi^\epsilon(x, t) = \phi(x, y, t)$ (refer to Di Stefano et al., 2020 for a comprehensive discussion regarding this notation). The sys- tem of Eqs. (3a)-(3c) has to be supplemented with boundary conditions on $\partial\Omega \times \mathbb{R}$ and initial conditions in $\Omega \times \{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 vis- coelastic material so that (Christensen, 1982)

$$ \sigma^{\epsilon}(x, t) = \int_{0}^{t} R^{\epsilon}(x, t - \tau) : \frac{\partial \xi(\mathbf{u}^{\epsilon}(x, \tau))}{\partial \tau} d\tau, \qquad (4) $$

where $R^\epsilon$ 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 $R^\epsilon$ and, thus, its components are assumed to be smooth functions of $x$ in $(\Omega \setminus \Gamma^\epsilon) \times \mathbb{R}$, but discontinuous on $\Gamma^\epsilon \times \mathbb{R}$. Furthermore, we notice that, in Eq. (4), $\xi$ denotes the second-order strain tensor for small displacements, namely

$$ \xi(\mathbf{u}^{\epsilon}(x, t)) = \frac{1}{2}(\nabla\mathbf{u}^{\epsilon}(x, t) + (\nabla\mathbf{u}^{\epsilon}(x, t))^T), \qquad (5) $$

and we require both minor and major symmetry properties for $\mathbb{R}$, i.e.,

$$ R_{ijkl}^{\epsilon} = R_{jkl}^{\epsilon} = R_{ijkl}^{\epsilon} = R_{klj}^{\epsilon}. \qquad (6) $$

The integral in Eq. (4), standing for the stress-strain relationship for non-aging, viscoelastic materials, can be manipulated by means of in- tegral transforms. In particular, the Laplace–Carson transform, which is given by

$$ \hat{\phi}^{\epsilon}(x, p) = p \int_{0}^{\infty} e^{-p\tau} \phi^{\epsilon}(x, t) dt, \quad \forall t \ge 0, \qquad (7) $$

where $p$ 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 orig- inally proposed by Hashin (1965) and known as the elastic-viscoelastic

correspondence principle, continues to gain interest in the scientific lit- erature (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

$$ \nabla \cdot [\hat{\mathbf{R}}^{\epsilon}(x, p) : \xi(\hat{\mathbf{u}}^{\epsilon}(x, p))] = \mathbf{0} \qquad \text{in } (\Omega \setminus \Gamma^{\epsilon}) \times [0, +\infty) \qquad (8a) $$

$$ [\hat{\mathbf{u}}^{\epsilon}(x, p)] = \mathbf{0} \qquad \text{on } \Gamma^{\epsilon} \times [0, +\infty) \qquad (8b) $$

$$ [\hat{[\mathbf{R}}^{\epsilon}(x, p) : \xi(\hat{\mathbf{u}}^{\epsilon}(x, p))]\mathbf{n}^{(y)}] = \mathbf{0} \qquad \text{on } \Gamma^{\epsilon} \times [0, +\infty). \qquad (8c) $$

## 3. 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 follow- ing relation for the spatial derivatives holds

$$ \frac{\partial \hat{\phi}_l^\epsilon(x, p)}{\partial x_j} = \frac{\partial \hat{\phi}_l^\epsilon(x, y, p)}{\partial x_j} + \frac{1}{\epsilon} \frac{\partial \hat{\phi}_l^\epsilon(x, y, p)}{\partial y_j}. \qquad (9) $$

Moreover, Eq. (5) becomes,

$$ \xi_{kl}(\hat{\phi}^{\epsilon}(x, p)) = \xi_{kl}(\hat{\phi}(x, y, p)) + \frac{1}{\epsilon} \xi_{kl}^{(y)}(\hat{\phi}(x, y, p)), \qquad (10) $$

where we have introduced the notation

$$ \xi_{kl}^{(y)}(\hat{\phi}(x, y, p)) = \frac{1}{2} \left( \frac{\partial \hat{\phi}_k(x, y, p)}{\partial y_l} + \frac{\partial \hat{\phi}_l(x, y, p)}{\partial y_k} \right). \qquad (11) $$

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 $\epsilon$. In the Laplace– Carson domain it reads

$$ \hat{\mathbf{u}}^{\epsilon}(x, p) = \sum_{i=0}^{+\infty} \epsilon^{i} \hat{\mathbf{u}}^{(i)}(x, y, p), \qquad (12) $$

where the coefficients $\hat{\mathbf{u}}^{(i)}(x, y, p)$ are assumed to be periodic in the micro- scopic variable $y$. 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 orig- inal problem (8a)-(8c) and by equating the result in the same powers of $\epsilon$, we obtain that, in the limit $\epsilon \to 0$,

$$ \begin{aligned} \hat{u}_m^\epsilon(x, p) &= \hat{u}_m^{(0)}(x, y, p) + \hat{u}_m^{(1)}(x, y, p)\epsilon + o(\epsilon) \\ &= \hat{v}_m(x, p) + \hat{\chi}_{klm}(y, p)\xi_{kl}^{(y)}(\hat{v}(x, p))\epsilon + o(\epsilon), \end{aligned} \qquad (13) $$

where $\hat{v}$ is a smooth vector function of $x$ and $p$, and the third order tensor $\hat{\chi}_{klm}$ is the solution of the $\epsilon$-local problem given by

$$ \frac{\partial}{\partial y_j} [\hat{\mathbf{R}}_{ijsq}(y, p) \xi_{sq}^{(y)}(\hat{\chi}_{kl}(y, p)) + \hat{\mathbf{R}}_{ijkl}(y, p)] = 0 \qquad \text{in } (Y \setminus \Gamma_Y) \times [0, +\infty), \qquad (14a) $$

$$ [\hat{\chi}_{klm}(x, y, p)] = 0 \qquad \text{on } \Gamma_Y \times [0, +\infty), \qquad (14b) $$

$$ [\hat{[\mathbf{R}}_{ijsq}(y, p) \xi_{sq}^{(y)}(\hat{\chi}_{kl}(y, p)) + \hat{\mathbf{R}}_{ijkl}(y, p)] \mathbf{n}_j^{(y)}] = 0 \qquad \text{on } \Gamma_Y \times [0, +\infty), \qquad (14c) $$

$$ \hat{\chi}_{klm}(y, p) = 0 \qquad \text{in } Y \times \{0\}, \qquad (14d) $$

where

$$ \xi_{sq}^{(y)}(\hat{\chi}_{kl}(y, p)) = \frac{1}{2} \left( \frac{\partial \hat{\chi}_{kls}(y, p)}{\partial y_q} + \frac{\partial \hat{\chi}_{klq}(y, p)}{\partial y_s} \right). \qquad (15) $$
The uniqueness of the solution of the local problem (14a)–(14d) is guaranteed by enforcing either, the condition $\langle \hat{\chi}_{klm} \rangle_y = 0$ or by fixing the value of $\hat{\chi}_{klm}$ at one point of the reference periodic cell Y (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 $\epsilon$-local problem has to be supplemented with an initial condition in $Y \times \{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 $\epsilon^0$. Specifically, this can be written as

$$ \hat{R}_{ijkl}^{(*)}(p) \frac{\partial}{\partial x_j} \xi_{kl}(\hat{\phi}(x, p)) = 0 \quad \text{in } \Omega^h \times [0, +\infty), \qquad (16a) $$

where

$$ \hat{R}_{ijkl}^{(*)}(p) := \langle \hat{R}_{ijkl}(y, p) + \hat{R}_{ijmn}(y, p) \rangle_{\xi_{mn}^{(y)}}(\hat{\chi}_{kl}(y, p)) \rangle_y. \qquad (17) $$

is the effective relaxation modulus in the Laplace–Carson space. In Eq. (17), the notation $\langle \phi \rangle_y$ denotes the cell average operator and is defined by the expression

$$ \langle \phi \rangle_y = \frac{1}{|Y|} \int_Y \phi \, dy, \qquad (18) $$

with $|Y|$ being the volume of the periodic cell Y.

## 4. Calculation of the effective properties

For simplicity in our calculations, we consider that the relaxation modulus, $\hat{R}_{ijkl}$, is y-constant in each constituent of the periodic cell Y, i.e.

$$ \hat{R}_{ijkl}(y, p) = \begin{cases} \hat{R}_{ijkl}^{(1)}(p), & \text{if } y \in Y_1, \\ \hat{R}_{ijkl}^{(2)}(p), & \text{if } y \in Y_2, \end{cases} \qquad (19) $$

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,

$$ \frac{\partial}{\partial y_j} \left[ \hat{R}_{ijsq}^{(1)}(p) \xi_{sq}^{(y)} \left( \hat{\chi}_{kl}^{(1)}(y, p) \right) \right] = 0 \quad \text{in } Y_1 \times [0, +\infty), \qquad (20a) $$

$$ \frac{\partial}{\partial y_j} \left[ \hat{R}_{ijsq}^{(2)}(p) \xi_{sq}^{(y)} \left( \hat{\chi}_{kl}^{(2)}(y, p) \right) \right] = 0 \quad \text{in } Y_2 \times [0, +\infty), \qquad (20b) $$

$$ \hat{\chi}_{klm}^{(1)}(y, p) = \hat{\chi}_{klm}^{(2)}(y, p) \quad \text{on } \Gamma_Y \times [0, +\infty), \qquad (20c) $$

$$ \left[ \hat{R}_{ijsq}^{(1)}(p) \xi_{sq}^{(y)} \left( \hat{\chi}_{kl}^{(1)}(y, p) \right) \right] n_j^{(y)} - \left[ \hat{R}_{ijsq}^{(2)}(p) \xi_{sq}^{(y)} \left( \hat{\chi}_{kl}^{(2)}(y, p) \right) \right] n_j^{(y)} = f_{ikl}^{(y)}(p) \quad \text{on } \Gamma_Y \times [0, +\infty), \qquad (20d) $$

$$ \hat{\chi}_{klm}(y, p) = 0 \quad \text{in } Y \times \{0\}. \qquad (20e) $$

We notice that, in Eq. (20d), the stress jump conditions lead to the occurrence of interface loads, i.e.,

$$ f_{ikl}^{(y)}(p) = \left[ \hat{R}_{ijkl}^{(2)}(p) - \hat{R}_{ijkl}^{(1)}(p) \right] n_j^{(y)}, \qquad (21) $$

which arise as a consequence of the discontinuities in the coefficients of $\hat{R}$ between the host medium and the sub-phases, and represent the driving force to obtain nontrivial solutions of the six elastic-type local problems ((k, l), k ≥ l) (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 $f_{ikl}^{(y)}$ read

$$ f_{11}^{(y)}(p) = [\hat{R}_{11}^{(2)}(p) - \hat{R}_{11}^{(1)}(p)]n_1^{(y)}\mathbf{e}_1 + [\hat{R}_{21}^{(2)}(p) - \hat{R}_{21}^{(1)}(p)]n_2^{(y)}\mathbf{e}_2 \\ + [\hat{R}_{31}^{(2)}(p) - \hat{R}_{31}^{(1)}(p)]n_3^{(y)}\mathbf{e}_3, \qquad (22a) $$

$$ f_{22}^{(y)}(p) = [\hat{R}_{12}^{(2)}(p) - \hat{R}_{12}^{(1)}(p)]n_1^{(y)}\mathbf{e}_1 + [\hat{R}_{22}^{(2)}(p) - \hat{R}_{22}^{(1)}(p)]n_2^{(y)}\mathbf{e}_2 \\ + [\hat{R}_{32}^{(2)}(p) - \hat{R}_{32}^{(1)}(p)]n_3^{(y)}\mathbf{e}_3, \qquad (22b) $$

$$ f_{33}^{(y)}(p) = [\hat{R}_{13}^{(2)}(p) - \hat{R}_{13}^{(1)}(p)]n_1^{(y)}\mathbf{e}_1 + [\hat{R}_{23}^{(2)}(p) - \hat{R}_{23}^{(1)}(p)]n_2^{(y)}\mathbf{e}_2 \\ + [\hat{R}_{33}^{(2)}(p) - \hat{R}_{33}^{(1)}(p)]n_3^{(y)}\mathbf{e}_3, \qquad (22c) $$

$$ f_{23}^{(y)}(p) = f_{32}^{(y)}(p) = [\hat{R}_{44}^{(2)}(p) - \hat{R}_{44}^{(1)}(p)]n_3^{(y)}\mathbf{e}_2 + [\hat{R}_{44}^{(2)}(p) - \hat{R}_{44}^{(1)}(p)]n_2^{(y)}\mathbf{e}_3, \qquad (22d) $$

$$ f_{13}^{(y)}(p) = f_{31}^{(y)}(p) = [\hat{R}_{55}^{(2)}(p) - \hat{R}_{55}^{(1)}(p)]n_3^{(y)}\mathbf{e}_1 + [\hat{R}_{55}^{(2)}(p) - \hat{R}_{55}^{(1)}(p)]n_1^{(y)}\mathbf{e}_3, \qquad (22e) $$

$$ f_{12}^{(y)}(p) = f_{21}^{(y)}(p) = [\hat{R}_{66}^{(2)}(p) - \hat{R}_{66}^{(1)}(p)]n_2^{(y)}\mathbf{e}_1 + [\hat{R}_{66}^{(2)}(p) - \hat{R}_{66}^{(1)}(p)]n_1^{(y)}\mathbf{e}_2, \qquad (22f) $$

where $\{e_i\}_{i=1}^3$ represents the standard vector basis.

### 4.1. Numerical approach

At this point, it is possible to solve numerically the set of elastic-type 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™ for Matlab® scripting (see Cruz-González et al., 2020b; Penta and Gerisch, 2016; Penta and Gerisch, 2017). Particularly, once $\hat{R}^{(*)}$ is known, the effective creep compliance $\hat{J}^{(*)}$ in the Laplace–Carson space can be computed through the relationship

$$ \hat{R}_{ijmn}^{(*)}(p) \hat{J}_{mnkl}^{(*)}(p) = I_{ijkl}, \qquad (23) $$

where $I_{ijkl}$ 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,

(a) Discretize the time interval $t = [t_1, t_2, ..., t_N]$

(b) For each $t_i$, obtain the components $p_j := \alpha_j/t_i$ and $B_j := \beta_j/t_i$ for $j = 1, ..., (ns + nd + 1)$, where ns and nd are implicit parameters, and $\alpha$ and $\beta$ are defined in Valsa's method.

(c) Calculate $\hat{R}^{(*)}(p_j)$ and $\hat{J}^{(*)}(p_j)$ for $j = 1, ..., (ns + nd + 1)$.

(d) Use the last step of Valsa's method to determine the effective coefficients in the time domain

$$ \mathcal{R}^{(*)}(t_i) = \sum_{j=1}^{ns+nd+1} \text{Re}[B_j \hat{R}^{(*)}(p_j)/p_j] \quad \text{for } i = 1, ..., N, \qquad (24a) $$

$$ \mathcal{J}^{(*)}(t_i) = \sum_{j=1}^{ns+nd+1} \text{Re}[B_j \hat{J}^{(*)}(p_j)/p_j] \quad \text{for } i = 1, ..., N, \qquad (24b) $$

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 $\mathcal{R}^{(*)}$ and $\mathcal{J}^{(*)}$ by an angle $\theta$ and obtain $\mathcal{R}_{\theta}^{(*)}$ and $\mathcal{J}_{\theta}^{(*)}$. In this respect, the following transformations are useful,

$$ \mathcal{R}_{\theta}^{(*)} = \mathcal{K} \mathcal{R}^{(*)} \mathcal{K}^{T}, \qquad (25a) $$

$$ \mathcal{J}_{\theta}^{(*)} = (\mathcal{K}^{-1})^{T} \mathcal{J}^{(*)} \mathcal{K}^{-1}, \qquad (25b) $$
<figure><img src="image_5.png" /><figcaption>Fig. 2. Methodology sketch of AHMFE.</figcaption></figure>

where

$$[\mathcal{K}] = \begin{bmatrix} Q_{11}^2 & Q_{12}^2 & Q_{13}^2 & 2Q_{12}Q_{13} & 2Q_{13}Q_{11} & 2Q_{11}Q_{12} \\ Q_{21}^2 & Q_{22}^2 & Q_{23}^2 & 2Q_{22}Q_{23} & 2Q_{23}Q_{21} & 2Q_{21}Q_{22} \\ Q_{31}^2 & Q_{32}^2 & Q_{33}^2 & 2Q_{32}Q_{33} & 2Q_{33}Q_{31} & 2Q_{31}Q_{32} \\ Q_{21}Q_{31} & Q_{22}Q_{32} & Q_{23}Q_{33} & Q_{22}Q_{33} + Q_{23}Q_{32} & Q_{23}Q_{31} + Q_{21}Q_{33} & Q_{21}Q_{32} + Q_{22}Q_{31} \\ Q_{31}Q_{11} & Q_{32}Q_{12} & Q_{33}Q_{13} & Q_{32}Q_{13} + Q_{33}Q_{12} & Q_{33}Q_{11} + Q_{31}Q_{13} & Q_{31}Q_{12} + Q_{32}Q_{11} \\ Q_{11}Q_{21} & Q_{12}Q_{22} & Q_{13}Q_{23} & Q_{12}Q_{23} + Q_{13}Q_{22} & Q_{13}Q_{21} + Q_{11}Q_{23} & Q_{11}Q_{22} + Q_{12}Q_{21} \end{bmatrix}$$

and $Q_{ij}$ ($i, j = 1, 2, 3$) are the coefficients of the orthogonal rotation tensor (see Ramírez-Torres et al., 2018; 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.

# 5. Results

## 5.1. 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-
<figure><img src="image_6.png" /><figcaption>Fig. 3. Mesh discretizations A, B and C for the hexagonal periodic cell with $V_f=0.6$.</figcaption></figure>

<table><caption>**Table 1**

Elastic properties of the constituents.</caption><thead><tr><th>Materials</th><th>E<sub>A</sub> (GPa)</th><th>E<sub>T</sub> (GPa)</th><th>μ<sub>A</sub> (GPa)</th><th>μ<sub>T</sub> (GPa)</th><th>ν<sub>A</sub></th></tr></thead><tbody><tr><td>AS4 graphite fiber</td><td>225</td><td>15</td><td>15</td><td>7</td><td>0.20</td></tr><tr><td>3501-6 epoxy</td><td>4.2</td><td>4.2</td><td>1.567</td><td>1.567</td><td>0.34</td></tr><tr><td>E-glass fiber</td><td>69.0</td><td>69.0</td><td>28.28</td><td>28.28</td><td>0.22</td></tr><tr><td>Boron fiber</td><td>420</td><td>420</td><td>175</td><td>175</td><td>0.20</td></tr><tr><td>Aluminum</td><td>69.0</td><td>69.0</td><td>25.94</td><td>25.94</td><td>0.33</td></tr></tbody></table>

<table><caption>**Table 2**

Maximum relative error (%).</caption><thead><tr><th>AHMFE → LEHT</th><th>E<sup>(*)</sup><sub>T</sub>/E<sub>m</sub></th><th>ν<sup>(*)</sup><sub>T</sub></th><th>μ<sup>(*)</sup><sub>A</sub>/μ<sub>m</sub></th></tr></thead><tbody><tr><td>glass/epoxy</td><td>0.5844</td><td>0.3969</td><td>0.0767</td></tr><tr><td>graphite/epoxy</td><td>0.1971</td><td>0.1311</td><td>0.0604</td></tr><tr><td>boron/aluminum</td><td>0.3476</td><td>0.2088</td><td>0.0604</td></tr><tr><td>aluminum/porosity</td><td>1.4575</td><td>1.1796</td><td>0.0677</td></tr></tbody></table>

niques. In this context, the elastic limit case is reached by considering $t=0$ in Eq. (4), which implies that $R_{ijmn}^{(*)}$ 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 $E_A$ ($E_T$) and $\mu_A$ ($\mu_T$) is used for the axial (transverse) Young's and shear moduli, respectively, and $\nu_A$ 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 $E_T^{(*)}$ since those related with the transverse Poisson's ratio $\nu_T^{(*)}$ and the axial shear modulus $\mu_A^{(*)}$ 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 $E_T^{(*)}$, normalized by the corresponding matrix Young's modulus $E_m$. 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 $V_f$. 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

$$ \text{Relative Error}(\phi) := \frac{\phi^{(\text{FVDAM})} - \phi^{(\text{AHMFE})}}{\phi^{(\text{FVDAM})}} \times 100\%, \quad (26) $$

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 FVDAM, 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 $V_f = 0.05$ to $V_f = 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 $E_T^{(*)}/E_m$, $\nu_T^{(*)}$ and $\mu_A^{(*)}/\mu_m$ 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%.
<figure><img src="image_7.png" /><figcaption>Fig. 4. Chart (a) shows the computation of the normalized effective Young's modulus of a graphite/epoxy system with hexagonal architecture. Chart (b) displays the relative error of AHMFE in relation to FVDAM. Chart (c) shows the total computing time and the three charts (d), (e) and (f) provide the specific computing time for each mesh, respectively.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Calculation of the effective moduli for unidirectional composites with different contrast in the constituents and hexagonal periodic cell. The comparisons are performed with Fig. 9 of Wang and Pindera (2016b).</figcaption></figure>

## 5.2. 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
<table><caption>**Table 3**
Maximum relative error (%) in the time interval under study.</caption><thead><tr><th rowspan="2">AHMFE → LEHT</th><th colspan="3">Constant Poissons ratio</th><th colspan="3">Constant bulk modulus</th></tr><tr><th>J<sub>22</sub><sup>(*)</sup></th><th>J<sub>44</sub><sup>(*)</sup></th><th>J<sub>66</sub><sup>(*)</sup></th><th>J<sub>22</sub><sup>(*)</sup></th><th>J<sub>44</sub><sup>(*)</sup></th><th>J<sub>66</sub><sup>(*)</sup></th></tr></thead><tbody><tr><td>Hexagonal cell</td><td>0.2840</td><td>0.2605</td><td>0.2622</td><td>0.2707</td><td>0.2345</td><td>0.2901</td></tr><tr><td>Square cell</td><td>1.7953</td><td>0.0994</td><td>0.3588</td><td>2.2884</td><td>0.2034</td><td>0.3230</td></tr></tbody></table>

the LEHT. We mention that in the upcoming simulations, we employ Mesh B.

### 5.2.1. 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 ($v_m$) or the bulk modulus ($K_m$) 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 $v_m = 0.38$ as the value for the matrix and therefore the time-dependent bulk modulus is given by $K_m(t) = E_m(t)/(3(1 - 2v_m))$, where $E_m(t)$ 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 $K_m = E_0/(3(1 - 2v_m))$ to determine its value, and then, the time-dependent Poisson's ratio $v_m(t)$ arises from the equation $v_m(t) = 1/2 - E_m(t)/(6K_m)$. Here, $E_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 $E_f = 68.77$ GPa and Poisson's ratio $v_f = 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,

$$E_m(t) = G_1 \exp\left(-\frac{t}{\eta_{\sigma,1}}\right) + G_2 \exp\left(-\frac{t}{\eta_{\sigma,2}}\right), \quad (27)$$

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, $G_1 = 1.12511$ GPa, $G_2 = 2.14489$ GPa, $\eta_{\sigma,1} = 6999.34$ h and $\eta_{\sigma,2} = 58.2551$ h, where $G_n$ ($n = 1, 2$) represents the elastic modulus of the spring and $\eta_{\sigma,n}$ is a relaxation time Mainardi and Spada (2011). It is worth to remark that from eq. (27), we obtain $E_0 = G_1 + G_2$.

In Fig. 6, we show the curves corresponding to the effective creep compliances $J_{22}^{(*)}$, $J_{44}^{(*)}$ and $J_{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.

### 5.2.2. 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><caption>**Table 4**
Elastic properties of the transversely isotropic T300 graphite fiber at room (22°C) and elevated (121°C) temperature.</caption><thead><tr><th>Temperature</th><th>E<sub>A</sub> (GPa)</th><th>E<sub>T</sub> (GPa)</th><th>μ<sub>A</sub> (GPa)</th><th>ν<sub>A</sub></th><th>ν<sub>T</sub></th></tr></thead><tbody><tr><td>22°C</td><td>202.82</td><td>25.30</td><td>44.12</td><td>0.443</td><td>0.05</td></tr><tr><td>121°C</td><td>214.33</td><td>14.82</td><td>68.18</td><td>0.450</td><td>0.05</td></tr></tbody></table>

Wang and Pindera (2016a). In Yancey and Pindera (1990), the authors observed that at 22°C and 121°C the T300 graphite fiber present an elastic behavior, whereas the creep response of the 934 epoxy matrix is fitted by the power-law

$$S_m(t) = \frac{1}{E_0} + Ct^n, \quad (28)$$

where C and n are experimentally measured parameters and $E_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 $\hat{E}_m(p)$ which is needed in our simulations.

With this purpose, we introduce the notation $\mu = E_0/2$, $\tau = \beta^{-1/n}$, $\beta = E_0 C \Gamma(1 + n)$ and $\alpha = n$, so that Eq. (28) can be equivalently rewritten as

$$S_m(t) = S_M(t) = \frac{1}{2\mu} \left[ 1 + \frac{(t/\tau)^\alpha}{\Gamma(1+\alpha)} \right], \quad (29)$$

where $\Gamma$ denotes the Gamma function, $\alpha \in [0, 1]$, and $S_M(t)$ 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

$$E_M(t) = 2\mu\mathcal{E}_\alpha(-\langle t/\tau \rangle^\alpha), \quad (30)$$

where $\mathcal{E}_\alpha$ denotes the Mittag-Leffler function of order $\alpha$, which is defined by the expression (Gorenflo et al., 2014)

$$\mathcal{E}_{\alpha}(-\langle t/\tau \rangle^{\alpha}) = \sum_{p=0}^{\infty} (-1)^p \frac{(t/\tau)^{\alpha p}}{\Gamma(1 + \alpha p)}, \quad 0 < \alpha < 1, \quad \tau > 0, \quad (31)$$

and for $\alpha = 1$ reduces to $\mathcal{E}_\alpha((-t/\tau)^\alpha) = \exp(-t/\tau)$. We notice that, from Eq. (30) and by virtue of the parameter identifications established for obtaining Eq. (29), obtain

$$E_m(t) = E_0 \mathcal{E}_n(-E_0 C \Gamma(1 + n)t^n). \quad (32)$$

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 $\hat{E}_m(p)$, 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

$$\hat{E}_M(p) = \frac{2\mu(p\tau)^\alpha}{1 + (p\tau)^\alpha}, \quad (33)$$

which implies that

$$\hat{E}_m(p) = \frac{E_0 p^n}{p^n + E_0 C \Gamma(1 + n)}. \quad (34)$$

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.
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Charts (a)-(c) show the results under the considerations of constant Poisson's ratio while charts (d)-(f) assume constant bulk modulus. We consider square and hexagonal arrays of fibers. Our results are compared with Fig. 7 of Wang and Pindera (2016a).</figcaption></figure>

<table><caption>**Table 5**

Material properties of the epoxy matrix at room (22°C) and elevated (121°C) temperature.</caption><thead><tr><td>Temperature</td><td>E<sub>0</sub> (GPa)</td><td>ν</td><td>C (1/(GPa×min))</td><td>n</td></tr></thead><tbody><tr><td>22°C</td><td>4.51</td><td>0.311</td><td>0.0135</td><td>0.17</td></tr><tr><td>121°C</td><td>3.36</td><td>0.317</td><td>0.0250</td><td>0.20</td></tr></tbody></table>

Fig. 7 shows the effective creep response of the viscoelastic composite material given by the coefficients $J_{11}^{(*)}$ and $J_{22}^{(*)}$. 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°C (room) and 121°C (elevated). The results correspond to composites with fibers that are rotated 10° and 90° counterclockwise about the y₂-axis. In this case, the fiber volume fraction of the fiber is fixed to $V_f = 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.

### 5.2.3. Modeling of short fiber reinforcement

The locally-exact homogenization theory (LEHT) is based on a two-dimensional formulation which is only capable of taking into account
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Calculation of the effective creep response in viscoelastic composites with hexagonal array of transversely isotropic fiber with 10° and 90° off-axis specimens about the x₂-axis. We consider two different temperatures, i.e. 22°C (room) and 121°C (elevated). In addition, Charts (a)-(b) shows the cases of constant Poissons ratio and charts (c)-(d) the constant bulk modulus.</figcaption></figure>

<table><caption>**Table 6**
Maximum relative error (%) in the time interval under study.</caption><thead><tr><th rowspan="2">Hexagonal cell</th><th colspan="2">Constant Poissons ratio</th><th colspan="2">Constant bulk modulus</th></tr><tr><th>θ = 10°</th><th>θ = 90°</th><th>θ = 10°</th><th>θ = 90°</th></tr></thead><tbody><tr><td>AHMFE → LEHT</td><td>J<sub>11</sub><sup>(*)</sup></td><td>J<sub>22</sub><sup>(*)</sup></td><td>J<sub>11</sub><sup>(*)</sup></td><td>J<sub>22</sub><sup>(*)</sup></td></tr><tr><td>22°C</td><td>0.0798</td><td>0.3287</td><td>0.0870</td><td>0.3799</td></tr><tr><td>121°C</td><td>0.1340</td><td>0.3847</td><td>0.1781</td><td>0.1666</td></tr><tr><td>AHMFE → Experiments</td><td>J<sub>11</sub><sup>(*)</sup></td><td>J<sub>22</sub><sup>(*)</sup></td><td>J<sub>11</sub><sup>(*)</sup></td><td>J<sub>22</sub><sup>(*)</sup></td></tr><tr><td>22°C</td><td>4.2505</td><td>3.3196</td><td>4.5503</td><td>5.2098</td></tr><tr><td>121°C</td><td>5.6338</td><td>4.7281</td><td>6.7843</td><td>1.7558</td></tr></tbody></table>

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 := h_1/H_1$ and $γ_2 := h_2/H_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 y₂-axis of $0 ≤ θ < π$.

In Fig. 8 (b) and (c), we show our findings in the calculation of the effective moduli $E_1^{(*)}$ and $ν_{32}^{(*)}$ for square and hexagonal periodic cell. In particular, we assume the room temperature (22°C) and the constant bulk modulus approach discussed in the Section 5.2.2. The fiber volume fraction and the ratio $γ$ are fixed to $V_f = 0.1$ and $γ = γ_1 = γ_2 =$
<figure><img src="image_11.png" /><figcaption>Fig. 8. Chart (a) shows the $\theta$ counterclockwise rotation about $y_2$-axis, and the square and hexagonal periodic cell for short fibers. Charts (b) and (c) show the results corresponding to $E_1^{(*)}$ and $v_{32}^{(*)}$.</figcaption></figure>

$\theta \in \{0, \pi/2\}$.

## 6. Discussion and conclusions

In this work, we studied the effective properties of three-dimensional, 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 $\gamma$ 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 $\theta$ in Fig. 8(b) and (c), when the rotation increases, we obtain that the values of the effective moduli $E_1^{(*)}$ and $v_{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 $ns + nd + 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 con-
sideration 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 condi-
tions between the interface of the constituents (see Daridon et al.,
2016; Escarpini Filho and Marques, 2016; Rizzoni et al., 2021; Ser-
pilli et al., 2019). Besides, another development is related to the mod-
eling 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 reinforce-
ment in hierarchical composite materials. Another challenging exten-
sion 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 Ramírez-Torres et al. (2018) and typically found in relevant bio-
logical 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 nat-
ural 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.

Acknowledgements

We gratefully acknowledge to Dr. Zhelong He and Prof. Marek-Jerzy Pindera for the valuable assistance in providing us with the data for the comparisons. OLCG kindly thanks to Ecole Doctorale no. 353 de L'Université de Aix Marseille and L'équipe Matériaux & Structures du Laboratoire de Mécanique et d'Acoustique LMA - UMR 7031 AMU - CNRS - Centrale Marseille 4 impasse Nikola Tesla CS 40006 13453 Marseille Cedex 13, France. RRR thanks to Departamento de Matemáticas y Mecánica, IIMAS and PREI-DGAPA at UNAM, for its support to his research project and the Aix-Marseille Université. RP is partially funded by EPSRC grants EP/S030875/1 and EP/T017899/1.

References

Allaire, G., Briane, M., 1996. Multiscale convergence and reiterated homogenisation. Proc. R. Soc. Edinburgh Sec. A Math. 126 (2), 297-342. doi:10.1017/S0308210500022757.

Amiri-Rad, A., Pastukhov, L., Govaert, L., van Dommelen, J., 2019. An anisotropic viscoelastic-viscoplastic model for short-fiber composites. Mech. Mater. 137, 103141. doi:10.1016/j.mechmat.2019.103141.

Atanacković, T.M., Konjik, S., Pilipović, S., Zorica, D., 2016. Complex order frac-
tional derivatives in viscoelasticity. Mech. Time-Depend. Mater. 20 (2), 175-195.
doi:10.1007/s11043-016-9290-3.

Atanacković, T.M., Pilipović, S., Stanković, B., Zorica, D., 2014. Fractional Calculus with Applications in Mechanics: Vibrations and Diffusion Processes. John Wiley & Sons.

Atthapreyangkul, A., Hoffman, M., Pearce, G., 2021. Effect of geometrical structure variations on the viscoelastic and anisotropic behaviour of cortical bone using multi-scale finite element modelling. J. Mech. Behav. Biomed. Mater. 113, 104153. doi:10.1016/j.jmbbm.2020.104153.

Auriault, J.L., Boutin, C., Geindreau, C., 2009. Homogenization of Coupled Phenomena in Heterogenous Media. ISTE doi:10.1002/9780470612033.

Bakhvalov, N.S., Panasenko, G., 1989. Homogenisation: Averaging Processes in Periodic Media: Mathematical Problems in the Mechanics of Composite Materials. In: Mathematics and its Applications. Kluwer Academic Publishers doi:10.1007/978-94-009-2247-1.

Beltempo, A., Bonelli, A., Bursi, O.S., Zingales, M., 2019. A numerical integration approach for fractional-order viscoelastic analysis of hereditary-aging structures. Int. J. Numer. Methods Eng. 121 (6), 1120-1146. doi:10.1002/nme.6259.

Beltempo, A., Bursi, O.S., Cappello, C., Zonta, D., Zingales, M., 2017. A viscoelastic model for the long-term deflection of segmental prestressed box girders. Comput.-Aided Civil Infrastruct. Eng. 33 (1), 64-78. doi:10.1111/mice.12311.

Bensoussan, A., Papanicolau, G., Lions, J.-L., 1978. Asymptotic Analysis for Periodic Structures, 5, 1st ed. North-Holland.

Bouras, Y., Zoric, D., Atanacković, T.M., Vrcelj, Z., 2018. A non-linear thermo-viscoelastic rheological model based on fractional derivatives for high temperature creep in concrete. Appl. Math. Model. 55, 551-568.

Cavalcante, M.A.A., Marques, S.P.C., 2014. Homogenization of periodic materials with viscoelastic phases using the generalized FVDAM theory. Comput. Mater. Sci. 87, 43-53. doi:10.1016/j.commatsci.2014.01.053.

Cepero-Méjías, F., Curiel-Sosa, J., Blázquez, A., Yu, T., Kerrigan, K., Phadnis, V., 2020. Review of recent developments and induced damage assessment in the modelling of the machining of long fibre reinforced polymer composites. Compos. Struct. 240, 112006. doi:10.1016/j.compstruct.2020.112006.

Chen, Q., Chatzigeorgiou, G., Meraghni, F., 2021. Extended mean-field homogenization of viscoelastic-viscoplastic polymer composites undergoing hybrid progressive degradation induced by interface debonding and matrix ductile damage. Int. J. Solids Struct. 210-211, 1-17. doi:10.1016/j.ijsolstr.2020.11.017.

Chen, Q., Wang, G., Chen, X., Geng, J., 2017. Finite-volume homogenization of elastic/viscoelastic periodic materials. Compos. Struct. 182, 457-470. doi:10.1016/j.compstruct.2017.09.044.

Christensen, R.M., 1982. Theory of Viscoelasticity - 2nd Edition An Introduction. Academic Press, Cambridge, MA.

Cioranescu, D., Donato, P., 1999. An Introduction to Homogenization. Oxford University Press.

Cruz-González, O.L., Ramírez-Torres, A., Rodríguez-Ramos, R., Penta, R., Bravo-Castillero, J., Guinovart-Díaz, R., Merodio, J., Sabina, F.J., Lebon, F., 2020a. A hierarchical asymptotic homogenization approach for viscoelastic composites. Mech. Adv. Mater. Struct. doi:10.1080/15376494.2020.1722872.

Cruz-González, O.L., Rodríguez-Ramos, R., Otero, J.A., Ramírez-Torres, A., Penta, R., Lebon, F., 2020b. On the effective behavior of viscoelastic composites in three dimensions. Int. J. Eng. Sci. 157, 103377. doi:10.1016/j.ijengsci.2020.103377.

Daridon, L., Licht, C., Orankitjaroen, S., Pagano, S., 2016. Periodic homogenization for Kelvin-Voigt viscoelastic media with a Kelvin-Voigt viscoelastic interphase. Eur. J. Mech. A/Solids 58, 163-171. doi:10.1016/j.euromechol.2015.12.007.

Destrade, M., Saccomandi, G., 2007. Creep, recovery, and waves in a nonlinear fiber-reinforced viscoelastic solid. SIAM J. Appl. Math. 68 (1), 80-97. doi:10.1137/060664483.

Di Stefano, S., Miller, L., Grillo, A., Penta, R., 2020. Effective balance equations for electrostrictive composites. Zeitschrift für angewandte Mathematik und Physik 71 (5), 1-36. doi:10.1007/s00033-020-01365-x.

Escarpini Filho, R. d. S., Marques, S.P.C., 2016. A model for homogenization of linear viscoelastic periodic composite materials with imperfect interface. Latin Am. J. Solids Struct. 13, 2706-2735. doi:10.1590/1679-78252749.

Gorenflo, R., Kilbas, A.A., Mainardi, F., Rogosin, S.V., 2014. Mittag-Leffler Functions, Related Topics and Applications. Springer Publishing Company, Incorporated.

Hashin, Z., 1965. Viscoelastic behavior of heterogeneous media. J. Appl. Mech. Trans. ASME. 32 (3), 630-636. doi:10.1115/1.3627270.

Hashin, Z., 1972. Theory of fiber reinforced materials. NASA contractor report. NASA CR-1974.

Huang, H.-B., Huang, Z.-M., 2020. Micromechanical prediction of elastic-plastic behavior of a short fiber or particle reinforced composite. Compos. Part A: Appl. Sci. Manuf. 134, 105889. doi:10.1016/j.compositesa.2020.105889.

Juraj, 2020. Numerical inversion of laplace transforms in matlab https://www.mathworks.com/matlabcentral/fileexchange/32824-numerical-inversion-of-laplace-transforms-in-matlab

Kern, L.S., Hine, P.J., Gusev, A.A., 2019. Optimizing the damping properties of unidirectional composites by incorporating carbon fibers with a thin viscoelastic coating. Compos. Struct. 208, 879-890. doi:10.1016/j.compstruct.2018.10.043.

Lakes, R., 2009. Viscoelastic Materials. Cambridge University Press doi:10.1017/CBO9780511626722.

Liu, S., Cheng, Y., Liang, B., Cheng, H., Luo, B., Zhang, K., 2020. Sequential homogenization in Laplace domain for viscoelastic properties of composites with functionally graded interphase. Compos. Struct. 113266. doi:10.1016/j.compstruct.2020.113266.

Mainardi, F., 2010. Fractional Calculus and Waves in Linear Viscoelasticity. Imperial College Press doi:10.1142/p614.

Mainardi, F., Spada, G., 2011. Creep, relaxation and viscosity properties for basic fractional models in rheology. Eur. Phys. J. Spec. Top. 193, 133-160. doi:10.1140/epist/e2011-01387-1.

Mishurova, T., Rachmatulin, N., Fontana, P., Oesch, T., Bruno, G., Radi, E., Sevostianov, I., 2018. Evaluation of the probability density of inhomogeneous fiber orientations by computed tomography and its application to the calculation of the effective properties of a fiber-reinforced composite. Int. J. Eng. Sci. 122, 14-29. doi:10.1016/j.ijengsci.2017.10.002.

Muliana, A., Rajagopal, K., 2012. Modeling the response of nonlinear viscoelastic biodegradable polymeric stents. Int. J. Solids Struct. 49 (7-8), 989-1000. doi:10.1016/j.ijsolstr.2011.12.007.

Nonato Da Silva, C., Ciambella, J., Barros, J., Costa, I., 2020. A model for optimizing hooked end steel fibre reinforcements in cracked cement composites. Appl. Eng. Sci. 3, 100011. doi:10.1016/j.apples.2020.100011.

Ojanen, X., Tanska, P., Malo, M., Isaksson, H., Väätänen, S., Koistinen, A., Grassi, L., Magnusson, S., Ribel-Madsen, S., Korhonen, R., Jurvelin, J., Töyräs, J., 2017. Tissue viscoelasticity is related to tissue composition but may not fully predict the apparent-level viscoelasticity in human trabecular bone an experimental and finite element study. J. Biomech. 65, 96-105. doi:10.1016/j.jbiomech.2017.10.002.

Ornaghi Jr., H.L., Neves, R.M., Monticeli, F.M., Almeida Jr., J.H.S., 2020. Viscoelastic characteristics of carbon fiber-reinforced epoxy filament wound laminates. Compos. Commun. 21, 100418. doi:10.1016/j.coco.2020.100418.
Otero, J.A., Rodríguez-Ramos, R., Guinovart-Díaz, R., Cruz-González, O.L., Sabina, F.J., Berger, H., Böhlke, T., 2020. Asymptotic and numerical homogenization methods applied to fibrous viscoelastic composites using Pronys series. Acta Mechanica doi:10.1007/s00707-020-02671-1.

Park, S.W., Kim, Y.R., 1999. Interconversion between relaxation modulus and creep compliance for viscoelastic solids. J. Mater. Civil Eng. 11, 76–82.

Parnell, W.J., De Pascalis, R., 2019. Soft metamaterials with dynamic viscoelastic functionality tuned by pre-deformation. Philos. Trans. R. Soc. A 377 (2144), 20180072. doi:10.1098/rsta.2018.0072.

Penta, R., Gerisch, A., 2016. Investigation of the potential of asymptotic homogenization for elastic composites via a three-dimensional computational study. Comput. Vis. Sci. 17. doi:10.1007/s00791-015-0257-8.

Penta, R., Gerisch, A., 2017. The asymptotic homogenization elasticity tensor properties for composites with material discontinuities. Continuum Mech. Thermodyn. 29 (1), 187–206. doi:10.1007/s00161-016-0526-x.

Pereira, A.S., Thompson, R.L., Mompean, G., 2020. Persistence of straining and polymer alignment in viscoelastic turbulence. Appl. Eng. Sci. 4, 100026. doi:10.1016/j.apples.2020.100026.

Ramírez-Torres, A., Di Stefano, S., Grillo, A., Rodríguez-Ramos, R., Merodio, J., Penta, R., 2018. An asymptotic homogenization approach to the microstructural evolution of heterogeneous media. Int. J. Non-Linear Mech. 106, 245–257. doi:10.1016/j.ijnonlinmec.2018.06.012.

Ramírez-Torres, A., Penta, R., Grillo, A., 2021a. Two-scale, non-local diffusion in homogenised heterogeneous media. Arch. Appl. Mech. doi:10.1007/s00419-020-01880-3.

Ramírez-Torres, A., Di Stefano, S., Grillo, A., 2021b. Influence of non-local diffusion in avascular tumour growth. Math. Mech. Solids doi:10.1177/1081286520975086.

Ramírez-Torres, A., Penta, R., Rodríguez-Ramos, R., Merodio, J., Sabina, F.J., Bravo-Castillero, J., Guinovart-Díaz, R., Preziosi, L., Grillo, A., 2018. Three scales asymptotic homogenization and its application to layered hierarchical hard tissues. Int. J. Solids and Struct. 130-131, 190–198. doi:10.1016/j.ijsolstr.2017.09.035.

Ramírez-Torres, A., Penta, R., Rodríguez-Ramos, R., Grillo, A., 2019. Effective properties of hierarchical fiber-reinforced composites via a three-scale asymptotic homogenization approach. Math. Mech. Solids 24, 3554–3574. doi:10.1177/1081286519847687.

Rizzoni, R., Dumont, S., Lebon, F., Sacco, E., 2021. Higher order adhesive effects in composite beams. Eur. J. Mech. - A/Solids 85, 104108. doi:10.1016/j.euromechsol.2020.104108.

Rodríguez-Ramos, R., Otero, J.A., Cruz-González, O.L., Guinovart-Díaz, R., Bravo-Castillero, J., Sabina, F.J., Padilla, P., Lebon, F., Sevostianov, I., 2020. Computation of the relaxation effective moduli for fibrous viscoelastic composites using the asymptotic homogenization method. Int. J. Solids Struct. 190, 281–290. doi:10.1016/j.ijsolstr.2019.11.014.

Sanchez-Palencia, E., 1980. Non-Homogeneous Media and Vibration Theory. Springer-Verlag doi:10.1007/3-540-10000-8.

Serpilli, M., Rizzoni, R., Lebon, F., Dumont, S., 2019. An asymptotic derivation of a general imperfect interface law for linear multiphysics composites. Int. J. Solids Struct. 180-181, 97–107. doi:10.1016/j.ijsolstr.2019.07.014.

Sevostianov, I., Levin, V., Radi, E., 2016. Effective viscoelastic properties of short-fiber reinforced composites. Int. J. Eng. Sci. 100, 61–73. doi:10.1016/j.ijengsci.2015.10.008.

Sherman, V.R., Tang, Y., Zhao, S., Yang, W., Meyers, M.A., 2017. Structural characterization and viscoelastic constitutive modeling of skin. Acta Biomaterialia 53, 460–469. doi:10.1016/j.actbio.2017.02.011.

Tang, T., Felicelli, S.D., 2015. Computational evaluation of effective stress relaxation behavior of polymer composites. Int. J. Eng. Sci. 90, 76–85. doi:10.1016/j.ijengsci.2015.02.003.

Ting, T.C.T., 1996. Anisotropic Elasticity: Theory and Applications. Oxford University Press.

Tran, A., Yvonnet, J., He, Q.-C., Toulemonde, C., Sanahuja, J., 2011. A simple computational homogenization method for structures made of linear heterogeneous viscoelastic materials. Comput. Methods Appl. Mech. Eng. 200, 2956–2970. doi:10.1016/j.cma.2011.06.012.

Valsa, J., Brančik, L., 1990. Micromechanical analysis of the creep response of unidirectional composites. J. Eng. Mater. Technol. 112 (2), 157–163. doi:10.1115/1.2903302.

Vilchevskaya, E., Levin, V., Seyedkavoosi, S., Sevostianov, I., 2019. Replacement relations for a viscoelastic material containing multiple inhomogeneities. Int. J. Eng. Sci. 136, 26–37. doi:10.1016/j.ijengsci.2018.12.006.

Wang, B., Fang, G., Tan, X., Liang, J., Ge, J., Liu, S., 2020. Investigation on the longitudinal compressive strength of unidirectional carbon fiber/nanoparticles reinforced polymer composites using fft-based method. Compos. Struct. 247, 112448. doi:10.1016/j.compstruct.2020.112448.

Wang, G., Pindera, M.-J., 2016a. Locally-exact homogenization of viscoelastic unidirectional composites. Mech. Mater. 103, 95–109. doi:10.1016/j.mechmat.2016.09.009.

Wang, G., Pindera, M.-J., 2016b. Locally-exact homogenization theory for transversely isotropic unidirectional composites. Mech. Res. Commun. 78, 2–14. doi:10.1016/j.mechrescom.2015.09.011.

Wang, Z., Smith, D.E., 2019. Numerical analysis on viscoelastic creep responses of aligned short fiber reinforced composites. Compos. Struct. 229, 111394. doi:10.1016/j.compstruct.2019.111394.

Yancey, R.N., Pindera, M.-J., 1990. Micromechanical analysis of the creep response of unidirectional composites. J. Eng. Mater. Technol. 112 (2), 157–163. doi:10.1115/1.2903302.

Yang, Z., Sun, Y., Cui, J., Ge, J., 2019. A three-scale asymptotic expansion for predicting viscoelastic properties of composites with multiple configuration. Eur. J. Mech. - A/Solids 76, 235–246. doi:10.1016/j.euromechsol.2019.04.016.

Yi, Y.-M., Park, S.-H., Youn, S.-K., 1998. Asymptotic homogenization of viscoelastic composites with periodic microstructures. Int. J. Solids Struct. 35 (17), 2039–2055. doi:10.1016/S0020-7683(97)00166-2.

Yu, H., Potter, K., Wisnom, M., 2014. A novel manufacturing method for aligned discontinuous fibre composites (high performance-discontinuous fibre method). Compos. Part A: Appl. Sci. Manuf. 65, 175–185. doi:10.1016/j.compositesa.2014.06.005.

Zakian, V., 1969. Numerical inversion of laplace transform. Electron. Lett. 5, 120–121.