Effective Properties of Homogenised Nonlinear Viscoelastic Composites

We develop a general approach for the computation of the effective properties of nonlinear viscoelastic composites. For this purpose, we employ the asymptotic homogenisation technique to decouple the equilibrium equation into a set of local problems. The theoretical framework is then specialised to the case of a strain energy density of the Saint-Venant type, with the second Piola–Kirchhoff stress tensor also featuring a memory contribution. Within this setting, we frame our mathematical model in the case of infinitesimal displacements and employ the correspondence principle which results from the use of the Laplace transform. In doing this, we obtain the classical cell problems in asymptotic homogenisation theory for linear viscoelastic composites and look for analytical solutions of the associated anti-plane cell problems for fibre-reinforced composites. Finally, we compute the effective coefficients by specifying different types of constitutive laws for the memory terms and compare our results with available data in the scientific literature.


Introduction
Viscoelasticity is a unique property exhibited by certain materials, which encompasses a combination of viscous and elastic behaviours during deformation. Unlike purely viscous materials, viscoelastic bodies exhibit time-dependent behaviour when subjected to constant stress, deforming at a consistent rate.
The study of the mechanical properties of viscoelastic composites has been of great interest due to their unique characteristics. There is a large amount of scientific literature addressing the homogenised properties of viscoelastic composites using different homogenisation approaches [1][2][3][4][5], including those within a variational setting [6][7][8][9]. In the context of the asymptotic homogenisation technique, investigations have mainly focused on linear viscoelastic materials since one can exploit the correspondence principle introduced by the use of the Laplace transform. For instance, in [10], the Authors obtain the effective properties of linear viscoelastic composites and focus on the role of memory effects. Moreover, in [11] the two-scale asymptotic homogenisation technique together with FEM simulations are used to study the effective properties of non-ageing linear viscoelastic properties with different cell geometrical arrangements. Besides, in [12], a three-scale asymptotic homogenisation is used to model hierarchical viscoelastic materials. In [13], the Authors show the benefit of computing the viscoelastic and thermoviscoelastic properties of heterogeneous materials by means of an asymptotic homogenisation approach which treats the temperature as a parameter. On the other hand, in [14], the influence of fibre

Kinematics
Let B be a body composed by two different constituents B 1 and B 2 such that Γ is the interface between B 1 and B 2 , and B 1 ∪ B 2 = B and B 1 ∩ B 2 = B 1 ∩ B 2 = ∅, where the upper bar denotes the closure of the set. If χ : B × [t 0 , t f [→ S denotes the motion, where S is the three-dimensional Euclidean space and, t 0 and t f corresponds to the initial and final times, respectively; the current configuration of B at time t is described by B t := χ(B, t). Then, by taking the reference configurations of B 1 and B 2 to be B R1 and B R2 , respectively, we set the current configurations of B 1 and B 2 to be B 1t := χ 1 (B R1 , t) and B 2t := χ 2 (B R2 , t) where χ 1 : B R1 × [t 0 , t f [→ S and χ 2 : B R2 × [t 0 , t f [→ S represent the restrictions of the motion to B R1 and B R2 . In particular, we denote by Γ R the interface between B R1 and B R2 , and consider that B R1 ∪ B R2 = B R and B R1 ∩ B R2 = B R1 ∩ B R2 = ∅, where B R is the reference configuration of B. Then, the interface between B 1t and B 2t is Γ t := χ(Γ R , t) and is such that B 1t ∪ B 2t = B t and B 1t ∩ B 2t = B 1t ∩ B 2t = ∅. For our purposes, we consider B R2 to be characterised by the disjoint union N k=1 B k R2 . We specify by X a (a = 1, 2, 3) the coordinates of a material point in the reference configuration B R , and each spatial point with coordinates x a (a = 1, 2, 3) is given by x a = χ a (X, t). We notice that we are tacitly assuming that the coordinates x and X are referred to the same orthogonal Cartesian system. In this setting, the deformation gradient tensor, F(X, t), can be defined as where I denotes the second-order identity tensor, u is the displacement vector and Grad u is the gradient of u with components [Grad u] ab = ∂u a /∂X b .

Scales Separation
We denote the length scales characterising the composite medium and its internal structure by L c and , respectively, and introduce the smallness parameter, ε, as Adopting the considerations in [28], if Φ represents a scalar field, or a component of a vector or tensor field, we can write Φ in its two-scale version as where X = X/L c and Y = X/ =X/ε are dimensionless variables usually referred to as the slow or macroscale variable and the fast or microscale variable, respectively. Consequently, by means of the chain rule, the partial derivative of Φ with respect to X a can be expressed as In the following, for the sake of simplicity in our notation, we will simply write X and Y when referring to X and Y, respectively.

Macroscopic Uniformity and Periodicity
The scales separation induced by the macroscopic variable X and the microscopic variable Y permits selecting at the lower scale a unit cell Y , such that Y = Y 1 ∪ Y 2 , where, in the context of this work, we consider that Y 1 surrounds the inclusion Y 2 with Γ Y being the interface between Y 1 and Y 2 . We notice that, in particular, ∂Y 1 = ∂Y ∪ Γ Y and ∂Y 2 = Γ Y , where ∂Y is the external boundary of the cell. In this work, we adopt the assumption of macroscopic uniformity (see, for instance [29][30][31] for further discussions). That is, we consider that the topological properties of the microscopic cell are independent of the macroscopic variable ( Figure 1A). There are two main consequences of the consideration of macroscopic uniformity. The first one is that we can select the unit cell independently of the position at the macroscopic domain and, thus, it is representative of the whole microstructure. The second consequence is that it allows to interchange the operator of differentiation with respect to the macroscopic variable with the integration over the repeating unit cell, namely, While the above considerations are of great help when dealing with calculations, the assumption of macroscopic uniformity can be relaxed by paying the price of elevated computational cost (see, for example, [32,33]).

Multiscale Formulation of the Problem
We set our problem within a purely mechanical framework by considering the equilibrium equation in the absence of body forces for each constituent, which, in its material form, can be written as and P η denotes the first Piola-Kirchhoff stress tensor for each constituent, defined by η = 1, 2. Furthermore, Equation (6) is complemented with conditions of ideal contact at the interface Γ between B R1 and B R2 , namely, where u η represents the displacement vector and N is the normal vector pointing from B R2 to B R1 . In Equations (7a) and (7b) the notations Φ 1 (X − , t) and Φ 2 (X + , t) are defined as where X ± ∈ Γ ± R with Γ − R being the part of the interface Γ R in contact with B R1 and Γ + R the part of Γ R in contact with B R2 .
Recalling the scale separation introduced in Section 2.2 and the assumption of macroscopic uniformity, the problem specified by Equations (6), (7a) and (7b) takes the form The first Piola-Kirchhoff stress tensor, P ε η , is related to the second Piola-Kirchhoff stress tensor, S ε η , through the formula which, for nonlinear viscoelastic constituents, can be written in the form [35] where ψ ε η denotes the strain energy density for each constituent and E ε η := 1 2 (C ε η − I) is the Green-Lagrange strain tensor, with C ε η := (F ε η ) T F ε η being the right Cauchy-Green deformation tensor. Furthermore, by embracing the terminology in [35], L η is the fourth order tensor referred to as the relaxation tensor which we assume to have left and right minor symmetries and major symmetry.

Asymptotic Homogenisation Procedure
By taking inspiration from [22,24], we propose a formal series representation for the the displacement u ε η and the first Piola-Kirchhoff stress tensor P ε η in powers of the smallness parameter ε, namely, The substitution of (11a) in the expression of the deformation gradient tensor (1) leads to the following power expansion where Consequently, the series expansion of the Green-Lagrange strain tensor E ε η in powers of ε is given by where, for each k = 1, 2, . . ., and Thus, the combination and substitution of the above results in the multiscale problem defined by Equations (8a)-(8c), together with the fact that the first and second Piola-Kirchhoff stress tensors are related by the formula P ε η = F ε η S ε η , leads to a set of cell problems obtained by equating the coefficients of the powers of ε. Precisely, at the leading order, i.e., ε 0 , we have that with Furthermore, for increasing powers of ε, namely, ε k with k = 1, 2, 3, . . ., we can write where In particular, for k = 1, we have that

Constitutive Considerations
To exemplify the results obtained so far, we assume that the strain energy density ψ ε η is of Saint-Venant type. That is, by denoting with C ε η the elasticity tensor for each constituent B η , with η = 1, 2, we have that where C ε η is endowed with both the left and right minor symmetries, and with the major symmetry, namely with a, b, c, d = 1, 2, 3. The introduction of the constitutive law (23) allows to compute the derivative of ψ ε η with respect to E ε η in the form where the terms E (k) η are given in Equations (15a) and (15b). Therefore, we can individuate the coefficients in the expansion of the second Piola-Kirchhoff stress tensor as follows which under the hypothesis of causal histories, that is by considering E , with H(t) denoting the Heaviside step function, we have that where G ε η is the fourth order tensor defined as G ε η := C ε η + L ε η . In the following, we simplify the notation by writing E

The Cell and Homogenised Problems for Infinitesimal Displacements
In order to take advantage of the elastic-viscoelastic correspondence principle (see, e.g., [12,15]), we frame the results of the previous sections in the setting of infinitesimal displacements. To this end, we introduce the infinitesimal displacement gradient H η (X, t) = Grad u η (X, t), where, with abuse of notation, u η denotes the infinitesimal displacement. In this context, the leading order terms of the deformation gradient tensor are given by where H for k = 0, 1. Thus, substituting the above results in the expansion of the right Cauchy-Green deformation tensor (16) and retaining only the linear terms, we can write Therefore, the expressions (15a) and (15b) for the Green-Lagrange strain tensors are, within this linearised setting, given by with k = 0, 1. Hence, by employing the right minor symmetry of G ε η , we can deduce that Thus, after substitution in (18a) and (21a), and retaining only the linear terms, we can deduce that

The First Cell Problem
We notice that the integral in (30a) represents the convolution of G ε η andḢ η , so that the application of the Laplace transform leads tõ where s is a complex variable dual to time andĜ ε η denotes the Laplace-Carson transform of G ε η , that is,Ĝ ε η = sG ε η = C ε η + sL ε η withG ε η andL ε η being the Laplace transform of G η and L ε η . Since the elasticity tensor, C ε η , is time independent, we can also writeĜ ε η =Ĉ ε η +L ε η . Therefore, by introducing the ansatz whereχ η represents the Laplace transform of the Y-periodic, third order tensor χ η andω η is a Y-constant vector field, we find, from Equation (17a), the cell or local problem We notice that Equations (33a)-(33c) determine the classical cell problem resulting from the use of the asymptotic homogenisation in linear viscoelasticity (see, for instance, [10,11,15]). Following the discussion in [24], in Equation (33a), the components of

The Second Cell Problem
The second cell problem is given by Equations (19a)-(19c) for k = 1. In view of the results obtained in the previous sections, the first equation of the second cell problem, that is Equation (19a), can be equivalently rewritten as Thus, by introducing, for a generic field Φ ε η (X, Y, s), the integral operators • Y η over the portion Y η of the unit cell Y , namely where |Y | represents the measure of the cell, and applying it to Equation (35), we can deduce that where we have taken into account the consideration of macroscopic uniformity. Therefore, since the second addend at the left-hand side of (37) reduces to the zero vector because of the assumption of local periodicity and the ideal contact conditions at the interface (refer to [28] for further details), we find the homogenised equation where the effective coefficientĜ eff is defined aŝ

Uniaxially Fibre-Reinforced Composites
We consider the case of a fibre reinforced composite with uniaxially oriented cylindrical fibres oriented in the direction specified by i 3 , where {i k } 3 k=1 denotes the standard Cartesian vector basis. In this framework, we represent the cross section of the microstructure as a matrix with evenly distributed circular inclusions and, especially, we choose the unit cell to be characterised by a square with a single inclusion Y 2 . That is, the portion of the cell representing the matrix is Y 1 = Y \ Y 2 as shown in Figure 2. Within this context, the unknowns and fields of interest depend only on the macroscale spatial coordinates X 1 and X 2 , on the microscale coordinates Y 1 and Y 2 , and on s. Consequently, we have that the leading order-term of the infinitesimal displacement gradient tensor, H (0) η , can be additively decomposed as follows where the decomposition (40), together with the linearity of the problem at hand and the superposition principle, permit us to split the cell problem (33a) into an in-plane and an anti-plane contribution. Here, we focus on solving only the anti-plane cell problem, which means that we work only with the second addend in (40), so that the anti-plane shear stresses are where the subindex "a" inP (0) η a is used to specify that only the anti-plane contribution is being taken into account.

Anti-Plane Cell Problems and Effective Coefficients
For generic monoclinic constituents, the relevant constitutive relation is [36,37], with b = 1, 2. Thus, because of the above considerations, the relevant components of the effective coefficient associated with the anti-plane problem are given by with b, d = 1, 2, which can be found by solving the cell problems, Equations (44a)-(44c), for d = 1, 2, represent two anti-plane cell problems and each one of them help to compute specific effective coefficients (see, e.g., [16] [15], and that the fourth-order tensorsĜ ε η are of the formĜ ε η (X, Y, s) =Ĝ ε η (s) (i.e., we assume that the material properties of each constituent depend only on s), the expressions for the effective coefficients corresponding to the anti-plane problem can be written as and the cell problem reduces to

Solution of the Anti-Plane Cell Problem
To find the solution of the the cell problem (46a)-(46c), we adhere to the procedure discussed in [15] for a square arrangement of cells with an embedded circular inclusion. The methodology, which is based on the the theory of harmonic functions and the Kolosov-Muskhelishvili complex potentials [38] has been extensively investigated in different scenarios (see, for instance, [16,39]). Here, we just report the main steps and refer the Reader to [15,16]. Therefore, we set where, for d = 1, 2, the complex potentials can be written as In particular, the symbol "o" indicates that the sum is performed over odd indices, a d l (s) and c d l (s) are complex coefficients and A d l (s) := ∑ ∞ m=1 o mΛ ml a d m (s). In the latter expression, the elements of the matrix Λ are defined as where the term S m+l denotes the reticular sums, namely S m+l = ∑ w∈L * w −(m+l) for w = rw 1 + sw 2 with r, s ∈ Z and L * representing the lattice excluding w = 0 (see [17] for further details). The substitution of the above expressions in the interface condition (46b) leads to the following algebraic equations involving the complex coefficients a d l and c d l where d = 1, 2, R is the radius of the circular inclusion and Z = Re iθ , with θ ∈ [0, 2π], characterises the interface between the inclusion and the matrix. Furthermore, following similar procedure on the second interface condition (46c) yields where ξ(s) is given by the expression In particular, rescaling the problem (51) with a d l (s) = [b d l (s)R l ]/ √ l and considering the real and imaginary parts, the system of linear Equation (51) can be equivalently rewritten in the form where (·) and (·) extract the real and imaginary parts of the complex quantities they are applied to. We notice that, Equations (53a) and (53b) can be represented in the equiva- where we have introduced the notations ξ r and ξ i for the real and imaginary parts of ξ, and Furthermore, we denote with I and O the infinite identity and zero matrices, respectively. We remark that in the limit case of linear elastic constituents, Equation (54) reduces to the one provided in [24].

Effective Coefficients
To find closed-form expressions for the effective coefficients, we follow, with slight modifications, the approach depicted in [16,17]. We notice that, for the problem at hand, the effective coefficients given in (43) can be specified as so that, due to the local periodicity ofχ η , we can write and Thus, after substitution of (47) in the above expressions, it is possible to write In the case of tetragonal symmetry of the fourth-order tensorĜ ε η , we know that [Ĝ eff ] 3132 = [Ĝ eff ] 3231 = 0 (see, for instance, [15]). Consequently, (59a) and (59b) reduce to

Benchmark Problems
Let us consider that both constituents B 1 and B 2 are isotropic. In the present framework, this consideration implies that the elastic and relaxation tensors have the following representations where κ e η and µ e η represent the elastic bulk modulus and the second Lamé's parameter, respectively, and κ v η and µ v η are memory functions which need to be constitutively chosen. Furthermore, the fourth order tensor K := 1 3 I ⊗ I extracts the spherical part of a secondorder tensor, while M := 1 2 [I⊗I + I⊗I] − K extracts its deviatoric part. Thus, recalling thatĜ ε η =Ĉ ε η +L ε η , we can deduce that To compute the effective coefficients, we truncate the system of Equation (54) to different orders of m = l and find a d 1 (s) in (62a) and (62b). We notice that, to obtain the effective coefficients as functions of time, we need to apply the inverse Laplace-Carson transform. For this purpose, we adopt part of the scheme used in, e.g., [12] which employs the MATLAB function INVLAP [40].

Fibre Reinforced Composite with Elastic Inclusion and Kelvin-Voight Matrix
As a first specialisation of the developed theory, we consider the case of a composite with viscoelastic matrix of Kelvin-Voigt type and elastic inclusion. Specifically, for the matrix constituent, we assume that the memory function µ v 1 is given by the expression In Figure 3, we report the numerical values of the effective coefficient [G eff ] 3131 with respect to time (that is, after computing the inverse Laplace-Carson transform) for different orders of truncation of the linear system. In particular, we choose the set of parameters µ e 1 = 0.75 GPa, µ e 2 = 0.375 GPa and ν 1 = 0.375 GPa·s [41]. In particular, in Figure 3, we notice that for the truncation orders l = m = 1 and l = m = 3, the effective coefficients are almost indistinguishable, which suggests that the numerical computations converge rapidly, making them cost-efficient. Because of the good agreement shown in Figure 3, we continue our analysis with the truncation order l = m = 1, unless indicated otherwise. We also advert that as the fibre volumetric fraction decreases, the effective coefficient [G eff ] 3131 tends to the shear modulus of the viscoelastic matrix, that is µ e 1 = 0.75 GPa. This behaviour agrees with the fact that, physically, the material becomes a viscoelastic homogeneous medium.

Fibre Reinforced Composite Characterised by Prony Series
In this second benchmark problem, we consider that both constituents are characterised by a Prony series. That is, we set where µ e ηi denote, for each i = 1, 2, . . . , N, the shear elastic modulus and τ ηi are the corresponding relaxation times. Consequently, the effective coefficients given in (62a) and (62b) can be written in the form To perform the numerical simulations, we consider the same parameters specified in [11], which for convenience are reported in Table 1. Specifically, the matrix is chosen as a soft viscoelastic material and the fibre as a stiff viscoelastic material. The results of our simulations are shown in Figure 5 and we compare them with those obtained in [11], in which the Authors used numerical simulations based on finite elements. Specifically, we select the numerical results reported in [11] for V 2 = 0.1 and V 2 = 0.3. These are represented as red dots in Figure 5 and agree with the values computed in this work. Additionally, in Figure 6, we also report the time-dependent behaviour of the effective coefficient [G eff ] 3131 for different values of the fibre volumetric fraction. While the differences in variations of the fibre diameter may not be as significant as in the previous benchmark problem, we still observe the influence of V 2 . Predictably, the loss of the stiffness provided by the fibre produces a softening of the composite.

Conclusions
In this work, we proposed a general formulation for the study of nonlinear viscoelastic composites using the asymptotic homogenisation technique in a purely mechanical framework involving the equilibrium equation in the absence of body forces. After obtaining the local problems associated with the original formulation, we showed the potential of the scheme by considering the Saint-Venant strain energy density and formulating the cell and homogenised problems in the limit of infinitesimal deformations. This choice was mainly due to the numerical complications arising in a finite theory and the convenience of the use of the correspondence principle in linear problems. Furthermore, we analytically solved the anti-plane local problems associated with uniaxially fibre-reinforced materials and our results showed good agreement with published data. We remark that the model developed here is versatile and that, although only the effective coefficients for Kelvin-Voigt and Prony series models were reported, it can be easily adapted to different selections of viscoelastic laws.
As we mentioned before, the selection of the Saint-Venant strain energy density was mostly driven by its simplicity. Nevertheless, the methodology presented here could be adapted to other types of strain energy densities. However, this will limit the use of the correspondence principle and, consequently, appropriate numerical schemes need to be considered for the time integral (see, e.g., [42]). We further remark that the method to find the solution to the cell problems can be extended, for example, as done in [15], to other topological settings.
Further advancements in this work encompass a wider analysis of diverse types of nonlinear viscoelastic materials and contextualising the general theory in relevant biological scenarios. For instance, in the case of biological fibrous tissues, such as muscles and connective tissues. Further research in this area is expected to lead to new research questions in materials science and biomathematics.