# On the effective behavior of viscoelastic composites in three dimensions

O.L. Cruz-González, R. Rodríguez-Ramos, J.A. Otero, A. Ramírez-Torres, R. Penta, Frédéric Lebon

► To cite this version:

O.L. Cruz-González, R. Rodríguez-Ramos, J.A. Otero, A. Ramírez-Torres, R. Penta, et al. On the effective behavior of viscoelastic composites in three dimensions. International Journal of Engineering Science, 2020, 157, pp.103377. 10.1016/j.ijengsci.2020.103377. hal-03085165

HAL Id: hal-03085165

https://hal.science/hal-03085165

Submitted on 26 Oct 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
# On the effective behavior of viscoelastic composites in three dimensions

O.L. Cruz-Gonzálezª, R. Rodríguez-Ramosᵇ,ᶜ,*, J.A. Oteroᵈ, A. Ramírez-Torresᵉ, R. Pentaf, F. Lebonª

ªAix Marseille Univ, CNRS, Centrale Marseille, LMA UMR Marseille, 7031, France

ᵇFacultad de Matemática y Computación, Universidad de La Habana, San Lázaro y L, Vedado, La Habana CP 10400 Cuba

ᶜEscuela de Ingeniería y Ciencias, Tecnológico de Monterrey, Campus Puebla Atlixcáyotl 5718, Reserva Territorial Atlixcáyotl, Puebla 72453, Mexico

ᵈEscuela de Ingeniería y Ciencias, Tecnológico de Monterrey, Campus Estado de México, EM CP 52926 Atizapán de Zaragoza, Mexico

ᵉDipartimento di Scienze Matematiche “G. L. Lagrange”, Politecnico di Torino, Turin 10129, Italy

ᵉ'School of Mathematics and Statistics, Mathematics and Statistics Building, University of Glasgow, University Place, Glasgow G12 8QQ, UK

We address the calculation of the effective properties of non-aging linear viscoelastic composite materials. This is done by solving the microscale periodic local problems obtained via the Asymptotic Homogenization Method (AHM) by means of finite element three-dimensional simulations. The work comprises the investigation of the effective creep and relaxation behavior for a variety of fiber and inclusion reinforced structures (e.g. polymeric matrix composites). As starting point, we consider the elastic-viscoelastic correspondence principle and the Laplace-Carson transform. Then, a classical asymptotic homogenization approach for composites with discontinuous material properties and perfect contact at the interface between the constituents is performed. In particular, we reach to the stress jump conditions from local problems and obtain the corresponding interface loads. Furthermore, we solve numerically the local problems in the Laplace-Carson domain, and compute the effective coefficients. Moreover, the inversion to the original temporal space is also carried out. Finally, we compare our results with those obtained from different homogenization approaches, such as the Finite-Volume Direct Averaging Micromechanics (FVDAM) and the Locally Exact Homogenization Theory (LEHT).

## 1. Introduction

The study of creep and relaxation in viscoelastic heterogeneous materials, in which the phases generally involve both instant elastic and time-dependent viscous behavior, has gone some way towards enhancing our understanding of these kinds of composites. Several examples of man-made materials consist of viscoelastic constituents and are designed with practical purposes (see Liu, Tang, Yu, & Pipes, 2018; Dastjerdi, Akgz, & Civalek, 2020). Others are found in nature such as biological tissues (see Carpenter, Gholipour, Ghayesh, Zander, & Psaltis, 2020; Gholipour, Ghayesh, Zander, & Psaltis, 2020; Kumar, Khader, Pai, Khan, & Kyriacou, 2020; Malhotra et al., 2019). From the point of view of modeling and manufacturing,

* Corresponding author at: Escuela de Ingeniería y Ciencias, Tecnológico de Monterrey, Campus Puebla Atlixcáyotl 5718, Reserva Territorial Atlixcáyotl, 72453 Puebla, Mexico.

E-mail address: reinaldo@matcom.uh.cu (R. Rodríguez-Ramos).
the high cost associated to experimental evaluations of the effective viscoelastic properties of composites motivates the
development of predictive modeling techniques.

Recently, several studies have addressed the problem of modeling viscoelastic composites. For instance, in
Sevostianov, Levin, and Radi (2016), the calculation of the effective viscoelastic properties of a short fiber reinforced com-
posite is studied, wherein the orientation distribution of the fibers varies from perfectly aligned to randomly oriented ones.
Moreover, in Vilchevskaya, Levin, Seyedkavoosi, and Sevostianov (2019), the overall properties of viscoelastic porous mate-
rials obtained by infilling of pores with either a viscoelastic solid or fluid are investigated. In Smirnov, Vilchevskaya, and
Sevostianov (2019), the Authors studied materials containing multiple inhomogeneities of oblate shape. In particular, they
use this framework to model the viscoelastic properties of blood. In these three works, the viscoelastic properties are ex-
pressed in terms of the fraction-exponential operators of Scott BlairRabotnov (Sevostianov, Mogilevskaya, & Kushch, 2019).
On the other hand, in Ebrahimi and Barati (2016), a nonlocal higher-order refined magneto-electro-viscoelastic beam model
for vibration analysis of smart nanostructures is considered. Moreover, a RVE-based Finite Element algorithm for evaluating
the effective viscoelastic creep behaviors of aligned short fiber composites is developed in Wang and Smith (2019). Further-
more, in Nguyen, Lee, Han, and Cho (2020), the higher-order zigzag theory is used to analyze the hygrothermo-mechanical
behavior of viscoelastic composite laminated plates. In addition, a thermodynamic framework for viscoelastic polymers un-
dergoing microstructural changes is studied in Song, Muliana, and Rajagopal (2019), and an approach for predicting the
five independent viscoelastic constants of a unidirectional carbon fibre epoxy resin composite is experimentally validated in
Hine and Gusev (2019). Besides, in Rodríguez-Ramos et al. (2020), the two-scale asymptotic homogenization method is ex-
tended to the viscoelastic case and analytical closed form expressions, in Laplace-Carson space, for the effective coefficients
of fibrous viscoelastic composites with square and hexagonal cells are reported. In Farajpour, Ghayesh, and Farokhi (2019),
the effects of viscoelasticity as well as geometrical imperfections on the mechanics of carbon nanotubes with large de-
formations are analysed via a modified nonlocal elasticity model. Finally, in Omari, Almagableh, Sevostianov, Ashhab, and
Yaseen (2020) and Rahmanpanah, Mouloodi, Burvill, Gohari, and Davies (2020), the Authors use an artificial neural network
to model the viscoelastic properties of a thermoset vinyl ester nanocomposite and predict the load-displacement curve in a
long bone, respectively.

The AHM has proved to be advantageous in the description of the coarse scale mechanics of composite materials. Many
studies have dealt with the theoretical bases of the method, see e.g., Bensoussan, Papanicolau, and Lions (1978), Sanchez-
Palencia (1980), Bakhvalov and Panasenko (1989), Allaire and Briane (1996), Auriault, Boutin, and Geindreau (2009). In gen-
eral, the AHM permits to obtain an effective characterization of the heterogeneous system or phenomenon under study
by encoding the information available at the microscale into the so-called effective coefficients. This, in turn, reduces dra-
matically the computational complexity of the resulting boundary problems. However, the main disadvantage of AHM is
that the analytical solution of the local problems has been derived for only a few composite structures (see Guinovart-
Díaz et al. (2011), Willoughby, Parnell, Hazel, & Abrahams, 2012, Kalamkarov, Andrianov, Pacheco, Savi, & Starushenko, 2016).
For instance, in Cruz-González et al. (2020) and Otero et al. (2020), the Authors report a set of formulae for laminated and
fiber reinforced viscoelastic composites, respectively. Therefore, in order to handle complex microstructures, numerical ap-
proaches based on the finite element method (FEM) offer a potential alternative to solve the local problems.

The aim of this work is to apply a semi-analytical technique that combines the theoretical strengths of the AHM with
numerical computations based on FEM. The novelty of this study lies in the computations of the effective properties of
non-aging linear viscoelastic composites in three dimensions. This is done via suitable generalization of the theoretical and
computational platform presented in Penta and Gerisch (2016, 2017) for linear elastic composites. To the best of our knowl-
edge, there are no studies addressing numerical computations of the properties of viscoelastic composites by varying the
microstructure within the periodic cell.

In addition, the manuscript presents some differences with respect to previous works of some of the authors (Otero
et al., 2020; Penta & Gerisch, 2016; Rodríguez-Ramos et al., 2020) and with others in the literature (Chen, Wang, Chen, &
Geng, 2017; Chen et al., 2020; Li, Chen, Liu, & Wang, 2019; Wang & Pindera, 2016a). Specifically, the numerical solution of
the local problem is tackled by means of the finite element software COMSOL Multiphysics® and LiveLink™ for Matlab®
scripting, and the solution's convergence is analyzed through three types of mesh discretization. We compute the effective
coefficients for several types of fiber and inclusion reinforced structures with isotropic components and perfect contact at
the interfaces. In particular, Maxwell's model, four parameter model, as well as Prony's series are employed to describe the
viscoelastic effects in the sub-phases and the results are tabulated in Section 4.5.

The manuscript is organized as follows. In Section 2, we give a brief overview of the linear viscoelastic heterogeneous problem. In Section 3, we describe the mathematical and computational methodology. In Section 4, we compare our numerical results with different methods existing in the literature. By doing this, we offer an alternative way to compute the effective coefficients of non-aging linear viscoelastic composites.

## 2. Linear viscoelastic heterogeneous problem

Let us denote by $\Omega \in \mathbb{R}^3$ a linear, viscoelastic composite material with well-separated length scales $\ell$ and $L$. Particularly, $\ell$ refers to the characteristic length scale of the microstructure, while $L$ is the characteristic length of the composite medium
<figure><img src="image_2.png" /><figcaption>Fig. 1. (a) Macroscale: viscoelastic heterogeneous material. (b) ε-structural level. (c) Microscale: periodic cell. The inclusions do not intersect the boundaries and the matrix and the sub-phases are perfectly bonded.</figcaption></figure>

(see Fig. 1). Then, we introduce the dimensionless, scaling parameter

$$ \varepsilon := \frac{\ell}{L} \ll 1. \qquad (1) $$

Moreover, by denoting with x the macroscopic, or slow, variable, and using Eq. (1), we define the microscopic, or fast, variable y as

$$ y := \frac{x}{\varepsilon}. \qquad (2) $$

The heterogeneous material Ω is assumed to be a two-phase composite where $\bar{\Omega} = \bar{\Omega}_1^\varepsilon \cup \bar{\Omega}_2^\varepsilon$ and $\bar{\Omega}_1^\varepsilon \cap \Omega_2^\varepsilon = \Omega_1^\varepsilon \cap \bar{\Omega}_2^\varepsilon = \emptyset$. In particular, we consider $\Omega_2^\varepsilon = \cup_{i=1}^N \Omega_2^\varepsilon$, where N represents the number of inclusions, and the interface between $\Omega_1^\varepsilon$ and $\Omega_2^\varepsilon$ is denoted by $\Gamma^\varepsilon$ (see Fig. 1(b)). Furthermore, Y stands for the unitary periodic cell and its constituents, $Y_1$ and $Y_2$, satisfy the constraint $\bar{Y} = \bar{Y}_1 \cup \bar{Y}_2$ and $\bar{Y}_1 \cap Y_2 = Y_1 \cap \bar{Y}_2 = \emptyset$. Finally, $\Gamma_Y$ represents the interface between $Y_1$ and $Y_2$ (see Fig. 1(c)).

We remark that in the remainder of this work, the notation $\phi^\varepsilon(x, t) = \phi(x, y, t)$ is considered, where $\phi$ is assumed to be periodic in y.

## 2.1. Statement of the problem

Here, we consider that the constitutive response of the constituents of Ω to be linear viscoelastic. Then, the balance of linear momentum, neglecting inertial terms, reads

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

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

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

Boundary conditions on $\partial\Omega \times \mathbb{R}$,

$$ \mathbf{u}^{\varepsilon}(x, t) = \tilde{\mathbf{u}}, \qquad (3d) $$

$$ \boldsymbol{\sigma}^{\varepsilon}(x, t) \cdot \mathbf{n} = \tilde{\mathbf{S}}. \qquad (3e) $$

Initial condition in $\Omega \times \{0\}$,

$$ \mathbf{u}^{\varepsilon}(x, t) = \mathbf{0}. \qquad (3f) $$

In Eqs. (3a)-(3f), $\boldsymbol{\sigma}^\varepsilon$ represents the second-order stress tensor and it is given by the expression, Christensen (1982)

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

where the term $\mathbf{u}^\varepsilon$ stands for the displacement field and $\mathcal{R}^\varepsilon$ refers to the relaxation modulus, which is assumed to be a smooth function of x in $(\Omega \setminus \Gamma^\varepsilon) \times \mathbb{R}$. In general, however, we consider discontinuous material properties on $\Gamma^\varepsilon \times \mathbb{R}$.
Moreover, the symmetry properties $R_{ijkl}^{\epsilon} = R_{jikl}^{\epsilon} = R_{ijlk}^{\epsilon} = R_{klij}^{\epsilon}$ hold true. Furthermore, $\xi$ represents the second-order strain tensor and is defined by the formula

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

We notice that in Eqs. (3b)-(3c), we impose continuity conditions for displacements and tractions on $\Gamma^{\epsilon}$, in which the term $\mathbf{n}^{(y)}$ stands for the outward unit vector to the surface $\Gamma^{\epsilon}$ and the operator $\llbracket \phi^{\epsilon} \rrbracket$ denotes the “jump” of $\phi^{\epsilon}$ across the interface between the constituents.

Since the scale-dependent constitutive law (4) corresponds to a non-aging linear viscoelastic material, we can employ the elastic-viscoelastic correspondence principle (see Lakes, 2009; Maghous & Creus, 2003 for further details) and state the problem (3a)-(3f) in the Laplace-Carson space. Before doing this, we introduce the Laplace-Carson transform, for all real numbers $t \ge 0$, of $\phi^{\epsilon}$ namely, Christensen (1982)

$$ \hat{\phi}^{\epsilon}(x, p) := p \int_{0}^{\infty} e^{-pt} \phi^{\epsilon}(x, t) dt, $$

where the “hat” symbol denotes the Laplace-Carson transform of $\phi^{\epsilon}$ and $p$ represents the variable in the Laplace-Carson space.

Then, taking into account the above considerations, we rewrite the system (3a)-(3f) in the Laplace-Carson space, as follows

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

$$ \llbracket \hat{\mathbf{u}}^{\epsilon}(x, p) \rrbracket = \mathbf{0}_{\mathrm{on}} \Gamma^{\epsilon} \times [0, +\infty), \qquad (6b) $$

$$ \llbracket \hat{\mathcal{R}}^{\epsilon}(x, p) : \xi(\hat{\mathbf{u}}^{\epsilon}(x, p)) \cdot \mathbf{n}^{(y)} \rrbracket = \mathbf{0}_{\mathrm{on}} \Gamma^{\epsilon} \times [0, +\infty). \qquad (6c) $$

Boundaryconditionson $\partial\Omega \times [0, +\infty)$,

$$ \hat{\mathbf{u}}^{\epsilon}(x, p) = \bar{\mathbf{u}}, \qquad (6d) $$

$$ \hat{\mathcal{R}}^{\epsilon}(x, p) : \xi(\hat{\mathbf{u}}^{\epsilon}(x, p)) \cdot \mathbf{n} = \bar{\mathbf{S}}. \qquad (6e) $$

Initialconditionin $\Omega \times \{0\}$,

$$ \hat{\mathbf{u}}^{\epsilon}(x, p) = \mathbf{0}. \qquad (6f) $$

### 3. Description of the semi-analytical technique

In this section, we describe the semi-analytical technique adopted in this work. Specifically, we highlight the main results of the AHM, and show the procedure for the computation of the effective properties of the composite. It is worth to recall that, in this study, we adapt the methodology described in Penta and Gerisch (2016, 2017) to our purposes, and extend the results to the case of non-aging, linear viscoelastic composite materials.

#### 3.1. Two-scale asymptotic homogenization approach

We notice that, according to the chain rule and Eq. (2), the following relation is obtained

$$ \frac{\partial \hat{\phi}_{i}^{\epsilon}(x, p)}{\partial x_{j}} = \frac{\partial \hat{\phi}_{i}(x, y, p)}{\partial x_{j}} + \frac{1}{\epsilon} \frac{\partial \hat{\phi}_{i}(x, y, p)}{\partial y_{j}}. \qquad (7) $$

Therefore, 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 (8) $$

with

$$ \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 (9) $$

Now, we propose the solution of the heterogeneous problem (6a)-(6f) as a formal multi-scale expansion of the displacement $\hat{\mathbf{u}}^{\epsilon}$ in powers of $\epsilon$, Bakhvalov and Panasenko (1989), Cioranescu and Donato (1999)

$$ \hat{\mathbf{u}}^{\epsilon}(x, p) = \sum_{l=0}^{+\infty} \epsilon^{l} \hat{\mathbf{u}}^{(l)}(x, y, p). \qquad (10) $$
After replacing (10) into (6a)-(6f) and using the relations (7) and (8), the procedure follows by grouping the resulting equations in powers of the smallness parameter $\epsilon$. Specifically, by performing some algebraic manipulations, we reach the following sequence of problems, (we refer the reader to Rodríguez-Ramos et al. (2020), Cruz-González et al. (2020) for further details)

**Problem for $\epsilon^{-2}$**

$$ \frac{\partial}{\partial y_j} [\hat{\mathcal{R}}_{ijkl}(y, p) \xi_{kl}^{(y)}(\hat{\mathbf{u}}^{(0)}(x, y, p))] = 0 \quad \text{in } (Y \setminus \Gamma_Y) \times [0, +\infty), \quad (11a) $$

$$ [\hat{\mathbf{u}}_i^{(0)}(x, y, p)] = 0 \quad \text{in } \Gamma_Y \times [0, +\infty), \quad (11b) $$

$$ [\hat{\mathcal{R}}_{ijkl}(y, p) \xi_{kl}^{(y)}(\hat{\mathbf{u}}^{(0)}(x, y, p)) n_j^{(y)}] = 0 \quad \text{in } \Gamma_Y \times [0, +\infty). \quad (11c) $$

Initial condition in $Y \times \{0\}$,

$$ \hat{\mathbf{u}}_i^{(0)}(x, y, p) = 0. \quad (11d) $$

Since the right hand side of (11a) is zero, the problem can be solved using the solvability condition reported in Bakhvalov and Panasenko (1989), Persson, Persson, Svanstedt, & Wyller, 1993, wherein $\hat{\mathbf{u}}^{(0)}$ is a solution of (11a)-(11c) if and only if it is constant in relation to the variable $y$, i.e.

$$ \hat{\mathbf{u}}^{(0)}(x, y, p) = \hat{\mathbf{v}}(x, p), \quad (12) $$

where $\hat{\mathbf{v}}$ is a smooth function.

**Problem for $\epsilon^{-1}$**

$$ \frac{\partial}{\partial y_j} [\hat{\mathcal{R}}_{ijkl}(y, p)(\xi_{kl}^{(y)}(\hat{\mathbf{u}}^{(1)}(x, y, p)) + \xi_{kl}(\hat{\mathbf{u}}^{(0)}(x, y, p)))] = 0 \quad \text{in } (Y \setminus \Gamma_Y) \times [0, +\infty), \quad (13a) $$

$$ [\hat{\mathbf{u}}_i^{(1)}(x, y, p)] = 0 \quad \text{in } \Gamma_Y \times [0, +\infty), \quad (13b) $$

$$ [\hat{\mathcal{R}}_{ijkl}(y, p)(\xi_{kl}^{(y)}(\hat{\mathbf{u}}^{(1)}(x, y, p)) + \xi_{kl}(\hat{\mathbf{u}}^{(0)}(x, y, p)))n_j^{(y)}] = 0 \quad \text{in } \Gamma_Y \times [0, +\infty). \quad (13c) $$

Initial condition in $Y \times \{0\}$,

$$ \hat{\mathbf{u}}_i^{(1)}(x, y, p) = 0. \quad (13d) $$

The solvability condition also applies in this case, ensuring the existence and uniqueness, up to a $y$-constant function, of the solution to the problem (13a)-(13d). An expression for the aforementioned solution can be derived by means of the method of separation of variables, namely (Bakhvalov & Panasenko, 1989; Persson, Persson, Svanstedt, & Wyller, 1993),

$$ \hat{\mathbf{u}}_m^{(1)}(x, y, p) = \hat{\chi}_{klm}(y, p) \xi_{kl}(\hat{\mathbf{v}}(x, p)), \quad (14) $$

where the third order tensor $\hat{\chi}_{klm}$, which is periodic in $y$, is the solution of the following problem referred to as the $\epsilon$-local problem,

$$ \frac{\partial}{\partial y_j} [\hat{\mathcal{R}}_{ijpq}(y, p) \xi_{pq}^{(y)}(\hat{\mathbf{X}}_{kl}(y, p)) + \hat{\mathcal{R}}_{ijkl}(y, p)] = 0 \quad \text{in } (Y \setminus \Gamma_Y) \times [0, +\infty), \quad (15a) $$

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

$$ [\langle \hat{\mathcal{R}}_{ijpq}(y, p) \xi_{pq}^{(y)}(\hat{\mathbf{X}}_{kl}(y, p)) + \hat{\mathcal{R}}_{ijkl}(y, p) \rangle n_j^{(y)}] = 0 \quad \text{in } \Gamma_Y \times [0, +\infty). \quad (15c) $$

Initial condition in $Y \times \{0\}$,

$$ \hat{\chi}_{klm}(y, p) = 0. \quad (15d) $$

Where,

$$ \xi_{pq}^{(y)}(\hat{\mathbf{X}}_{kl}(y, p)) := \frac{1}{2} \left( \frac{\partial \hat{\chi}_{klp}(y, p)}{\partial y_q} + \frac{\partial \hat{\chi}_{klq}(y, p)}{\partial y_p} \right). \quad (16) $$

According to the problem (15a)-(15d), a further condition is required to obtain uniqueness of the solution, for instance, by requiring that $\langle \hat{\chi}_{klm} \rangle_y = 0$ (Penta & Gerisch, 2016; 2017) or alternately, by fixing the value of $\hat{\chi}_{klm}$ at one point of the reference cell $Y$. Here, the notation $\langle \phi \rangle_y$ stands for the cell average operator and is defined as,

$$ \langle \phi \rangle_y := \frac{1}{|Y|} \int_Y \phi \, dy, \quad (17) $$

where $|Y|$ represents the volume of the periodic cell $Y$.

**Problem for $\epsilon^0$**
Now, by considering relation (14), the homogenized problem can be written in the form

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

Boundaryconditionson $\partial\Omega^h \times [0, +\infty)$,

$$ \hat{v}_i(x, p) = \bar{u}_i, \quad (18b) $$

$$ \hat{\mathcal{R}}_{ijkl}^{(*)}(p) \xi_{kl}(\hat{\mathbf{v}}(x, p)) n_j = \bar{S}_i. \quad (18c) $$

Initialconditionin $\Omega^h \times \{0\}$,

$$ \hat{v}_i(x, p) = 0. \quad (18d) $$

In Eq. (18a), $\hat{\mathcal{R}}^{(*)}$ represents the effective relaxation modulus in the Laplace-Carson space and is defined through the expression

$$ \hat{\mathcal{R}}_{ijkl}^{(*)}(p) := \langle \hat{\mathcal{R}}_{ijkl}(y, p) + \hat{\mathcal{R}}_{ijpq}(y, p) \xi_{pq}^{(y)}(\hat{\mathbf{X}}_{kl}(y, p)) \rangle_y. \quad (19) $$

## 3.2. Computational approach

In this section, we report some of the main aspects of the computational approach used in this work. In particular, for simplicity in our computations, we assume a y-constant expression for the relaxation modulus ($\hat{\mathcal{R}}_{ijkl}$) in each phase of the periodic cell, i.e.

$$ \hat{\mathcal{R}}_{ijpq}(y, p) = \begin{cases} \hat{\mathcal{R}}_{ijpq}^{(1)}(p), & \text{if } y \in Y_1, \\ \hat{\mathcal{R}}_{ijpq}^{(2)}(p), & \text{if } y \in Y_2, \end{cases} \quad (20) $$

where the superscripts (1) and (2) indicate the matrix and the inclusion phases, respectively (see Fig. 1). Therefore, the local problem (15a)-(15d) can be rewritten as,

$$ \frac{\partial}{\partial y_j} \left[ \hat{\mathcal{R}}_{ijpq}^{(1)}(p) \xi_{pq}^{(y)}(\hat{\mathbf{X}}_{kl}^{(1)}(y, p)) \right] = 0 \quad \text{in } Y_1 \times [0, +\infty), \quad (21a) $$

$$ \frac{\partial}{\partial y_j} \left[ \hat{\mathcal{R}}_{ijpq}^{(2)}(p) \xi_{pq}^{(y)}(\hat{\mathbf{X}}_{kl}^{(2)}(y, p)) \right] = 0 \quad \text{in } Y_2 \times [0, +\infty), \quad (21b) $$

$$ \hat{\mathbf{X}}_{klm}(y, p) = \hat{\mathbf{X}}_{klm}^{(2)}(y, p) \quad \text{on } \Gamma_Y \times [0, +\infty), \quad (21c) $$

$$ \begin{aligned} & (\hat{\mathcal{R}}_{ijpq}^{(1)}(p) \xi_{pq}^{(y)}(\hat{\mathbf{X}}_{kl}^{(1)}(y, p))) n_j^{(y)} - (\hat{\mathcal{R}}_{ijpq}^{(2)}(p) \xi_{pq}^{(y)}(\hat{\mathbf{X}}_{kl}^{(2)}(y, p))) n_j^{(y)} \\ & = (\hat{\mathcal{R}}_{ijkl}^{(2)}(p) - \hat{\mathcal{R}}_{ijkl}^{(1)}(p)) n_j^{(y)} \quad \text{on } \Gamma_Y \times [0, +\infty). \end{aligned} \quad (21d) $$

Initial condition in $Y \times \{0\}$,

$$ \hat{\mathbf{X}}_{klm}(y, p) = 0. \quad (21e) $$

We notice that, in Eq. (21d), the interface loads $f_{kl}^{(y)} := (\hat{\mathcal{R}}_{ijkl}^{(2)} - \hat{\mathcal{R}}_{ijkl}^{(1)}) n_j^{(y)}$ represent the driving forces leading to nontrivial solutions of the six elastic-type local problems ($(k, l), k \ge l$) when discontinuities in the coefficients (20), between the matrix and inclusions, arise (see Penta & Gerisch, 2016; Penta & Gerisch, 2017).

Now, in order to compute the numerical solution of the six elastic-type local problems (21a)-(21e) and the effective viscoelastic properties of the given composite material, we employ a combination of the finite element software COMSOL Multiphysics® with LiveLink™ for Matlab® scripting. The procedure can be summarized as follows:

(I) To define the geometry of the periodic cell Y, introduce the parameters, set the properties of the materials and implement the six elastic-type local problems (21a)-(21e) in COMSOL Multiphysics®. In this first step, the scheme proposed in Penta and Gerisch (2016) is adapted for our purposes.

(II) To integrate into a Matlab script the codes for working with LiveLink™ for Matlab®, the viscoelastic material properties and the algorithm for the numerical inversion of the Laplace-Carson transform.

(III) To solve, for every value of p, the set of elastic local problems (21a)-(21e) in Laplace-Carson space and to calculate the effective relaxation modulus $\hat{\mathcal{R}}^{(*)}$ using (19), updating in each iteration the corresponding parameters and the properties of the materials.

(IV) To invert to the original time domain. The function INVLAP, developed in Juraj (2020) is used for this purpose. The algorithm needs to calculate 40 points in Laplace-Carson's space in order to find one point for the effective coefficient in the time domain.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Three types of mesh discretization for the periodic cells: Mesh A, B and C.</figcaption></figure>

For completeness of the analysis in the calculation of the effective creep compliance $J^{(*)}$, the mathematical relationship between both magnitudes $R^{(*)}$ and $J^{(*)}$ is presented in Laplace-Carson space as follows,

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

where $I$ is the fourth order identity tensor (Hashin, 1972). We notice that expression (22) is used to compute $J^{(*)}$ once $\hat{R}^{(*)}$ is known.

# 4. Numerical results and discussion

In this section, we compute the effective properties of a composite for several periodic cell configurations, report the results of our computations and compare them with other approaches found in the literature. By doing this, we validate the methodology described in Section 3 and highlight its major advantages. We notice that in the remainder of this manuscript the acronym PM stands for present model.

To start with, in Sub-Section 4.1 we address the instant elastic response of a composite material reinforced with long circular fibers aligned along the axis $Ox_3$ (see Fig. 9(a)), and we study the solution's convergence through three types of mesh discretization. Moreover, in Sub-Section 4.2, we compute the effective relaxation modulus and the effective creep compliance for unidirectionally-reinforced glass/epoxy composites. Specifically, the four parameter model is employed to describe the viscoelastic polymeric matrix. In Sub-Section 4.3, we compute the effective relaxation modulus for viscoelastic composites reinforced with spherical inclusions. In this case, the matrix viscoelastic properties are expressed in terms of a Maxwell's model. In Sub-Section 4.4, we consider three different inclusion shapes for the calculation of the effective relaxation modulus. Here, the Prony's series are used to describe both viscoelastic constituents. Finally, in Sub-Section 4.5, we take into account all the inclusions shapes displayed in Fig. 9 for the computation of the effective relaxation modulus.

## 4.1. Instant elastic response

We start our analysis by considering a composite with periodically distributed and non-overlapping long fibers with circular cross-sections and perfectly bonded to the matrix (see Fig. 9(a)). Besides, we assume that the both constituents of the heterogeneous material are linear elastic and isotropic. We notice that, by referring to Eq. (4), the instant elastic response of the composite can be obtained by considering $t = 0$. In this limit case, the composite describe a linear elastic material, and $R_{ijkl}^{(*)}$ will represents the effective elastic tensor.

In the present context, we determine the instant elastic effective response in a wide range of volume fractions, and assume that $v_1 = 0.39$, $v_2 = 0.2$. Particularly, we focus on the porous ($\mu_2 = 0$) and moderate contrast ($\mu_2/\mu_1 = 6$) cases. We notice that the notation $v_1(v_2)$ stands for Poisson's ratio and $\mu_1(\mu_2)$ for the shear modulus of the matrix (fibers). The values of these parameters have been obtained from Mogilevskaya, Kushch, and Nikolskiy (2014).

It is worth to remark that, here, we keep track of the solution's convergence through three types of mesh discretization for the periodic cells (see Fig. 2). In particular, the mesh of the periodic cell $Y$ is built such that, first, we create surface meshes of free triangular type for the interface $\Gamma_Y$, and then, the latter are extended into a three-dimensional mesh of free tetrahedral type, which covers the whole periodic cell (see Fig. 2). The increasing mesh refinements are described in Table 1. We notice that the mesh on the interface between the matrix and the inclusion is more refined to better capture the boundary loads contribution. The mesh statistic outlines that, for instance, the number of tetrahedra subvolumes for a volume fractions of 0.3 are 1015 (Mesh A), 8460 (Mesh B) and 37 067 (Mesh C), respectively.

The values of the normalized effective elastic coefficients, denoted by $(R_{11}^{(*)} - R_{12}^{(*)})/2\mu_0$, $R_{55}^{(*)}/\mu_0$, $R_{13}^{(*)}/R_{13}^{(1)}$ and $R_{33}^{(*)}/R_{33}^{(1)}$, as functions of the volume fraction ($V_f$), are reported in Fig. 3(a)-(d) for both porous and moderate contrast
<figure><img src="image_4.png" /></figure>

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

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

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

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

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

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

<figure><img src="image_11.png" /><figcaption>**Fig. 3.** Charts (a)–(d) show the calculation of the instant elastic response for normalized effective coefficients. The comparisons are carried out with Tables 1–5 of Mogilevskaya et al. (2014). In addition, charts (e)–(h) display the relative errors with respect to RUC, for both porous and moderate contrast of the three mesh discretizations.</figcaption></figure>
<table><caption>**Table 1**

Mesh discretization of the unit cells.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Element Size</th></tr><tr><th>Interface (Γ<sub>γ</sub>)</th><th>Remaining (Y<sub>1</sub>Γ<sub>γ</sub>)</th></tr></thead><tbody><tr><td>Mesh A</td><td>Coarser</td><td>Extremely Coarse</td></tr><tr><td>Mesh B</td><td>Fine</td><td>Coarse</td></tr><tr><td>Mesh C</td><td>Finer</td><td>Fine</td></tr></tbody></table>

cases. These figures display the results obtained by the present model (PM), for each of the meshes A, B and C, with respect to the volume fraction. Moreover, in Fig. 3 we compare with different sets of simple, closed-form and approximate formulae given by Kantor and Bergman (K&B), Luciano and Barbero (L&B), non-interacting Maxwell's approximations (non-int) and exact periodic solutions (RUC) Mogilevskaya et al. (2014).

A closer look at the data in Fig. 3(a)–(d) reveals that there is good agreement in the comparisons between PM and RUC and no significant differences are found in relation to K&B and L&B. The most notable discrepancy arises from the data corresponding to the non-int approach in Fig. 3(a), nonetheless the remaining coefficients are consistent within the comparisons. Additionally, we note the three types of mesh discretizations are barely distinguishable.

The relative errors with respect to RUC, for both porous and moderate contrast of the three mesh discretizations are shown in Fig. 3(e)–(h). The results provide further evidence of the good correlation between PM and RUC. A remarkable pattern is that with the increased mesh refinement the relative error is improved, however, even with the coarser mesh (Mesh A) the maximum observed error is less than 3%, highlighting the good convergence of the present model.

Then, based on the good performance of the three mesh discretizations and the execution time, the remaining calculations are performed using Mesh B (see Fig. 2).

## 4.2. Viscoelastic response: four-parameter model

Now, the potential of PM in the calculation of the effective relaxation modulus and the effective creep compliance for unidirectionally-reinforced glass/epoxy composites (see Fig. 9(a)) is shown. Particularly, the linear isotropic constituents of the heterogeneous structure are elastic glass fibers embedded into a viscoelastic polymeric matrix. The latter is described using the four-parameter model or Burger's model (see Mainardi & Spada, 2011), where the expression of the creep compliance is given as follows,

$$J(t) = \frac{1}{E_0} + \frac{1}{E_1} \left(1 - \exp\left(-\frac{E_1}{\eta_1}t\right)\right) + \frac{t}{\eta_2} \quad (23)$$

The material constants for this model are $E_0 = 3.27$ GPa, $E_1 = 1.8$ GPa, $\eta_1 = 300$ GPa·h and $\eta_2 = 8000$ GPa·h, where $E$ and $\eta$ represent the elastic modulus of spring and the constitutive parameter of dashpot, respectively. Furthermore, the matrix Poisson's ratio is assumed to be $\nu_1 = 0.38$, and the mechanical properties of the fibers are $E_2 = 68.77$ GPa and $\nu_2 = 0.21$. The data has been taken from Cavalcante and Marques (2014), Chen et al. (2017).

In Fig. 4, we compare the results of PM with those obtained by the Locally exact homogenization theory (LEHT) reported in Wang and Pindera (2016a,b), the finite-volume direct averaging micromechanics theory (FVDAM 2D) proposed in Cavalcante and Marques (2014), and the three-dimensional FVDAM developed in Chen et al. (2017). These results correlate fairly well in the comparison against LEHT and FVDAM for three fiber volume fractions 20%, 30% and 40%. Additionally, Fig. 5 shows the relative errors of PM with respect to both FVDAM and LEHT. We can notice that the maximum error is less than 1%.

Finally, in Fig. 6 we report the values of the effective creep compliance $J^{(*)}$ and compare them with those resulting from FVDAM 2D. The data for FVDAM 2D were collected from Fig. 12 of Cavalcante and Marques (2014) by using the free software G3Data.

## 4.3. Viscoelastic response: Maxwell's model

We now consider a spherical reinforced composite with repeating unit cell of the type shown in Fig. 9(h), and analyze the overall viscoelastic behavior under the assumption of isotropic phases. In particular, the spherical inclusions are assumed to be linear elastic with Young's modulus and Poisson's ratio being $E_2 = 168.4$ GPa and $\nu_2 = 0.443$, respectively, whereas the matrix is described as a linear viscoelastic constituent and the Maxwell's model (see Lakes, 2009; Mainardi & Spada, 2011) is used to simulate its properties, namely

$$S(t) = \frac{1}{E_0} + \frac{t}{\eta_0} \quad (24)$$

In the above equation, $E_0 = 4.082$ GPa and $\eta_0 = 159.81$ GPa·min. Moreover, the matrix Poisson's ratio is $\nu_1 = 0.311$. The input data can be found in Chen et al. (2017).
<figure><img src="image_12.png" /><figcaption>**Fig. 4.** Computation of the effective relaxation modulus for unidirectionally-reinforced glass/epoxy composites and the comparison against LEHT and FV-DAM. The comparisons are made with Fig. 11 of Chen et al. (2017).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 5. The estimation of relative errors of PM with respect to both FVDAM and LEHT for three different volume fractions.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 6. Computation of the effective creep compliance. The comparisons are performed with Fig. 12 of Cavalcante and Marques (2014).</figcaption></figure>
<figure><img src="image_15.png" /></figure>

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

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

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

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

<figure><img src="image_20.png" /><figcaption>**Fig. 7.** Charts (a)–(c) show the calculation of the effective relaxation modulus for spherical reinforced composites. Our results are compared with Fig. 6 of Chen et al. (2017). Moreover, charts (d)–(f) present the relative errors of PM with respect to FVDAM and FEM.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. 8. Charts (a)-(b) display comparisons with Figs. 2 and 3 of Kouri et al. (2016). Charts (c)-(d) present the results of a numerical experiment for the fiber's volumetric fraction equal to zero.</figcaption></figure>

As illustrated in Fig. 7(a)-(c), there is a good agreement between PM, FVDAM and FEM in the calculation of the effective relaxation moduli using three different volume fractions (10%, 20%, 30%). Moreover, due to low relative errors in relation to FVDAM (see Fig. 7(d)-(f)), we can confirm the complete consistency between PM and FVDAM. In contrast, as time passes, much higher values for the relative errors with respect to FEM are reported. A possible explanation for this is tackled in Chen et al. (2017), wherein it is suggested that the values of the effective relaxation moduli become small after 50 min.

### 4.4. Viscoelastic response: Prony's series

In this section, we take inspiration from Kouri, Bakkali, and Azrar (2016) and consider three different inclusion shapes for the calculation of the effective relaxation modulus i.e. fibrous, spherical and ellipsoidal (see Fig. 9(a),(h) and (i), respectively). Moreover, both constituents are assumed isotropic and linear viscoelastic.

The shear and bulk relaxation functions $G(t)$ and $K(t)$ can be computed in terms of Prony's series as follows,

$$G(t) = G_{\infty} + \sum_{i=1}^{N_G} G_i \exp\left(-\frac{t}{\xi_i^{(G)}}\right), \qquad (25a)$$
<figure><img src="image_22.png" /><figcaption>**Fig. 9.** Representation of the analyzed periodic cells. The acronyms LCF, LEF, LRF, LSF stands for long circular, elliptical, rhombic and square fibers, respectively. In addition, we consider Two LSF and Three LSF to denote two and three long square fibers and CI, SI, EI for cubic, spherical and ellipsoidal inclusions, respectively.</figcaption></figure>

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

Input data for Prony's series (see Table 1 of Kouri et al. (2016)).</caption><thead><tr><th colspan="4">Stiff material</th><th colspan="4">Soft material</th></tr><tr><th>Relaxation times (s)</th><th>Shear moduli (GPa)</th><th>Bulk moduli (GPa)</th><th>Asymptotic moduli (GPa)</th><th>Relaxation times (s)</th><th>Shear moduli (GPa)</th><th>Bulk moduli (GPa)</th><th>Asymptotic moduli (GPa)</th></tr></thead><tbody><tr><td>ξ<sub>i</sub><sup>(G)</sup></td><td>ξ<sub>j</sub><sup>(K)</sup></td><td>G<sub>i</sub></td><td>K<sub>j</sub></td><td>ξ<sub>i</sub><sup>(G)</sup></td><td>ξ<sub>j</sub><sup>(K)</sup></td><td>G<sub>i</sub></td><td>K<sub>j</sub></td></tr><tr><td>3</td><td>10000</td><td>3.162</td><td>40000</td><td>0.032</td><td>100</td><td>2.512</td><td>3000</td></tr><tr><td>10</td><td></td><td>17.1783</td><td></td><td>0.100</td><td>316.228</td><td>10</td><td>100</td></tr><tr><td>32</td><td></td><td>100</td><td></td><td>0.316</td><td></td><td>56.234</td><td></td></tr><tr><td>100</td><td></td><td>316.228</td><td></td><td>1</td><td></td><td>316.228</td><td></td></tr><tr><td>316</td><td></td><td>1000</td><td></td><td>3.162</td><td></td><td>1000</td><td></td></tr><tr><td>1000</td><td></td><td>5623.413</td><td></td><td>10</td><td></td><td>199.526</td><td></td></tr><tr><td>3162</td><td></td><td>10000</td><td></td><td>31.623</td><td></td><td>50.119</td><td></td></tr><tr><td>10000</td><td></td><td>562.341</td><td></td><td>100</td><td></td><td>19.953</td><td></td></tr><tr><td>31623</td><td></td><td>141.254</td><td></td><td>316.228</td><td></td><td>12.589</td><td></td></tr><tr><td>100000</td><td></td><td>56.234</td><td></td><td>1000</td><td></td><td>2.512</td><td></td></tr><tr><td>316228</td><td></td><td>17.783</td><td></td><td>3162.278</td><td></td><td>1.698</td><td></td></tr><tr><td>1000000</td><td></td><td>5.623</td><td></td><td>10000</td><td></td><td>1.202</td><td></td></tr><tr><td>3162278</td><td></td><td>3.162</td><td></td><td>31622.777</td><td></td><td>1.148</td><td></td></tr><tr><td>10000000</td><td></td><td>1.778</td><td></td><td>100000</td><td></td><td>1.096</td><td></td></tr></tbody></table>
<table><caption>**Table 3**
The input data for the sub-phases are the same as Section 4.3.</caption><thead><tr><th rowspan="2">TIME (min)</th><th colspan="10">Maxwell's model</th></tr><tr><th>R<sub>ij</sub><sup>(*)</sup> (V.F.)</th><th>LCF</th><th>LEF</th><th>LRF</th><th>LSF</th><th>Two LSF</th><th>Three LSF</th><th>CI</th><th>SI</th><th>EI</th></tr></thead><tbody><tr><td rowspan="10">0.01</td><td>R<sub>11</sub><sup>(*)</sup> (10%)</td><td>6.5782</td><td>6.7394</td><td>6.6068</td><td>6.6782</td><td>6.5821</td><td>13.3268</td><td>7.0848</td><td>6.8205</td><td>7.0177</td></tr><tr><td>R<sub>12</sub><sup>(*)</sup> (10%)</td><td>2.8986</td><td>2.8982</td><td>2.9290</td><td>2.8526</td><td>2.9972</td><td>3.6123</td><td>2.8709</td><td>2.9171</td><td>2.9199</td></tr><tr><td>R<sub>44</sub><sup>(*)</sup> (10%)</td><td>1.8818</td><td>1.8242</td><td>1.9176</td><td>1.9148</td><td>2.3344</td><td>2.0521</td><td>1.9040</td><td>1.8720</td><td>1.8504</td></tr><tr><td>R<sub>11</sub><sup>(*)</sup> (20%)</td><td>7.7657</td><td>8.2663</td><td>7.8575</td><td>8.0175</td><td>7.6637</td><td>22.9399</td><td>9.0257</td><td>8.4137</td><td>9.1976</td></tr><tr><td>R<sub>12</sub><sup>(*)</sup> (20%)</td><td>3.2669</td><td>3.2625</td><td>3.3423</td><td>3.1343</td><td>3.5389</td><td>5.7337</td><td>3.1548</td><td>3.2878</td><td>3.2896</td></tr><tr><td>R<sub>44</sub><sup>(*)</sup> (20%)</td><td>2.2842</td><td>2.1369</td><td>2.3830</td><td>2.3534</td><td>4.2381</td><td>3.0749</td><td>2.2788</td><td>2.2249</td><td>2.1643</td></tr><tr><td>R<sub>11</sub><sup>(*)</sup> (30%)</td><td>9.3610</td><td>10.6932</td><td>9.6522</td><td>9.8140</td><td>8.9964</td><td>34.6365</td><td>11.6231</td><td>10.7037</td><td>14.2783</td></tr><tr><td>R<sub>12</sub><sup>(*)</sup> (30%)</td><td>3.6610</td><td>3.6501</td><td>3.8094</td><td>3.3992</td><td>4.2584</td><td>9.0923</td><td>3.4451</td><td>3.6813</td><td>3.7241</td></tr><tr><td>R<sub>44</sub><sup>(*)</sup> (30%)</td><td>2.7962</td><td>2.4992</td><td>3.0373</td><td>2.9006</td><td>7.7179</td><td>4.9077</td><td>2.7188</td><td>2.6699</td><td>2.5231</td></tr><tr><td rowspan="10">50</td><td>R<sub>11</sub><sup>(*)</sup> (10%)</td><td>1.4305</td><td>1.4693</td><td>1.4376</td><td>1.4587</td><td>1.4346</td><td>8.4110</td><td>1.5724</td><td>1.4910</td><td>1.5416</td></tr><tr><td>R<sub>12</sub><sup>(*)</sup> (10%)</td><td>0.6246</td><td>0.6244</td><td>0.6318</td><td>0.6095</td><td>0.6459</td><td>1.2403</td><td>0.6078</td><td>0.6272</td><td>0.6271</td></tr><tr><td>R<sub>44</sub><sup>(*)</sup> (10%)</td><td>0.4136</td><td>0.3990</td><td>0.4238</td><td>0.4230</td><td>0.8066</td><td>0.5932</td><td>0.4217</td><td>0.4111</td><td>0.4056</td></tr><tr><td>R<sub>11</sub><sup>(*)</sup> (20%)</td><td>1.7012</td><td>1.8250</td><td>1.7246</td><td>1.7747</td><td>1.6832</td><td>17.5080</td><td>2.0687</td><td>1.8697</td><td>2.0849</td></tr><tr><td>R<sub>12</sub><sup>(*)</sup> (20%)</td><td>0.6983</td><td>0.6965</td><td>0.7158</td><td>0.6538</td><td>0.7593</td><td>2.9414</td><td>0.6402</td><td>0.6966</td><td>0.6928</td></tr><tr><td>R<sub>44</sub><sup>(*)</sup> (20%)</td><td>0.5105</td><td>0.4726</td><td>0.5395</td><td>0.5306</td><td>2.4171</td><td>1.3995</td><td>0.5130</td><td>0.4952</td><td>0.4796</td></tr><tr><td>R<sub>11</sub><sup>(*)</sup> (30%)</td><td>2.0726</td><td>2.4163</td><td>2.1495</td><td>2.2109</td><td>1.9914</td><td>28.6057</td><td>2.7635</td><td>2.4421</td><td>3.6273</td></tr><tr><td>R<sub>12</sub><sup>(*)</sup> (30%)</td><td>0.7697</td><td>0.7633</td><td>0.8020</td><td>0.6785</td><td>0.9108</td><td>5.7818</td><td>0.6504</td><td>0.7547</td><td>0.7350</td></tr><tr><td>R<sub>44</sub><sup>(*)</sup> (30%)</td><td>0.6368</td><td>0.5589</td><td>0.7104</td><td>0.6681</td><td>5.7110</td><td>3.0512</td><td>0.6213</td><td>0.6042</td><td>0.5656</td></tr><tr><td rowspan="10">150</td><td>R<sub>11</sub><sup>(*)</sup> (10%)</td><td>0.1439</td><td>0.1485</td><td>0.1448</td><td>0.1481</td><td>0.1451</td><td>7.1801</td><td>0.1649</td><td>0.1514</td><td>0.1578</td></tr><tr><td>R<sub>12</sub><sup>(*)</sup> (10%)</td><td>0.0620</td><td>0.0619</td><td>0.0628</td><td>0.0594</td><td>0.0640</td><td>0.6516</td><td>0.0578</td><td>0.0619</td><td>0.0618</td></tr><tr><td>R<sub>44</sub><sup>(*)</sup> (10%)</td><td>0.0423</td><td>0.0405</td><td>0.0438</td><td>0.0437</td><td>0.4135</td><td>0.2212</td><td>0.0439</td><td>0.0420</td><td>0.0413</td></tr><tr><td>R<sub>11</sub><sup>(*)</sup> (20%)</td><td>0.1732</td><td>0.1881</td><td>0.1761</td><td>0.1844</td><td>0.1725</td><td>16.1373</td><td>0.2302</td><td>0.1950</td><td>0.2250</td></tr><tr><td>R<sub>12</sub><sup>(*)</sup> (20%)</td><td>0.0683</td><td>0.0680</td><td>0.0702</td><td>0.0605</td><td>0.0744</td><td>2.2496</td><td>0.0546</td><td>0.0669</td><td>0.0656</td></tr><tr><td>R<sub>44</sub><sup>(*)</sup> (20%)</td><td>0.0536</td><td>0.0488</td><td>0.0579</td><td>0.0566</td><td>1.9354</td><td>0.9668</td><td>0.0548</td><td>0.0517</td><td>0.0497</td></tr><tr><td>R<sub>11</sub><sup>(*)</sup> (30%)</td><td>0.2146</td><td>0.2586</td><td>0.2248</td><td>0.2370</td><td>0.2068</td><td>27.0712</td><td>0.3285</td><td>0.2663</td><td>0.4798</td></tr><tr><td>R<sub>12</sub><sup>(*)</sup> (30%)</td><td>0.0730</td><td>0.0715</td><td>0.0760</td><td>0.0563</td><td>0.0885</td><td>4.9615</td><td>0.0433</td><td>0.0674</td><td>0.0574</td></tr><tr><td>R<sub>44</sub><sup>(*)</sup> (30%)</td><td>0.0688</td><td>0.0587</td><td>0.0804</td><td>0.0736</td><td>5.1674</td><td>2.5659</td><td>0.0679</td><td>0.0647</td><td>0.0597</td></tr></tbody></table>

$$K(t) = K_{\infty} + \sum_{j=1}^{N_K} K_j \exp\left(-\frac{t}{\xi_j^{(K)}}\right), \qquad (25b)$$

where $N_G$ and $N_K$ are the number of terms for the corresponding modulus, $\xi_i^{(G)}$ and $\xi_j^{(K)}$ stand for the corresponding relaxation times, $G_\infty$ and $K_\infty$ denote the asymptotic moduli and, $G_i$ and $K_j$ represent the shear elastic modulus and bulk elastic modulus, respectively. Furthermore, by applying Laplace-Carson transform to Eqs. (25a)-(25b) the following expressions are obtained

$$G(s) = G_{\infty} + \sum_{i=1}^{N_G} \frac{sG_i}{\frac{1}{\xi_i^{(G)}} + s}, \qquad (26a)$$

$$K(s) = K_{\infty} + \sum_{j=1}^{N_K} \frac{sK_j}{\frac{1}{\xi_j^{(K)}} + s}. \qquad (26b)$$

The parameters for the viscoelastic constituents used in the computations are reported in Table 2. In particular, we consider the fiber as the stiff material and the matrix as the soft one.

In Fig. 8(a)-(b), we display the effective relaxation moduli $R_{12}^{(*)}$ and $R_{44}^{(*)}$ considering PM and the asymptotic homogenization method (AHM) developed in Rodríguez-Ramos et al. (2020), Otero et al. (2020). The formulae for AHM involve two orders of truncation in the calculation of the effective properties, which are denoted by AHMT1 and AHMT2, respectively (see Bravo-Castillero, Guinovart-Díaz, Rodríguez-Ramos, Sabina, & Brenner, 2012). We notice that the curves in Fig. 8(a)-(b) show similar trends to those reported in Fig. 1 and 2 of Kouri et al. (2016) except for an order of 10. This apparent lack of correlation can be attributed to the assumed dimensions in relation to shear, bulk and asymptotic modulus sum-
<table><caption>**Table 4**
The input data for the sub-phases are the same as Section 4.2.</caption><thead><tr><th rowspan="2">TIME (hrs)</th><th colspan="10">The four-parameter model (Burgers model)</th></tr><tr><th>R<sup>(*)</sup><sub>ij</sub></th><th>LCF</th><th>LEF</th><th>LRF</th><th>LS</th><th>Two LSF</th><th>Three LSF</th><th>CI</th><th>SI</th><th>EI</th></tr></thead><tbody><tr><td rowspan="6">0.01</td><td>R<sup>(*)</sup><sub>22</sub> (10%)</td><td>6.9280</td><td>6.8556</td><td>6.9408</td><td>7.0033</td><td>10.3621</td><td>9.3053</td><td>7.2770</td><td>7.0949</td><td>7.0470</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (10%)</td><td>1.4236</td><td>1.3823</td><td>1.4481</td><td>1.4462</td><td>1.6315</td><td>1.4867</td><td>1.4416</td><td>1.4238</td><td>1.4088</td></tr><tr><td>R<sup>(*)</sup><sub>66</sub> (10%)</td><td>1.3701</td><td>1.3652</td><td>1.4097</td><td>1.3652</td><td>1.3926</td><td>1.4860</td><td>1.4416</td><td>1.4238</td><td>1.4306</td></tr><tr><td>R<sup>(*)</sup><sub>22</sub> (20%)</td><td>7.9713</td><td>7.7699</td><td>8.0122</td><td>8.1660</td><td>15.2759</td><td>13.0746</td><td>8.8627</td><td>8.4295</td><td>8.2657</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (20%)</td><td>1.7177</td><td>1.6128</td><td>1.7847</td><td>1.7648</td><td>2.5510</td><td>1.9982</td><td>1.7160</td><td>1.6865</td><td>1.6442</td></tr><tr><td>R<sup>(*)</sup><sub>66</sub> (20%)</td><td>1.5665</td><td>1.5561</td><td>1.6536</td><td>1.5568</td><td>1.6254</td><td>1.9984</td><td>1.7160</td><td>1.6865</td><td>1.7071</td></tr><tr><td rowspan="6">100</td><td>R<sup>(*)</sup><sub>22</sub> (30%)</td><td>9.3467</td><td>8.9045</td><td>9.4856</td><td>9.7073</td><td>20.8178</td><td>17.4240</td><td>10.9678</td><td>10.3008</td><td>9.8516</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (30%)</td><td>2.0884</td><td>1.8788</td><td>2.2494</td><td>2.1589</td><td>4.1179</td><td>2.8360</td><td>2.0357</td><td>2.0120</td><td>1.9096</td></tr><tr><td>R<sup>(*)</sup><sub>66</sub> (30%)</td><td>1.7921</td><td>1.7722</td><td>1.9535</td><td>1.7775</td><td>1.8954</td><td>2.8353</td><td>2.0357</td><td>2.0118</td><td>2.0878</td></tr><tr><td>R<sup>(*)</sup><sub>22</sub> (10%)</td><td>3.2059</td><td>3.1673</td><td>3.2145</td><td>3.2509</td><td>6.9821</td><td>5.8175</td><td>3.4141</td><td>3.2980</td><td>3.2716</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (10%)</td><td>0.6589</td><td>0.6377</td><td>0.6726</td><td>0.6716</td><td>0.8409</td><td>0.7273</td><td>0.6706</td><td>0.6590</td><td>0.6512</td></tr><tr><td>R<sup>(*)</sup><sub>66</sub> (10%)</td><td>0.6315</td><td>0.6290</td><td>0.6547</td><td>0.6294</td><td>0.6475</td><td>0.7268</td><td>0.6706</td><td>0.6590</td><td>0.6625</td></tr><tr><td rowspan="6">600</td><td>R<sup>(*)</sup><sub>22</sub> (20%)</td><td>3.7383</td><td>3.6285</td><td>3.7659</td><td>3.8580</td><td>11.8724</td><td>9.4425</td><td>4.2871</td><td>3.9976</td><td>3.9039</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (20%)</td><td>0.8042</td><td>0.7497</td><td>0.8426</td><td>0.8310</td><td>1.6098</td><td>1.1235</td><td>0.8071</td><td>0.7880</td><td>0.7660</td></tr><tr><td>R<sup>(*)</sup><sub>66</sub> (20%)</td><td>0.7262</td><td>0.7209</td><td>0.7778</td><td>0.7219</td><td>0.7655</td><td>1.1238</td><td>0.8072</td><td>0.7880</td><td>0.7989</td></tr><tr><td>R<sup>(*)</sup><sub>22</sub> (30%)</td><td>4.4583</td><td>4.2105</td><td>4.5513</td><td>4.6886</td><td>17.4378</td><td>13.6729</td><td>5.4964</td><td>5.0246</td><td>4.7535</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (30%)</td><td>0.9906</td><td>0.8802</td><td>1.0853</td><td>1.0315</td><td>3.0636</td><td>1.8531</td><td>0.9674</td><td>0.9510</td><td>0.8968</td></tr><tr><td>R<sup>(*)</sup><sub>66</sub> (30%)</td><td>0.8358</td><td>0.8257</td><td>0.9345</td><td>0.8293</td><td>0.9009</td><td>1.8522</td><td>0.9675</td><td>0.9509</td><td>0.9934</td></tr><tr><td rowspan="6">600</td><td>R<sup>(*)</sup><sub>22</sub> (10%)</td><td>2.2114</td><td>2.1846</td><td>2.2174</td><td>2.2430</td><td>6.0825</td><td>4.8869</td><td>2.3576</td><td>2.2751</td><td>2.2568</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (10%)</td><td>0.4545</td><td>0.4398</td><td>0.4640</td><td>0.4633</td><td>0.6258</td><td>0.5230</td><td>0.4628</td><td>0.4546</td><td>0.4492</td></tr><tr><td>R<sup>(*)</sup><sub>66</sub> (10%)</td><td>0.4355</td><td>0.4338</td><td>0.4519</td><td>0.4341</td><td>0.4478</td><td>0.5226</td><td>0.4628</td><td>0.4545</td><td>0.4570</td></tr><tr><td>R<sup>(*)</sup><sub>22</sub> (20%)</td><td>2.5791</td><td>2.5031</td><td>2.5985</td><td>2.6636</td><td>10.9497</td><td>8.4549</td><td>2.9666</td><td>2.7593</td><td>2.6941</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (20%)</td><td>0.5549</td><td>0.5172</td><td>0.5815</td><td>0.5735</td><td>1.3449</td><td>0.8837</td><td>0.5572</td><td>0.5437</td><td>0.5284</td></tr><tr><td>R<sup>(*)</sup><sub>66</sub> (20%)</td><td>0.5009</td><td>0.4972</td><td>0.5372</td><td>0.4980</td><td>0.5297</td><td>0.8841</td><td>0.5572</td><td>0.5436</td><td>0.5512</td></tr><tr><td rowspan="6">600</td><td>R<sup>(*)</sup><sub>22</sub> (30%)</td><td>3.0773</td><td>2.9053</td><td>3.1427</td><td>3.2405</td><td>16.5020</td><td>12.6307</td><td>3.8148</td><td>3.4727</td><td>3.2830</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (30%)</td><td>0.6836</td><td>0.6072</td><td>0.7497</td><td>0.7122</td><td>2.7572</td><td>1.5785</td><td>0.6680</td><td>0.6563</td><td>0.6188</td></tr><tr><td>R<sup>(*)</sup><sub>66</sub> (30%)</td><td>0.5765</td><td>0.5696</td><td>0.6462</td><td>0.5721</td><td>0.6231</td><td>1.5775</td><td>0.6680</td><td>0.6562</td><td>0.6858</td></tr></tbody></table>

marized in Table 2. We take all this magnitudes with GPa dimension and pointed out this possible limitation in the paper Kouri et al. (2016).

Now, in order to support our results, we perform a numerical experiment only for the case of a fiber reinforced composite. That is, we consider the fiber's volumetric fraction equal to zero, so that the heterogeneous material is reduced to the homogeneous matrix with the properties of the soft material in Table 2. In this respect, the coefficients $R_{12}^{(*)}$ and $R_{44}^{(*)}$ for $t = 0$ can be obtained from (25a)-(25b) as,

$$G(0) = G_{\infty} + \sum_{i=1}^{N_G} G_i = 1677.979 \qquad (27a)$$

$$K(0) = K_{\infty} + \sum_{j=1}^{N_K} K_j = 3300 \qquad (27b)$$

$$R_{44}^{(*)}(0) = G(0) = 1677.979 \qquad (27c)$$

$$R_{12}^{(*)}(0) = K(0) - \frac{2}{3}G(0) = 2181.347 \qquad (27d)$$

In particular, Fig. 8(c)-(d) highlights the values of the effective coefficients $R_{12}^{(*)}$ and $R_{44}^{(*)}$ for $t = 0.01$ when the fiber's volumetric fraction is equal to zero. As predicted, they are very close to those calculated in Eqs. (27c)-(27d) for $t = 0$ and provide positive evidence in support of the present study.

### 4.5. Tabulation of results for different inclusions shapes

We conclude our discussion by presenting a tabulation of the effective relaxation modulus through three different numerical tests (see Tables 3–5). In particular, each case refers to one of the models analyzed in Sections 4.2–4.4 which describe
<table><caption>**Table 5**
The input data for the sub-phases are the same as Section 4.4.</caption><thead><tr><th rowspan="2">TIME (s)</th><th colspan="10">Prony's series</th></tr><tr><th>R<sup>(*)</sup><sub>ij</sub></th><th>LCF</th><th>LEF</th><th>LRF</th><th>LSF</th><th>Two LSF</th><th>Three LSF</th><th>CI</th><th>SI</th><th>EI</th></tr></thead><tbody><tr><td rowspan="4">0.01</td><td>R<sup>(*)</sup><sub>11</sub> (10%)</td><td>6330.9686</td><td>6463.9069</td><td>6353.2006</td><td>6395.5246</td><td>6319.6010</td><td>7922.7994</td><td>6679.2780</td><td>6529.3003</td><td>6677.3653</td></tr><tr><td>R<sup>(*)</sup><sub>12</sub> (10%)</td><td>2471.7239</td><td>2472.4374</td><td>2496.2939</td><td>2450.4499</td><td>2500.0625</td><td>2624.5756</td><td>2479.4462</td><td>2493.5456</td><td>2496.7699</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (10%)</td><td>1970.2266</td><td>1923.0107</td><td>1994.0628</td><td>1992.2400</td><td>2125.9549</td><td>2000.0913</td><td>1977.5356</td><td>1959.6305</td><td>1941.9464</td></tr><tr><td>R<sup>(*)</sup><sub>11</sub> (20%)</td><td>7354.3096</td><td>7744.2020</td><td>7421.8783</td><td>7506.3527</td><td>7241.6860</td><td>10784.4328</td><td>8167.4082</td><td>7850.5967</td><td>8376.1987</td></tr><tr><td rowspan="4">14.5</td><td>R<sup>(*)</sup><sub>12</sub> (20%)</td><td>2779.1465</td><td>2779.3829</td><td>2839.1636</td><td>2719.1023</td><td>2862.4893</td><td>3271.4743</td><td>2780.4459</td><td>2820.5936</td><td>2828.9589</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (20%)</td><td>2332.3499</td><td>2214.6194</td><td>2395.7865</td><td>2377.3201</td><td>2884.3603</td><td>2467.2618</td><td>2311.3543</td><td>2281.1879</td><td>2232.2256</td></tr><tr><td>R<sup>(*)</sup><sub>11</sub> (30%)</td><td>8675.7725</td><td>9625.4785</td><td>8873.7830</td><td>8931.4620</td><td>8348.5790</td><td>14134.7816</td><td>10038.2104</td><td>9611.7325</td><td>11428.1469</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (30%)</td><td>3117.5790</td><td>3115.7507</td><td>3234.0979</td><td>3006.9721</td><td>3312.1272</td><td>4171.26466</td><td>3126.8903</td><td>3194.9530</td><td>3258.5255</td></tr><tr><td rowspan="4">60</td><td>R<sup>(*)</sup><sub>11</sub> (10%)</td><td>3420.2613</td><td>3444.4360</td><td>3421.1834</td><td>3450.3520</td><td>3440.0138</td><td>5265.4117</td><td>3528.1570</td><td>3452.4691</td><td>3481.5114</td></tr><tr><td>R<sup>(*)</sup><sub>12</sub> (10%)</td><td>3080.4828</td><td>3082.1377</td><td>3090.7180</td><td>3060.7676</td><td>3022.6373</td><td>3132.2673</td><td>3056.9273</td><td>3083.1828</td><td>3082.7123</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (10%)</td><td>159.6480</td><td>153.5738</td><td>164.3080</td><td>163.9929</td><td>273.9138</td><td>210.6148</td><td>167.4352</td><td>164.0731</td><td>162.0413</td></tr><tr><td>R<sup>(*)</sup><sub>11</sub> (20%)</td><td>3882.0742</td><td>3967.2178</td><td>3883.6528</td><td>3974.8329</td><td>3913.4023</td><td>7949.9654</td><td>4219.6897</td><td>3993.3174</td><td>4143.3180</td></tr><tr><td rowspan="4">60</td><td>R<sup>(*)</sup><sub>12</sub> (20%)</td><td>3399.6095</td><td>3401.0045</td><td>3426.0061</td><td>3330.2289</td><td>3259.6196</td><td>3659.07117</td><td>3306.1255</td><td>3395.9505</td><td>3386.8793</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (20%)</td><td>198.8048</td><td>182.9519</td><td>212.0436</td><td>207.9328</td><td>740.7131</td><td>440.8431</td><td>205.7879</td><td>200.7157</td><td>194.7131</td></tr><tr><td>R<sup>(*)</sup><sub>11</sub> (30%)</td><td>4503.7518</td><td>4774.3710</td><td>4517.7570</td><td>4713.3316</td><td>4526.0909</td><td>11138.8076</td><td>5246.4682</td><td>4795.6902</td><td>5887.4423</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (30%)</td><td>3761.2571</td><td>3754.4058</td><td>3818.5060</td><td>3594.4532</td><td>3532.8398</td><td>4417.5313</td><td>3542.4684</td><td>3730.6670</td><td>3686.6281</td></tr><tr><td rowspan="4">60</td><td>R<sup>(*)</sup><sub>11</sub> (10%)</td><td>2203.2138</td><td>2211.4122</td><td>2205.1478</td><td>2213.7889</td><td>2211.8782</td><td>4027.5795</td><td>2237.2881</td><td>2212.8662</td><td>2221.0832</td></tr><tr><td>R<sup>(*)</sup><sub>12</sub> (10%)</td><td>2100.5041</td><td>2102.6491</td><td>2105.0875</td><td>2096.0546</td><td>2060.8204</td><td>2169.8632</td><td>2094.4675</td><td>2102.6113</td><td>2102.6042</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (10%)</td><td>49.0575</td><td>47.3497</td><td>50.3364</td><td>50.2619</td><td>150.6666</td><td>96.3096</td><td>51.0229</td><td>49.9909</td><td>49.4189</td></tr><tr><td>R<sup>(*)</sup><sub>11</sub> (20%)</td><td>2481.3152</td><td>2508.8586</td><td>2485.4144</td><td>2515.6742</td><td>2497.1469</td><td>6536.0159</td><td>2594.9182</td><td>2513.5329</td><td>2554.2208</td></tr><tr><td rowspan="4">60</td><td>R<sup>(*)</sup><sub>12</sub> (20%)</td><td>2340.9156</td><td>2345.3027</td><td>2351.5403</td><td>2321.8900</td><td>2237.5040</td><td>2634.8713</td><td>2309.0240</td><td>2342.6406</td><td>2338.7861</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (20%)</td><td>60.3727</td><td>55.9703</td><td>63.9319</td><td>62.8304</td><td>568.5508</td><td>298.0519</td><td>62.0211</td><td>60.4567</td><td>58.7335</td></tr><tr><td>R<sup>(*)</sup><sub>11</sub> (30%)</td><td>2844.1822</td><td>2932.0577</td><td>2855.1898</td><td>2932.8775</td><td>2871.3467</td><td>9527.5119</td><td>3136.2394</td><td>2927.9989</td><td>3402.6565</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (30%)</td><td>2633.4308</td><td>2639.9065</td><td>2657.7396</td><td>2574.2137</td><td>2444.0244</td><td>3326.9868</td><td>2536.3056</td><td>2624.1038</td><td>2602.0840</td></tr><tr><td>R<sup>(*)</sup><sub>44</sub> (30%)</td><td>75.0864</td><td>66.0692</td><td>83.9259</td><td>78.9781</td><td>1449.0097</td><td>725.3000</td><td>74.9251</td><td>73.6305</td><td>69.2229</td></tr></tbody></table>

the viscoelastic behavior of the constituents, i.e. four parameter model, Maxwell's model and Prony's series, respectively. Specifically, we closely examine, for fixed values of time, some of the numerical results shown in Figs. 4–8. Besides, we go further and perform the computations for all the inclusions shapes displayed in Fig. 9. Regarding the geometrical details, we consider a unit periodic cell. Moreover, we notice that the aspect ratio for the ellipsoidal inclusion is $b = c = \frac{4}{5}a$ and for the elliptical fibers $b = \frac{4}{5}a$, where $a$ is the largest semi axis. It is worth mentioning that these computations used Mesh B and the fibers are considered aligned along $Ox_3$.

## 5. Conclusions

In this paper, we use a semi-analytical technique that naturally appear as a combination of the two-scale asymptotic homogenization method and a three dimensional computational study, by means of the finite element software COMSOL Multiphysics®. Particularly, the methodology is framed to the investigation of the effective creep and relaxation behavior of non-aging linear viscoelastic composite materials. The results derived from the AHM's procedure are based on the elastic-viscoelastic correspondence principle and the Laplace-Carson transform, and comprised the local and homogenized problem, and the expression for the effective coefficients. Moreover, the stress jump conditions and the corresponding interface loads for composites with discontinuous material properties and perfect contact at the interfaces as well as significantly details of the three dimensional computational study are provided.

The strength of our work lies in achieving accurate results for the calculation of the effective properties of composite materials reinforced with several types of inclusions. This is evidenced by the good correlation and low relative errors in the comparisons with others methods reported in the literature. In this respect, a study of the solution's convergence through three different mesh discretizations for the particular case of linear elastic composites is performed. In the computational simulations, we describe the viscoelastic behavior of the constituents by having recourse to Maxwell's model, four parameter model and Prony's series.
One of the limitations of this work is related to the computational time involved to pass from the Laplace-Carson space to the time domain. One way to deal with this issue can be through parallel computing. This concern is part of our current investigations.

Finally, we take the occasion to notice that our work can be adapted to model biological and bioinspired materials with hierarchical structure, e.g. bones (Eberhardsteiner, Hellmich, & Scheiner, 2014; Nguyen, Vu, Vu, & To, 2018; Penta, Raum, Grimal, Schrof, & Gerisch, 2016). Furthermore, the present framework can also be suitable for the study of the effective behavior of piezoelectricpolymeric composites (see Jayendran & Arockiarajan, 2015; Zhang & Wu, 2010).

Future developments of the present work include the consideration of imperfect contact conditions (see Dumont, Rizzoni, Lebon, & Sacco, 2018; Ye et al., 2019) and the extension to ageing viscoelastic composite materials (see Barthlmy, Giraud, Sanahuja, & Sevostianov, 2019; Yang, Sun, Liu, Guan, & Dong, 2019).

## Declaration of Competing Interest

We confirm that the contribution entitled "On the effective behavior of viscoelastic composites in three dimensions" has not been published elsewhere and is not under consideration by another journal. All authors have furthermore approved the manuscript and agree to its publication in the International Journal of Engineering Science.

## Acknowledgements

We gratefully acknowledge to Ph.D. Qiang Chen 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é. ART kindly acknowledges the Dipartimento di Scienze Matematiche (DISMA) "G.L. Lagrange" of the Politecnico di Torino, "Dipartimento di Eccellenza 2018-2022", Project code: E11G18000350001. RP is partially supported by EPSRC grant (EP/S030875/1).

## References

* Allaire, G., & Briane, M. (1996). Multiscale convergence and reiterated homogenisation. *Proceedings of the Royal Society of Edinburgh Section A Mathematics*, 126(2), 297-342. https://doi.org/10.1017/S0308210500022757

* Auriault, J. L., Boutin, C., & Geindreau, C. (2009). Homogenization of coupled phenomena in heterogenous media. ISTE. https://doi.org/10.1002/9780470612033

* Bakhvalov, N. S., & Panasenko, G. (1989). Homogenisation: Averaging processes in periodic media: Mathematical problems in the mechanics of composite materials. Mathematics and its Applications. Kluwer Academic Publishers. https://doi.org/10.1007/978-94-009-2247-1

* Barthlmy, J.-F., Giraud, A., Sanahuja, J., & Sevostianov, I. (2019). Effective properties of ageing linear viscoelastic media with spheroidal inhomogeneities. *International Journal of Engineering Science*, 144, 103104. https://doi.org/10.1016/j.ijengsci.2019.05.015

* Bensoussan, A., Papanicolau, G., & Lions, J.-L. (1978). Asymptotic analysis for periodic structures(vol. 5). (1st ed.), North-Holland.

* Bravo-Castillero, J., Guinovart-Díaz, R., Rodríguez-Ramos, R., Sabina, F. J., & Brenner, R. (2012). Unified analytical formulae for the effective properties of periodic fibrous composites. *Materials Letters*, 73, 68-71. https://doi.org/10.1016/j.matlet.2011.12.106

* Carpenter, H. J., Gholipour, A., Ghayesh, M. H., Zander, A. C., & Psaltis, P. J. (2020). A review on the biomechanics of coronary arteries. *International Journal of Engineering Science*, 147, 103201. https://doi.org/10.1016/j.ijengsci.2019.103201

* Cavalcante, M. A. A., & Marques, S. P. C. (2014). Homogenization of periodic materials with viscoelastic phases using the generalized FVDAM theory. *Computational Materials Science*, 87, 43-53. https://doi.org/10.1016/j.commatsci.2014.01.053

* Chen, Q., Wang, G., Chen, X., & Geng, J. (2017). Finite-volume homogenization of elastic/viscoelastic periodic materials. *Composite Structures*, 182, 457-470. https://doi.org/10.1016/j.compstruct.2017.09.044

* Chen, Y., Yang, P., Zhou, Y., Guo, Z., Dong, L., & Busso, E. P. (2020). A micromechanics-based constitutive model for linear viscoelastic particle-reinforced composites. *Mechanics of Materials*, 140, 103228. https://doi.org/10.1016/j.mechmat.2019.103228

* Christensen, R. M. (1982). Theory of viscoelasticity. An introduction - 2nd edition, 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., et al. (2020). A hierarchical asymptotic homogenization approach for viscoelastic composites. *Mechanics of Advanced Materials and Structures, 0*(0), 1-12. https://doi.org/10.1080/15376494.2020.1722872

* Dastjerdi, S., Akgz, B., & Civalek (2020). On the effect of viscoelasticity on behavior of gyroscopes. *International Journal of Engineering Science*, 149, 103236. https://doi.org/10.1016/j.ijengsci.2020.103236

* Dumont, S., Rizzoni, R., Lebon, F., & Sacco, E. (2018). Soft and hard interface models for bonded elements. *Composites Part B: Engineering*, 153, 480-490. https://doi.org/10.1016/j.compositesb.2018.08.076

* Ebrahimi, F., & Barati, M. R. (2016). A nonlocal higher-order refined magneto-electro-viscoelastic beam model for dynamic analysis of smart nanostructures. *International Journal of Engineering Science*, 107, 183-196. https://doi.org/10.1016/j.ijengsci.2016.08.001

* Eberhardsteiner, L., Hellmich, C., & Scheiner, S. (2014). Layered water in crystal interfaces as source for bone viscoelasticity: Arguments from a multiscale approach. *Computer Methods in Biomechanics and Biomedical Engineering*, 17(1), 48-63. https://doi.org/10.1080/10255842.2012.670227

* Farajpour, A., Ghayesh, M. H., & Farokhi, H. (2019). Nonlocal nonlinear mechanics of imperfect carbon nanotubes. *International Journal of Engineering Science*, 142, 201-215. https://doi.org/10.1016/j.ijengsci.2019.03.003

* Gholipour, A., Ghayesh, M. H., Zander, A. C., & Psaltis, P. J. (2020). In vivo based biomechanics of right and left coronary arteries. *International Journal of Engineering Science*, 154, 103281. https://doi.org/10.1016/j.ijengsci.2020.103281

* Guinovart-Díaz, R., López-Realpozo, J. C., Rodríguez-Ramos, R., Bravo-Castillero, J., Ramírez, M., Camacho-Montes, H., et al. (2011). Influence of parallelogram cells in the axial behaviour of fibrous composite. *International Journal of Engineering Science*, 49(1), 75-84. https://doi.org/10.1016/j.ijengsci.2010.06.024

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

* Hine, P. J., & Gusev, A. A. (2019). Validating a micromechanical modelling scheme for predicting the five independent viscoelastic constants of unidirectional carbon fibre composites. *International Journal of Engineering Science*, 144, 103133. https://doi.org/10.1016/j.ijengsci.2019.103133

* Jayendran, R., & Arockiarajan, A. (2015). Numerical modelling and experimental characterization of temperature-dependent viscoelastic effect on the ferroelastic behaviour of 1-3 piezocomposites. *Sensors and Actuators A: Physical*, 226, 81-97. https://doi.org/10.1016/j.sna.2015.02.002
Juraj (2020). *Numerical inversion of laplace transforms in matlab*.

Kalamkarov, A. L., Andrianov, I. V., Pacheco, P. M. C. L., Savi, M. A., & Starushenko, G. A. (2016). Asymptotic analysis of fiber-reinforced composites of hexagonal structure. *Journal of Multiscale Modelling*, 07(03), 1650006. https://doi.org/10.1142/S1756973716500062.

Kouri, M. E., Bakkali, A., & Azrar, L. (2016). Mathematical modeling of the overall time-dependent behavior of non-ageing viscoelastic reinforced composites. *Applied Mathematical Modelling*, 40(7), 4302–4322. https://doi.org/10.1016/j.apm.2015.11.031.

Kumar, N., Khader, S. M. A., Pai, R., Khan, S. H., & Kyriacou, P. A. (2020). Fluid structure interaction study of stenosed carotid artery considering the effects of blood pressure. *International Journal of Engineering Science*, 154, 103341. https://doi.org/10.1016/j.ijengsci.2020.103341.

Lakes, R. (2009). *Viscoelastic materials*. Cambridge University Press. https://doi.org/10.1017/CBO9780511626722.

Li, Q., Chen, W., Liu, S., & Wang, J. (2019). A novel implementation of asymptotic homogenization for viscoelastic composites with periodic microstructures. *Composite Structures*, 208, 276–286. https://doi.org/10.1016/j.compstruct.2018.09.056.

Liu, X., Tang, T., Yu, W., & Pipes, R. B. (2018). Multiscale modeling of viscoelastic behaviors of textile composites. *International Journal of Engineering Science*, 130, 175–186. https://doi.org/10.1016/j.ijengsci.2018.06.003.

Maghous, S., & Creus, G. J. (2003). Periodic homogenization in thermoviscoelasticity: case of multilayered media with ageing. *International Journal of Solids and Structures*, 40(4), 851–870. https://doi.org/10.1016/S0020-7683(02)00549-8.

Malhotra, D., Pan, S., Rüther, L., Goudoulas, T., Schlippe, G., Voss, W., et al. (2019). Linear viscoelastic and microstructural properties of native male human skin and in vitro 3d reconstructed skin models. *Journal of the Mechanical Behavior of Biomedical Materials*, 90, 644–654. https://doi.org/10.1016/j.jmbbm.2018.11.013.

Mogilevskaya, S., Kushch, V., & Nikolskiy, D. (2014). Evaluation of some approximate estimates for the effective tetragonal elastic moduli of two-phase fiber-reinforced composites. *Journal of Composite Materials*, 48 (19), 2349–2362. https://doi.org/10.1177/0021998313498103.

Nguyen, S.-N., Lee, J., Han, J.-W., & Cho, M. (2020). A coupled hygrothermo-mechanical viscoelastic analysis of multilayered composite plates for long-term creep behaviors. *Composite Structures*, 242, 112030. https://doi.org/10.1016/j.compstruct.2020.112030.

Nguyen, S.-T., Vu, M.-B., Vu, M.-N., & To, Q.-D. (2018). A homogenization approach for the effective drained viscoelastic properties of 2d porous media and an application for cortical bone. *Journal of the Mechanical Behavior of Biomedical Materials*, 78, 134–142. https://doi.org/10.1016/j.jmbbm.2017.11.020.

Omari, M. A., Almableh, A., Sevostianov, I., Ashhab, M. S., & Yaseen, A. B. (2020). Modeling of the viscoelastic properties of thermoset vinyl ester nanocomposite using artificial neural network. *International Journal of Engineering Science*, 150, 103242. https://doi.org/10.1016/j.ijengsci.2020.103242.

Otero, J. A., Rodríguez-Ramos, R., Guinovart-Díaz, R., Cruz-González, O. L., Sabina, F. J., & Berger, H. (2020). Asymptotic and numerical homogenization methods applied to fibrous viscoelastic composites using Prony's series. *Acta Mechanica*. https://doi.org/10.1007/s00707-020-02671-1.

Penta, R., & Gerisch, A. (2016). Investigation of the potential of asymptotic homogenization for elastic composites via a three-dimensional computational study. *Computing and Visualization in Science*, 17. https://doi.org/10.1007/s00791-015-0257-8.

Penta, R., & Gerisch, A. (2017). The asymptotic homogenization elasticity tensor properties for composites with material discontinuities. *Continuum Mechanics and Thermodynamics*, 29(1), 187–206. https://doi.org/10.1007/s00161-016-0526-x.

Penta, R., Raum, K., Grimal, Q., Schrof, S., & Sevostianov, A. (2016). Can a continuous mineral foam explain the stiffening of aged bone tissue? A micromechanical approach to mineral fusion in musculoskeletal tissues. *Bioinspiration & Biomimetics*, 11(3). https://doi.org/10.1088/1748-3190/11/3/035004.

Persson, L., Persson, L., Svanstedt, N., & Wyller, J. (1993). *The homogenization method. An introduction*. Studentlitteratur.

Rahmanpanah, H., Mouloodi, S., Burvill, C., Gohari, S., & Davies, H. M. S. (2020). Prediction of load-displacement curve in a complex structure using artificial neural networks: A study on a long bone. *International Journal of Engineering Science*, 154, 103319. https://doi.org/10.1016/j.ijengsci.2020.103319.

Rodríguez-Ramos, R., Otero, J. A., Cruz-González, O. L., Guinovart-Díaz, R., Bravo-Castillero, J., Sabina, F. J., et al. (2020). Computation of the relaxation effective moduli for fibrous viscoelastic composites using the asymptotic homogenization method. *International Journal of Solids and Structures*, 190, 281–290. https://doi.org/10.1016/j.ijsolstr.2019.11.014.

Sanchez-Palencia, E. (1980). Non-Homogeneous media and vibration theory. *Lecture Notes in Physics*. Springer-Verlag. https://doi.org/10.1007/3-540-10000-8.

Sevostianov, I., Levin, V., & Radi, E. (2016). Effective viscoelastic properties of short-fiber reinforced composites. *International Journal of Engineering Science*, 100, 61–73. https://doi.org/10.1016/j.ijengsci.2015.10.008.

Sevostianov, I., Mogilevskaya, S. G., & Kushch, V. I. (2019). Maxwell's methodology of estimating effective properties: Alive and well. *International Journal of Engineering Science*, 140, 35–88. https://doi.org/10.1016/j.ijengsci.2019.05.001.

Smirnov, A., Vilchevskaya, E., & Sevostianov, I. (2019). Evaluation of the effective viscoelastic properties of a material containing multiple oblate inhomogeneities using fraction-exponential operators. *International Journal of Engineering Science*, 144, 103124. https://doi.org/10.1016/j.ijengsci.2019.103124.

Song, R., Muliana, A., & Rajagopal, K. (2019). A thermodynamically consistent model for viscoelastic polymers undergoing microstructural changes. *International Journal of Engineering Science*, 142, 106–124. https://doi.org/10.1016/j.ijengsci.2019.05.009.

Mainardi, F., & Spada, G. (2011). Creep, relaxation and viscosity properties for basic fractional models in rheology. *The European Physical Journal Special Topics*, 193, 133–160. https://doi.org/10.1140/epist/e2011-01387-1.

Vilchevskaya, E., Levin, V., Seyedkavoosi, S., & Sevostianov, I. (2019). Replacement relations for a viscoelastic material containing multiple inhomogeneities. *International Journal of Engineering Science*, 136, 26–37. https://doi.org/10.1016/j.ijengsci.2018.12.006.

Wang, G., & Pindera, M.-J. (2016a). Locally-exact homogenization of viscoelastic unidirectional composites. *Mechanics of Materials*, 103, 95–109. https://doi.org/10.1016/j.mechmat.2016.09.009.

Wang, G., & Pindera, M.-J. (2016b). Locally-exact homogenization theory for transversely isotropic unidirectional composites. *Mechanics Research Communications*, 78, 2–14. https://doi.org/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. *Composite Structures*, 229, 111394. https://doi.org/10.1016/j.compstruct.2019.111394.

Willoughby, N., Parnell, W. J., Hazel, A. L., & Abrahams, I. D. (2012). Homogenization methods to approximate the effective response of random fibre-reinforced composites. *International Journal of Solids and Structures*, 49(13), 1421–1433. https://doi.org/10.1016/j.ijsolstr.2012.02.010.

Yang, Z., Sun, Y., Liu, Y., Guan, T., & Dong, H. (2019). A three-scale asymptotic analysis for ageing linear viscoelastic problems of composites with multiple configurations. *Applied Mathematical Modelling*, 71, 223–242. https://doi.org/10.1016/j.apm.2019.02.021.

Ye, J., Hong, Y., Cai, H., Wang, Y., Zhai, Z., & Shi, B. (2019). A new three-dimensional parametric FVDAM for investigating the effective elastic moduli of particle-reinforced composites with interphase. *Mechanics of Advanced Materials and Structures*, 26(22), 1870–1880. https://doi.org/10.1080/15376494.2018.1452321.

Zhang, G., & Wu, M. (2010). Connectivity and shape effects on the effective properties of piezoelectricpolymeric composites. *International Journal of Engineering Science*, 48(1), 37–51. https://doi.org/10.1016/j.ijengsci.2009.06.006.