A hierarchical asymptotic homogenization approach for viscoelastic composites

Abstract Effective properties of non-aging linear viscoelastic and hierarchical composites are investigated via a three-scale asymptotic homogenization method. In this approach, we consider the assumption of a generalized periodicity in the different structural levels and their characterization through the so-called stratified functions. The expressions for the associated local and homogenized problems, and the effective coefficients are derived at each level of organization by using the correspondence principle and the Laplace-Carson transform. Considering isotropic components and a perfect contact at the interfaces between the constituents, analytical solutions, in the Laplace-Carson space, are found for the local problems and the effective coefficients are computed. An interconversion procedure between the effective relaxation modulus and the effective creep compliance is carried out for obtaining information about both viscoelastic properties. The numerical inversion to the original temporal space is also performed. Finally, we exploit the potential of the approach and study the overall properties of a hierarchical viscoelastic composite structure representing the dermis.


Introduction
In the scientific literature there exist several works focusing on the development of micromechanical techniques to predict the macroscopic properties of composite materials. An excellent review on these methods can be found in Sevostianov and Giraud [1]; Davit et al. [2]. In particular, multiscale asymptotic homogenization methods take advantage of the information available at the smaller scales of a given heterogeneous medium to predict the effective properties at its larger scales (see e.g. Lukkassen and Milton [3]; Penta and Gerisch [4]; Ram ırez-Torres et al. [5]; Ram ırez-Torres et al. [6]; Yang et al. [7]; Dong et al. [8]). This homogenization procedure requires the solution of cell problems with input data corresponding to the homogenized material properties resulting in previous steps. In addition, the use of more general periodic or stratification functions (see e.g. Tsalis, Chatzigeorgiou, and Charalambakis [9]; Tsalis et al. [10]) permits to describe the different length scales of the composite materials. They can be applied to homogenization problems related to hierarchical real-world composites as human tissue (see e.g. Penta et al. [11] for bone and tendons) that generally form geometries which cannot be produced by the repetition of one elementary volume.
In the present work, the three-scale Asymptotic Homogenization Method (AHM) introduced in Ram ırez-Torres et al. [5]; Ram ırez-Torres et al. [6]; Ram ırez-Torres et al. [12] is used to model a non-aging linear viscoelastic composite material with generalized periodicity and two hierarchical levels. This three-scale asymptotic approach has been used for modeling hierarchical laminated and fiber-reinforced elastic composites in Ram ırez-Torres et al. [5] and Ram ırez-Torres et al. [6]; Ram ırez-Torres et al. [12], respectively. Also, there are several other works dealing with the modeling of multiscale hierarchical heterogeneous media. For example, from the theoretical point of view we find the works Allaire and Briane [13]; Bensoussan, Papanicolau, and Lions [14]; Trucu, Chaplain, and Marciniak-Czochra [15]. Regarding these results, a rigorous study about the multiscale convergence is developed in Bensoussan, Papanicolau, and Lions [14]; here the authors consider that y ¼ x=e and z ¼ x=e 2 are the local variables. Subsequently, and applied to the heat equation for composites, the method of reiterate homogenization is well founded in Allaire and Briane [13]. An additional overview of reiterated homogenization is presented in Trucu, Chaplain, and Marciniak-Czochra [15] via a three-scale convergence approach where the asymptotic parameters independently approach zero. More recently, in Dong et al. [8]; Yang et al. [7,16], the authors study the properties of thermo-mechanical, non-aging and aging viscoelastic composites with multiple spatial scales by using a three-scale asymptotic expansion and a periodic layout of the heterogeneities in the structures. They also provide a finite element algorithm based on inverse Laplace transform and the three-scale asymptotic homogenization to obtain the numerical results.
The article presents some novelties with respect to previous works of ours [5,17] and, to the best of our knowledge, with others in the literature [7,16,[18][19][20]). Specifically, we generalize the results obtained in [5] for linear elasticity by extending them to a non-aging, linear viscoelasticity framework; we deal with the stratified functions in the homogenization procedure; an analytical solution for the local problems associated with each hierarchical level is obtained; the expressions for the effective coefficients for hierarchical laminated composites with generalized periodicity, isotropic components and perfect contact at the interfaces are provided; and the methodology is used to model the overall behavior of the dermis in the skin. A better understanding of this tissue has a real impact on biomedical applications and in the development of modern technology such as flexible instance electronics, soft robotics and prosthetics (see Ref. [21]). Many researches are focused on the study of the viscoelastic, non-linear hyperelastic and anisotropic behavior of the structural components of the skin (see Ref. [22,23]) and the relation with its hierarchical structure (see Ref. [21,24]). An inherent characteristic of the mechanical properties of the skin lies in the fact that its properties vary according to the skin type, age, gender, location, and other environmental factors Panchal et al. [25]. In particular, skin has three main layers: the epidermis, the dermis and the hypodermis, in addition to the stratum corneum that acts as the outermost barrier of human skin Muha et al. [26]. Our research is focused on the dermis which represents the most important mechanical and thermal unit of the skin, as it occupies around 90% of the skin's thickness (see Ref. [27]). In the present work, the dermis is assumed to behave as a non-aging linear viscoelastic and hierarchical composite material, and thus, the investigation of its effective properties is based on the correspondence principle and the Laplace transform (see Ref. [28]). The procedure, see for instance To et al. [29] and Yang et al. [16], consists in the change of the convolution constitutive law by a fictitious elastic one in the Laplace domain. Then, the inversion of the Laplace transform is considered to derive the effective behavior in the time domain.
The article is organized as follows. The mathematical framework of the problem and the homogenization steps are illustrated in Sections 2 and 3, respectively. In Section 4, we present the analytical solution of the local problems and the calculation of the effective coefficients. In Section 5, we highlight the potential of the current approach in the investigation of the effective relaxation and creep compliance properties of the dermis.

Separation of scales
Let us consider the dimensionless, scaling parameters e 1 and e 2 , defined as follows, (see Ref. [12]) In addition, the well-separated characteristic length scales ' 1 , ' 2 and L introduce three dimensionless spatial variables in the reference configuration, where x is said to be the physical spatial variable, whereas x, y and z represent the macroscopic, mesoscopic and the microscopic dimensionless spatial variables, respectively. By using Eq. (1), x, y and z can be related through the expressions, Now, by defining a field / over the region of interest, the separation of scales allows to rephrase the space dependence of / as /ðxÞ ¼ /ð xðxÞ, yðxÞ, zðxÞÞ, and taking into account Eqs. (1)-(4), the spatial derivative of / can be computed as follows, In the present work, we consider the generalized periodicity (see Ref. [9,10]) as an additional condition in our model, and thus, we rewrite the three dimensionless spatial variables in Eq. (2) as follows where q ðaÞ : R 3 ! R 3 , with a ¼ y, z and the property q ðaÞ 2 C 1 ðXÞ, represent the so-called stratified functions. It is worthy to note that the parametric equations q ðaÞ ðxÞ ¼ constant, describe the structural organization of the constituents at the e 1 and e 2 hierarchical levels, respectively. Therefore, we have that In addition, the spatial derivative of / is given by the formula By following this approach, all equations should be written in dimensionless form. In the literature, the switch to auxiliary variables x, y and z is often omitted. However, as shown for example in Auriault, Boutin, and Geindreau [30], both paths are equivalent. Therefore, in the article, the analysis is carried out directly in a system of physical variables x, y and z, where y ¼ q ðyÞ ðxÞ e 1 and z ¼ q ðzÞ ðxÞ e 2 : Similar ideas are developed in Ram ırez-Torres et al. [31] to the case of two scales.

The physical model
Let us denote by X 2 R 3 a linear multiscale viscoelastic composite material (see Figure 1) characterized by a generalized periodicity at the different levels of organization. First stage: e 1 -structural level: The domain X is considered to be a two-phase quasi-periodic composite such that where N 1 denotes the number of inclusions, and we denote with C e1 the interface between X e1 1 and X e1 2 : In this description, Y represents the unitary periodic cell and the its constituents fulfill the constraint Second stage: e 2 -structural level At this structural level we consider that, at the same time, each phase i X e1 1 (i ¼ 1, :::, N 1 ) is a two-phase quasi-periodic composite material, where i X Figure 1(c)). Besides, X e2 1 ¼ [ N2 j¼1j X e2 1 and C e2 represents the interface between X e2 1 and X e2 2 : Besides, Z is the corresponding unitary periodic cell, in which we enforce that the constituents satisfy the following relation At this point, we extend the concept of generalized or quasiperiodicity, given in Tsalis, Chatzigeorgiou, and Charalambakis [9]; Tsalis et al. [10], to the case of hierarchical composites with three scales. The definition is presented as follows, Definition 2.1. (Generalized periodicity) A composite material with two structural levels of organization is said to have generalized periodicity (or quasi-periodicity) if there exists a global spatial coordinate system x and functions q ðaÞ : R 3 ! R 3 , with a ¼ y, z (see Eq. (10)) such that the operator F which relates the second-order stress tensor r and the second-order strain tensor n, by the relation r ¼ Fðn; x, y, zÞ, is regular in x, and periodic in y and z.
This concept will permit us to study different structural effects in the composite.
Remark 1. From now on, we consider the following notation / e ðx, tÞ ¼ /ðx, y, z, tÞ, where / e is assumed to be periodic in y and z.

Formulation of the problem
Considering that the constitutive response of all the constituents of the composite body is linear viscoelastic, and ignoring inertial terms, the problem in X reads, where, by taking into account Definition 1, r is given by the following scale-dependent constitutive law (see Ref. [32]), R e ðx, t À sÞ : @nðu e ðx, sÞÞ @s ds: In Eq. (12), n denotes the second-order strain tensor determined by the formula n u e ðx, tÞ and u is the displacement field. Moreover, f denotes the action of external volume forces that satisfies f ðx, tÞ 2 L 2 ðX Â RÞ, u and S are the prescribed displacement and traction on the boundary @X ¼ @X d [ @X n , with @X d \ @X n ¼ ; and n is the outward unit vector normal to the surface @X: The relaxation modulus is denoted by R e with the following symmetry properties R e ijkl ¼ R e jikl ¼ R e ijlk ¼ R e klij : Furthermore, we assume that R e 2 L 1 ðX Â RÞ and that it is positively defined, i.e. the relation R e ijkl g ij g kl ! bg ij g kl is verified for all symmetric real valued second-order tensors g and some positive constant b.
We further impose continuity conditions for displacements and tractions on both C e 1 and C e 2 , i.e.
where the outward unit vectors to the surfaces C e 1 and C e 2 are represented by n ðyÞ and n ðzÞ , respectively. The operator / e ½ ½ denotes the "jump" of / e across the interface between the constituents.

Re-formulation of the problem in Laplace-Carson space
The scale-dependent constitutive law Eq. (12) corresponds to the special form of a non-aging linear viscoelastic materials where we emphasize the specific functional dependency on t À s (see Ref. [33]). Therefore, the viscoelastic problem given by Eqs. (11a-14) can be transformed into an elastic one using the Laplace-Carson transform. It is defined for all real numbers t ! 0 as (see Appendix B in Christensen [32]) / e ðx, pÞ ¼ p where the variable p denotes the Laplace-Carson space.
The aforementioned transformation is known as the correspondence principle and permits us to rewrite the system Eqs. (11a-11d) in the Laplace-Carson space, Furthermore, the interface contact conditions become, vu e ðx, pÞb ¼ 0, vR e ðx, pÞ : n u e ðx, pÞ Á n ðaÞ b ¼ 0 ða ¼ y, zÞ: In what follows the functions depending on the variable p (e.g. / e ðx, pÞ) are defined in the Laplace-Carson space and the homogenization process is developed in Laplace-Carson space.
3. The three-scale asymptotic homogenization approach In this section we recall the three-scale asymptotic homogenization approach introduced in Ram ırez-Torres et al. [5] and we adapt it to the investigation of the overall behavior of hierarchical, linear viscoelastic composites with generalized periodicity of its internal structure. We first notice that, since each field and material property / e ðx, pÞ are assumed to be regular in x and periodic in y and z, the separation of scales together with Eq. (10) imply that, and using a similar idea, the Eq. (13)

Three-scale homogenization procedure
Following Ram ırez-Torres et al. [5], the solution to the problem Eqs. (15a-16) is given as follows, whereũ ð0Þ is defined as u ð0Þ ðx, y, z, pÞ ¼ u ð0Þ ðx, y, z, pÞ þ X þ1 i¼1 u ðiÞ ðx, y, z, pÞe i 1 : (21) Here and subsequently (unless necessary), the variable dependence is dropped out for convenience. Replacing expansion Eq. (20) into Eq. (15a) and using the relations Eqs. (17) and (18), it is obtained and analogously, the interface conditions Eq. (16) read, i þ :: Before proceeding, the following cell average operators are introduced, where jZj and jYj represents the volume of the periodic cell Z and Y, respectively.
where the third order tensorṽ klm is periodic in y and z and U Then, substituting Eq. (29) into Eq. (28a-28c), and after some simplifications, the following auxiliary problem, which we call the e 2 -cell problem, is obtained The expression n ðzÞ pq ðṽ kl Þ in Eq. (31a) is defined taking into account Eq. (19b) as follows Averaging Eq. (33) and taking into account the divergence theorem and the periodicity of the involved functions in the variable z (see Ref. [9,34]), we obtain that where is the effective coefficient at the e 1 -structural level.

Second level of the hierarchical structure
In the previous Section, we have carried out an analysis for finding the homogenized properties at the e 1 -structural level. Now, in order to find the overall behavior of the hierarchical composite, we follow a similar procedure where the results found so far become the new input values for the new problems. Using relations Eqs. (21) and (30) Considering the representation Eq.
Working with Eqs. (36)(37)(38) and grouping in powers of e 1 , the following sequence of problems are analyzed.
Problem for e À2 1 @q ðyÞ n @x j @ @y n ½ R ijkl n ðyÞ kl u ð0Þ ¼ 0 in ðY n C y Þ Â 0, þ 1Þ, ½ where the following result is reached (see Ref. [9]) u ð0Þ ¼ u ð0Þ ðx, pÞ: From Eqs. (37)(38) and applying the cell average operator over Z, Problem for e À1 1 At this point, using Eqs. (36)(37)(38), the fact that n ðyÞ kl ðu ð0Þ Þ ¼ 0 and applying the cell average operator over the unit cell Z, we obtain that @q ðyÞ n @x j @ @y n R ijkl n Similarly to what was said in the previous Section, the existence and uniqueness of a solution for the problem given by  -42c), the third order tensor v klm , periodic in y, is the solution of the following problem referred to as the e 1 -cell problem, Problem for e 0 1 , Equating in powers of e 0 1 , Averaging Eq. (46) over the unit cell Y, the homogenized problem is obtained and it can be written in the form is the expression for the effective relaxation modulus of the hierarchical composite media. Furthermore, from Eqs. (15b-15d), the boundary conditions associated to the problem Eq. (47) arê and the initial condition is given by Remark 3. In the previous sections and as a consequence of the application of the three-scale Asymptotic Homogenization Method, we were able to obtain the homogenized problem Eqs.
(47-51) that describes the overall behavior of a hierarchical, non-aging, linear viscoelastic composite material with generalized periodicity; the respective cell problems Eqs. (31a-31c) and Eqs. (44a-44c) and the effective coefficients Eqs. (35) and (48) for each level of organization. The procedure can be used in order to solve a linear elastic problem (t ¼ 0) in composites with the same structural characteristics. Also, we can recover the classical results of two-scale as particular cases of the formulation (see e.g. Cruz-Gonz alez et al. [17]).

Effective coefficients for hierarchical laminated composites with generalized periodicity
Laminated composites materials can be described by using a stratified function such that q : R n ! R with n ¼ 2, 3 (see Ref. [9]). Here, we consider the general case when the stratified functions are given by q ðaÞ : R 3 ! R i.e. q ðaÞ q ðaÞ ðx 1 , x 2 , x 3 Þ with a ¼ y, z: Taking into account Voigt's notation, the symmetry properties of the relaxation modulus and expressions Eqs. (32) and (45), the effective coefficients Eqs. (35) and (48) can be written, respectively, as follows @y y , 0 @ (52b) for i, j ¼ 1, 2, :::, 6: Additionally, the local problem Eq. (31a) can be transformed into the following partial differential equations system by using the Voigt's notation where i, k ¼ 1, 2, 3, j ¼ 1, 2, :::, 6 and the expressions for the components arẽ Analogously, the local problem Eq. (44a) is expressed as follows, where Systems of Eqs. (53-54d) and (55-56d) can be solved analytically by integrating each equation with respect to the local variables and by determining the constants of integration using the corresponding cell average operator (see Ref. [35] for more details). The resolution scheme to obtain the effective viscoelastic properties of a hierarchical composite can be summarized as following: (i) To solve the system Eqs. (53-54d), and then the expressions for @ṽ ð1Þ bi =@z, with i ¼ 1, 2, 3, can be substituted into Eq. (52a).
(ii) At this point, it is possible to solve the system given in Eq. (55-56d) and to obtain the expressions for @v ð1Þ bj =@y, with j ¼ 1, 2, 3. (iii) Finally, the effective viscoelastic properties can be obtained by using formula Eq. (52b).
It is worthy to remark that the expressions to the local problems Eqs. (53-54d) and Eqs. (55-56d) are developed for anisotropic composites.

An approach for modeling the mechanical properties of the dermis
The dermis is considered as a multilayer collagen-reinforced structure (see Ref. [27]). Two of its principal components are collagen and ground substance (also known as supporting matrix). On the one hand, the collagen contributes to 75% of the fat-free dry mass, 18%-30% of the volume of the dermis and can be considered to behave as a linear elastic material (see Ref. [36,37]). On the other hand, the ground substance is responsible for the viscoelastic behavior of the dermis and comprises about 20% of the dry weight of skin and makes up to between 70% and 90% of the skin's volume. It is a gel-like substance containing a class of chemicals including glycosaminoglycans (GAG), proteoglycans and glycoproteins (see Ref. [37]). In this section, an approach for modeling the dermis as a hierarchical viscoelastic composite material is proposed. We follow the structural ideas depicted in Figure 3 of Sherman et al. [21]. It is stated that the stiffness of collagenous materials is a consequence of the properties, arrangement, and geometric distribution of collagen fibrils, which are embedded in the ground substance (see Ref. [36]) and forming the dermis. These structures are long strands with wavy effects and are arranged in wavy parallel fibers, which are flattened in the plane of the dermis and determine its properties.
In this work, we present the problematic from the point of view of laminated composite materials (as represented in Figure 2). We notice that, in doing this, we take inspiration in the work of Sherman et al. [21]. In fact, Figure 2(ii) corresponds to Figure 3(b) of Sherman et al. [21] and Figure  2(i) corresponds to Figure 3(d) of Sherman et al. [21], respectively.
Moreover, we consider the dermis as a two-phase composite made of collagen (elastic constituent) and ground substance (viscoelastic constituent) organized in a hierarchical form. The mechanical properties of the ground substance can be taken approximately from McBride et al. [38]; Panchal et al. [25]; Xu and Lu [37] and the collagen's mechanical properties are obtained from Ben ıtez and Mont ans [27]; Xu and Lu [37]. They are listed in Table 1.
The viscoelastic constituent (ground substance) is modeled using the normalized Prony series (see Ref. [36]), where G(t) is the time dependent relaxation shear modulus, s n denote the relaxation times, G n represent the modulus coefficients and G 1 is the long-term modulus (see Table 2). For the sake of simplicity in the model, only two term in the Prony series are considered (see Ref. [36]). The notation V i with i ¼ 1, 2 represents the volumetric fractions of each constituent (see Fig. 2). The values of the volume fraction in the e 1 -structural level for the fiber and ground substance (see Fig. 2(i)) are found in pages 10 and 14 of Xu 2011, respectively). However, in the case of the values of the volume fraction at the e 2 -structural level for the fibril and ground substance see Fig. 2(ii)), to the best of our knowledge and understanding, there is no experimental data available and therefore we take approximate values having into account the transmission electron micrograph of collagen fibrils shown in Figure 3(b) of Sherman et al. [21].
The stratification functions describe the wavy effect in the meso-and nano-structures, respectively (see Figure 2).

Results and discussion
The previously developed methodology and the corresponding data for the dermis constituents allow to estimate the effective relaxation module and creep compliance of this biological tissue. In the formulation, the dermis was considered as a non-aging linear viscoelastic and hierarchical composite with isotropic constituents. The presence of a generalized periodicity at the different length scales, determined by the stratified functions q ðaÞ : R 3 ! R, with a ¼ y, z, and the consideration of x 2 , x 3 as principal axes (see Fig. 2(ii)) have an influence in the anisotropic character of the structure. In this sense, the homogenized viscoelastic tensorR ab given in Eq. (52b) has a monoclinic symmetry (13 independent viscoelastic coefficients). This result is not true when the stratification function coincides with the identity function where it is obtained instead a material with transversal isotropic symmetry. For the completeness of the analysis, we present the mathematical relationship between the effective relaxation modulus and the effective creep compliance, given in the Laplace-Carson space, where I is the fourth order identity tensor. The expression Eq. (59) is reported in Hashin [39] and it is used to compute the effective creep compliance onceR is known. The final results in the temporal space are displayed in Table 3 and Figure 3. The MATLAB's functions INVLAP and GAVSTEH developed by Hollenbeck [40] and Srigutomo [41], respectively are used in the inversion of Laplace-Carson transform. The algorithms can transform functions of complex variable s a , where a is a real exponent. They can also transform functions which contain rational, irrational and transcendent expressions. As a negative aspect, they present problems close to zero.
In Figure 3, we focus on the behavior of some chosen coefficients:R 2233 ,R 2323 ,Ĵ 2233 ,Ĵ 2323 : As observed, the curves present an asymptotic behavior, which have a good agreement with the physical relaxation and creep response in viscoelastic materials for long times. Moreover, the monotony changes and the convexities forR andĴ in the same components, for instance,R 2233 andĴ 2233 are opposite. This behavior is a consequence of the Eq. (59) and the fact thatR ijmn ðpÞ andĴ mnkl ðpÞ are inverse tensors. Similar patterns have been remarked in other works (see Ref. 42,43]).

Conclusions
A general methodology in the use of the three-scale Asymptotic Homogenization Method (AHM) is considered for modeling non-aging, quasi-periodic, hierarchical and linear viscoelastic composite materials. The present work takes Table 1. Mechanical properties for the constituents of the dermis (see Ref. [27,37]  into account the stratified functions in the homogenization procedure allowing the study with more general periodic structures. The analytical solution for the local problems and the expressions of the effective coefficients for hierarchical laminated composites with anisotropic components and perfect contact at the interfaces are derived. Also, the interconnection between the effective relaxation modulus and the effective creep compliance is performed. The approach is applied to study the overall viscoelastic behavior of the dermis. The results shown in Table 3 and Figure 3 could have a special interest in clinical applications, for instance, in the study of the epidural anesthesia procedure and in the development of haptic devices used in surgical operations. Its study is important in order to provide a location as accurate as possible for the insertion of the needle and the mechanical response of the skin (see Ref. [44]). The results can be used in the automation of the needle insertion procedure (see Ref. [45]).  The theoretical framework that we have developed here can be applied to large variety of hierarchical tissues which exhibit a viscoelastic mechanical response and are characterized by a complex, tortuous geometry, such as solid tumors, see, e.g. Netti et al. [46]; Penta and Ambrosi [47]; Ram ırez-Torres et al. [48]. A further development also includes the generalization to a viscoplastic framework by considering remodeling aspects Ram ırez-Torres et al. [31]. Although in the most general case numerical simulations cannot be avoided, application of our generalized framework would enable us to encode most of the geometrical complexity analytically, rather than enforce it directly in the setup of the computational grid, thus leading to a substantial reduction of the computational cost.
Also, the article is open to several improvements; for instance, with the framework established to layered materials, it is possible to work with more realistic geometries (e.g. fiber-reinforced composites) by using the stratified functions in the general form q ðaÞ : R 3 ! R 3 , with a ¼ y, z (see Ref. [10]). It represents a novel point of view in order to make comparisons with framework entirely developed to fiber elastic [12] and viscoelastic [16,42] reinforced composites. Besides, the results can be extended to the case of imperfect contact conditions by following the methodology in Guinovart-Sanju an et al. [35] and considering viscoelastic interfaces (see e.g. Ref. [49]) for each structural level. Another natural step is to generalize the present work to aging viscoelastic solids, see e.g. Sanahuja [50]. This way, the model will be applicable to much larger variety of biophysical (diseased and healthy) scenarios of interest.