Jingwen Zhu

Emilie-Laure Zins

Mohammad Esmaïl Alikhani
email: esmail.alikhani@sorbonne-universite.fr

Mohammad Esmail Alikhani

Dehydrocoupling of dimethylamine borane by titanocene: Elucidation of ten years of inconsistency between theoretical and experimental descriptions

Introduction

The dehydrocoupling reactions between main group elements are currently the subject of numerous studies in organometallic synthesis laboratories. [1][2][3] These reactions result in the synthesis of new inorganic compounds in particular new polymers with varied properties.

Since these reactions lead to the elimination of dihydrogen, they are also studied for possible applications in the elaboration of hydrogen storage material. 4,5 The study of the catalytic dehydrocoupling of amine-boranes fits perfectly in this double context of generation of high added-value products such as white graphene, and the controlled release of dihydrogen. [6][7][START_REF] Johnson | Synthesis and Application of Organoboron Compounds[END_REF][START_REF] Bhunya | [END_REF] Organometallic catalysis has proved particularly well suited to promote and control the dehydrogenation of dimethylamine borane (DMAB). [10][11][12][13][14] In addition, Group IV metallocenes have been among the most used catalysts for more than 30 years due to their easy synthesis and use. Since titanium is the second most abundant transition metal on Earth, titanocene can be considered as a green catalyst in a sustainable economy. 15 It is probably all these reasons that explain the numerous studies carried out on the dehydrogenation of DMAB by titanocene. In 2006, Manners' team studied experimentally this reaction for the first time. 16 They have outlined the efficiency of this catalyst for the dehydrocoupling of both DMAB and diisopropylamine borane (D i PrAB) and demonstrated that the dehydrocoupling of DMAB involves a homogeneous catalyst and leads to the formation of a cyclic diborazane (Me 2 N-BH 2 ) 2 , with the loss of two dihydrogen molecules.

However, for more sterically encumbered substrates, such as D i PrAB, the aminoborane i Pr 2 N=BH 2 is the final product of the dehydrogenation reaction.

Taking up the challenge, Luo and Ohno 17 studied the reaction pathways leading to the catalytic dehydrogenation of DMAB using density functional theory (DFT), with the B3LYP functional. They proposed two possible reaction pathways: (i) a first process involving a single molecule of DMAB (intramolecular mechanism) resulting in the formation of the aminoborane Me 2 N=BH 2 , and (ii) a second, intermolecular process involving two DMAB molecules leading to the formation of the linear dimer. From their calculated potential energy surfaces, they concluded that the intramolecular mechanism was the preferable pathway compared to the intermolecular one. To come to an agreement with the experimental finding, they suggested a further off-metal dimerization of the Me 2 N=BH 2 species, leading to the formation of the cyclic diborazane. Three years later, Manners et al. 18 In 2013, Tao and Qi reinvestigated the reaction pathway of the catalytic dehydrocoupling of the DMAB within the DFT framework. 19 Their conclusion was fully in line with the previous theoretical work and again, at odds with the experimental finding. Thus, despite two theoretical investigations, the dehydrocoupling mechanism of DMAB catalyzed by titanocene remains incompletely understood. The inconsistency between theoretical and experimental results is summarized in Scheme 1.

Scheme 1: DMAB dehydrocoupling mechanism from both experimental (Exp.) 1,18 and theoretical (DFT) 17,19 points of view. This inconsistency has also been highlighted in two recent reviews on the catalyzed dehydrogenation of amine boranes, by Weller et al. [START_REF] Johnson | Synthesis and Application of Organoboron Compounds[END_REF] and Rossin and Peruzzini. 20 Understanding the full mechanism would be an essential step towards an improvement in tailoring processes establishing new chemical bonds between p-block elements, with concomitant loss of H 2 . The present study aims at re-investigating these reaction pathways

Off-metal

On-metal combining density functional theory (DFT) and interpretative methods. To this end, we reconsider the previously proposed mechanisms using the dispersion-corrected density functional theory (DFT-D) that enables a much suitable description of non-covalent interactions. 21 A detailed electronic and bonding description of each compound involved in the reaction path will be performed within the topological frameworks. We will pay a particular attention to the effects of non-covalent interactions on the bondformation/dissociation process along the reaction.

Computational methods

All calculations were performed using the Gaussian 09 Rev D.01 quantum chemical package. 22 The geometrical structures at all stationary points were fully optimized with the popular three-parameter functional B3LYP, to which has been added the so-called Grimme empirical dispersion correction with Becke-Johnson damping (GD3BJ). 23,24 This functional is referred to as B3LYP-GD3BJ in the following. There are two main reasons for the choice of this functional: first, the reliability of this functional to describe the non-covalent interactions, 25,26 and second, to facilitate comparisons with the previous theoretical works 17,19 which used the same, B3LYP functional. All optimizations were undertaken with the same basis set as the one used in the Luo and Ohno's theoretical work: 6-31G* for cyclopentadienyl (Cp) and methyl groups, 6-31++G** for all the atoms involved in the bondformation/dissociation, and the LanL2DZ relativistic pseudo-potential for the metal atom.

Vibrational frequency calculations have been performed to obtain the ZPE (vibrational zeropoint-energy) as well as Gibbs free energies and to check the nature of the stationary point (minimum or transition state, TS).

The effect of solvent in energy was investigated along the reaction pathways by performing calculations with the PCM (Polarizable continuum) model 27 and taking toluene as solvent.

Bonding characterization was discussed using both the Electron localization Function (ELF) 28 -30 and Quantum Theory of Atoms In Molecules (QTAIM) 31 , 32 topological frameworks with the TopMod 33 and AIMAll 34 softwares, respectively. In addition, an orbital analysis has been done with the NBO 6.0 program on the optimized geometry. [START_REF] Glendening | Natural bond orbital analysis program: NBO 6.0[END_REF] In order to check our theoretical results obtained with the B3LYP-GD3BJ functional, we performed a single point calculation using the high level ab initio CCSD(T) [START_REF] Purvis | [END_REF]37 method on the optimized geometry of the first η 1 -and η 2 -complexes formed in the entrance channel of the Cp 2 Ti + H 3 B-NHMe 2 reaction.

Results and discussion

-The [Cp 2 Ti(BH 3 •NHMe 2 )] adduct: a 3-center 2-electron species

In the entrance channel, the two-body interaction between one Cp 2 Ti and one DMAB molecule may lead to the formation of two different σ complexes, in which one or two hydrogen atoms of the BH 3 moiety interact via a 3-center 2-electron (3C/2e) bonding motif with Ti center. In Figure 1 are reported some relevant geometrical parameters for both adducts. with respect to the isolated DMAB (1.202 Å). The Ti-B distance is notably shorter in the η 2complex (2.341 Å) than in η 1 -complex (3.079 Å). It is also worth noting that our geometrical parameters are rather close to those obtained with the same B3LYP functional in previous studies. 17,19 Nevertheless, the binding energy calculated with respect to free partners varies much depending on whether or not the dispersion-correction component is included in the B3LYP functional (E(η 1 -complex) = -33 & -78 kJ/mol and E(η 2 -complex) = -53 & -104 kJ/mol, for B3LYP and B3LYP-GD3BJ levels respectively). A single-point calculation at the CCSD(T) level on the structures optimized at the B3LYP level gives a similar binding energy to that obtained with B3LYP-GD3BJ (E(η 1 -complex) = -76 and E(η 2 -complex) = -105 kJ/mol). Consideration of the energy difference between B3LYP and B3LYP-GD3BJ results leads to quite contradictory conclusion when Gibbs free energy is taken into account. As shown in Figure 2, both η 1 -and η 2 -complexes are endoergic at the B3LYP level, while they are exoergic with B3LYP-GD3BJ and CCSD(T) methods. The QTAIM analysis 38 of these complexes reveal one BCP (Bond Critical Point) and one BP (Bond Path) between Ti and one hydrogen atom of BH 3 in the case of the η 1 -complex

(ρ[BCP(Ti•••H)] = 0.039 a.u. and  2 [BCP(Ti•••H)] = +0.14 a.u.
). The ELF topological analysis 39 , 40 of the protonated V(Ti,B,H) basin clearly shows the presence of a 3C/2e interaction in the η 1 -complex (population of basin V(Ti,B,H) = 2.00 e, with 5% for Ti, 13%

for B and 82% for H). For the η 2 -complex, the close contact between boron and titanium As illustrated in Figure 3, the initial step of this dehydrogenation reaction starts with a formation of the two-body η 1 -complex (labelled as S 1 ) which proceeds, via a low barrier height ( 20 kJ/mol), to an N-H oxidative addition at the Ti center, forming the [Cp 2 TiH:H 3 B-NMe 2 ] species (labelled as S 2 in Figure 3) with a large stabilization energy ( 100 kJ/mol). 3). This alternative reaction path (Route II in Figure 3) goes actually downhill with a small energy gain of about 15 kJ/mol. Consequently, the Route II is energetically more favorable than the Route I. To the best of our knowledge, the S 4 complex is optimized here for the first time along the dehydrogenation path of DMAB by Cp 2 Ti. It should be noted that the optimized S 4 complex is actually a good candidate for the hypothetical structure suggested by R. Waterman in a comprehensive review 41 in 2013 to explain and complement the Manners' mechanistic model.

The QTAIM analysis of the bonding pattern in the S 4 complex clearly reveals the formation of a dihydrogen bond-like interaction (Figure 4), with ρ[BCP(H•••H)] = 0.005 a.u. and

 2 [BCP(H•••H)] = +0
.01 a.u. between the second DMAB molecule and the hydrogen atom directly linked to the Ti center. However, a thorough charge analysis demonstrates that the net atomic charge borne by each of the two hydrogen atoms engaged in the H•••H interaction is negative, q(H(B)) = -0.68 e and q(H(Ti))= -0.34 e. As a consequence, this interaction could be labelled as X-H δ-••• δ-HML n , unlike in a classical dihydrogen bond that correspond to X- 42 The concept of "homopolar dihydrogen bond", expressed as X-H 43- 45 It drew the attention of researchers, because of the dominant role of van der Waals attractive interactions rather than the classical electrostatic interactions in the stabilization of the intermediate complexes. However, the kind of interaction identified in S 4 appears to be slightly different. To the best of our knowledge of the literature, this is the first identification of an X-H δ-••• δ-HML n homopolar dihydrogen bond. Other dihydrogen bonds are also identified in this intermediate. B-H•••H-C interactions between the BH 3 moiety and methyl groups or Cp ligands (which are not presented in Figure 4) have electron density in the [0.004 a.u., 0.006 a.u.] ranges.

H δ+ ••• δ-H-Y pattern.
δ-••• δ- H-X or X-H δ+ ••• δ+ H-X, is a well-known concept for B-H•••H-B or C-H•••H-C interactions.
From S 4 one, the formation of the precursor of the linear dimer and titanium hydride (S 5 ) is an exothermic process, which releases approximately 25 kJ/mol. The formation of this intermediate involves the cleavage of the Ti-N bond and a concerted Ti-H-B-N ring opening process. Concomitantly, as shown in Figure 4, the electron density at BCP(Ti•••H) decreases from 0.063 to 0.036 a.u. in going from S 4 to S 5 . An ELF atomic contribution showed that the Ti•••H bonding should be considered as a 3C/2e interaction (V(Ti,B,H) = 1.93 e, with 3% for Ti, 13% for B and 84% for H). Furthermore, a dihydrogen bond is identified in S 5 (labelled as BCP(H•••H) in Figure 4), in which the net atomic charge of the hydrogen atoms are of opposite signs. The electron density at this bond critical point amounts to 0.045 a.u., slightly larger than that of the Ti•••H 3C/2e bond. In line with the outcomes of the discussion on the energy profile along the reaction path 1, we could therefore suggest a full mechanistic cycle for the catalytic formation of a linear diborazane and release of an equivalent H 2 (Figure 5). The suggest mechanism is in agreement with the experimental scheme proposed by Manners' team. Such a catalytic cycle involving a three-body interaction was already proposed in one of the previous theoretical work, 17 but had been ruled out, because of an endoergic reaction in the entrance channel with a high barrier height ( 130 kJ/mol) at the B3LYP level of theory.

In the present work, owing to the use of the dispersion-corrected DFT approach, this path has been calculated to be both energetically and kinetically favorable. As show in Figure 6, the first step corresponds to the formation of a three-body complex, consisting of an η 2 -complex and the addition of a second DMAB to the η 2 -complex, leading to a new complex with a binding energy about 164 kJ/mol (labelled as S 6 in Figure 6). From the QTAIM analysis of the S 6 complex, a weak intermolecular interaction is identified between the titanium atom (q(Ti) = +1.34 e) and the hydrogen atom (q(H(N)) = +0. the NBO (Natural Bond Orbital) analysis, which provides a measure of the electron donoracceptor delocalization. 46 As illustrated in Figure 7, this interaction should be considered as a donor-acceptor interaction from a lone pair orbital of titanium ( ) toward an anti-bonding orbital of N-H (σ*), producing a stabilization of 6.3 kcal/mol.

Figure 7: Donor-acceptor interaction in the S 6 complex.

Transformation of S 6 to S 7 is characterized by a hydrogen transfer from the nitrogen atom to the metal center, which leads to the formation of titanium hydride. From energetic point of view, the S 7 intermediate evolves towards a new complex (labelled as S 8 in Figure 6) without any energy demand. In the latter complex, dihydride titanium interacts with the linear diborazane fragment via two dihydrogen bonds (ρ[BCP(Ti-H•••H-C)] = 0.011 a.u. and

 2 [BCP(Ti-H•••H-C)] = +0.

-Path 3: From the linear diborazane to the cyclic diborazane

Experimental data demonstrated that the cyclization of the linear diborazane is a heterogeneously catalyzed process. We thus computationally investigated an "on-metal" cyclization of the linear dimer. The calculated energy profiles of this reaction pathway are presented in Figure 9. A similar mechanism was already investigated by Tao and Qi at the DFT level. 19 Their results suggested that this mechanism is highly endothermic. On the contrary, our results obtained within the DFT-D framework demonstrate that this dehydrocoupling pathway is an energetically feasible process under the experimental conditions (Figure 9). This reaction pathway involves three intermediates, named S 9 , S 10 and S 11 , as well as two transition states, named TS 9_10 and TS 10_11 . slightly smaller than its counterpart in S 8 (0.032 a.u., Figure 8).

Over the low barrier height, the cleavage of the N-H bond leads to a new intermediate (labelled as S 10 in Figure 9), in which we note an asymmetric η 2 3C/2e interactions between the titanocene and the amidoborane fragment, characterized by two distinct BCPs (ρ = 0.056 and 0.062 a.u.). The catalytic cycle is completed with a second hydrogen transfer from borane to the metal center, producing the so-called cyclic diborazane (labelled as S 11 in Figure 8).

The QTAIM and ELF analyses suggest that the cyclic diborazane is connected to the dihydride titanocene via a 3C/2e interaction (ρ = 0.037 a.u., V(Ti,B,H) = 1.97 e, with 3% for Ti, 15% for B and 82% for H) and five weaker dihydrogen bonds (0.004 < ρ < 0.008 a.u.). As interim first conclusion, we summed up the three paths studied above in Scheme 2, which consistently allowed us to describe and complement the experimental considerations.

Scheme 2: A bridged schematic representation of the catalytic formation of the cyclic diborazane.

-Substitution effect on the terminal product

Experimental investigations also showed a fundamentally different reaction path in case of a bulky substituant on the nitrogen, such as diisopropylamine borane (D i PrAB). Indeed, the dehydrocoupling of the D i PrAB leads to the formation of the so-called aminoborane i Pr 2 N=BH 2 product. Our computational results enable us to explain the singularity of D i PrAB dehydrogenation path in two ways.

First, a careful look to the energetics reported in Figure 10 reveals that in contrast to the case of DMAB, the "Route I" transforming the S 2 structure to the aminoborane (S 3 ) is more favorable than the "Route II" for the D i PrAB. In addition with the geometrical parameters reported in Figure 11, we would like to mention that in the S 4 compound, we identified several dihydrogen bonds between bulky groups of two D i PrAB units, while the dihydrogen bonds observed with DMAB molecules were found between the borane part of the second DMAB and the NMe 2 subunit of the amidoborane.

Every dihydrogen bond only presents a low electron density ( 0.005 a.u.) at the identified bond critical points for both DMAB and D i PrAB complexes.

Consequently, the dehydrogenation of sterically bulky groups at amine results in the formation of aminoborane rather than diborazane.

Conclusion

Theoretical calculations within the DFT-D framework make it possible to obtain a reaction pathway in perfect agreement with the experimental data: the cyclization of the intermediate occurs via an "on-metal" mechanism. A major conclusion of this study is that the formation of a van der Waals complex in which a second DMAB molecule interacts with amidoborane subunit via dihydrogen bonds appears to play a central role in the formation of the linear dimer intermediate. The topological (ELF and QTAIM) and orbital (NBO) analyses allowed us to characterize the natures of non-covalent interactions in the reaction intermediates.

Noticeably, we underlined the important role played by weak interactions, such as dihydrogen bonds and 3C/2e interactions, along the catalyzed dehydrocoupling process of DMAB. For a consistent description of these interactions, the coverage of dispersion contribution is fundamental.

Finally, it was also shown that the bulky substituents, as in the D i PrAB, prevents the formation of a van der Waals complex leaving open only the path to the formation of a diisopropyl aminoborane.

  published out an extensive experimental study to elucidate the reaction mechanism of titanocene-catalyzed dehydrogenation of amine-boranes. In the case of the dehydrogenation of DMAB, their study has highlighted the role of the HNMe 2 -BH 2 -Me 2 N-BH 3 linear diborazane as a reaction intermediate prior to the formation of the cyclic diborazane (Me 2 N-BH 2 ) 2 . They also demonstrated that this dehydrogenation process involves two consecutive catalytic cycles:  At first, they proposed a two-step cycle leading to a coordination complex [Cp 2 Ti(η 2 -BH 3 •NHMe 2 )] involving a 3-center 2-electron (3C/2e) Ti-H-B interaction. This complex further undergoes a proton transfer. When adding a second molecule of DMAB, another two-body reaction is expected to lead to the formation of the linear dimer, with the concomitant loss of one H 2 molecule.  The second cycle consists in the dehydrocoupling of the linear dimer resulting in the cyclic diborazane upon Cp 2 Ti catalytic process.

Figure 1 :

 1 Figure 1: Some important geometrical parameters of the [Cp 2 Ti(BH 3 •NHMe 2 )] complexes.

Figure 2 :

 2 Figure 2: Binding energies (relative Gibbs free energies) calculated at three different levels of theory.

  atoms results in the formation of a BCP (ρ[BCP(Ti-B) = 0.055 a.u. and  2 [BCP(Ti-B)] = +0.13 a.u.) with high ellipticity ( 16) indicating the easy and probable evolution of this point to the formation of two BCP connecting the Ti center to two hydrogen atoms of BH 3 . This feature is clearly revealed by the ELF atomic contribution of these two protonated basins.Owing the positive sign of the Laplacian charge density at bond critical point, the Ti-H and Ti-B interactions belong to the closed-shell interaction.

  Two bond critical points are characterized for this intermediate: a BCP located between Ti and N atoms and another between Ti and H of BH 3 subunit. The former indicates the formation of a Ti-N bond with ρ[BCP(Ti-N)] = 0.059 a.u. and  2 ρ[BCP(Ti-N)] = +0.17 a.u.. The second BCP is characteristic of a 3C/2e interaction (ρ[BCP(Ti•••H)] = 0.064 a.u. and  2 [BCP(Ti•••H)] = +0.13 a.u.).

Figure 3 :

 3 Figure 3: Energy profile along the first reaction pathway. Relative gas-phase zero-point corrected energies and solvent effect-corrections (values in parentheses) are given in kJ/mol.

Figure 4 :

 4 Figure 4: Bond Critical Points and AIM atomic charges in S 4 and S 5 complexes.

18

 18

Figure 5 : 2 :

 52 Figure 5: First catalytic cycle from [Cp 2 Ti(η 1 -BH 3 •NHMe 2 )] adduct to linear diborazane via a dihydrogen-bonded complex. Relative energies calculated in the gas-phase including ZPE correction in kJ/mol are mentioned.

Figure 6 :

 6 Figure 6: Energy profile along the second reaction pathway. Relative gas-phase zero-point corrected energies and solvent effect-corrections (values in parentheses) are given in kJ/mol.

  44 e) linked to the nitrogen atom of the second DMAB molecule (ρ[BCP(Ti•••H-N)] = 0.017 a.u. and  2 ρ[BCP(H•••H-N)] = +0.01 a.u.). The ELF analysis of the atomic contributions on the corresponding protonated basin clearly shows that the present interaction could not be considered as a 3C/2e one. We further investigated the nature of the interaction by means of

  03 a.u.). As illustrated in Figure8, we identified three bond critical points at the N-B bonds with electron density of 0.122, 0.140, and 0.122 a.u.. It is interesting to note the presence of a dihydrogen bond (labelled as N-H•••H-B in Figure8) in the linear diborazane fragment (ρ[BCP(N-H•••H-B)] = 0.032 a.u.), which seems to be the likely cause of the closing of the six-membered ring in a chair-like conformation.

Figure 8 :

 8 Figure 8: Bond Critical Points and AIM atomic charges in S 8 complex. Accordingly, this three-body mechanism beginning with the [Cp 2 Ti(η 2 -BH 3 •NHMe 2 )] adduct is calculated to be energetically similar to the first catalytic cycle induced by the formation of the [Cp 2 Ti(η 1 -BH 3 •NHMe 2 )] adduct. Thus, within the DFT-D framework, two competitive reaction pathways are found to lead to the linear diborazane intermediate that should further undergo an "on-metal" cyclization, according to the experimental findings.

8 Figure 9 :

 89 Figure 9: Energy profile for the titanocene-catalyzed dehydrocoupling of linear diborazane. Relative gas-phase zero-point corrected energies and energies taking into account of solvent effect (values in parentheses) are shown in kJ/mol.

Figure 10 :

 10 Figure10: Relative Gibbs free energy of four structures along the path 1 (see Figure3) for both DMAB (full squares, dotted lines) and D i PrAB (full circles, full lines) compounds.

Figure 11 :

 11 Figure 11: Optimized structure of van der Waals complex labelled as S 4 with (a) DMAB and (b) D i PrAB compounds.

[Cp 2 Ti] + DMAB Exp. DFT η 2 -Complex η 1 -Complex N-H cleavage N-H and B-H cleavage [Cp 2 TiH 2 ] H 2 B=NMe 2 [Cp 2 Ti] + DMAB + H 2 [Cp 2 TiH 2 ] [Cp 2 Ti] + H 2 Linear diborazane Cyclic diborazane -H 2


Oscar Luis

Cruz González
email: gonzalez@lma.cnrs-mrs.fretlebon@lma.cnrs-mrs.fr

Ariel Ramírez Torres

Reinaldo Rodríguez Ramos
email: reinaldo@matcom.uh.cu

Frédéric Lebon

Homogénéisation triple échelle de matériaux composites viscoélastiques Three scales asymptotic homogenization of viscoelastic composite materials

Keywords: Échelles multiples, Homogénéisation asymptotique, Viscoélasticité linéaire, Derme Multiple scales, Asymptotic homogenization, Linear viscoelasticity, Dermis 1

Une méthode d'homogénéisation asymptotique à plusieurs échelles est proposée pour estimer les propriétés viscoélastiques effectives d'un matériau composite hiérarchique. Les différents niveaux structurels sont analysés sous l'hypothèse d'une périodicité généralisée et se caractérisent par l'utilisation de fonctions stratifiées. Le matériau est modélisé comme un composite viscoélastique linéaire non vieillissant avec une structure hiérarchique en couches. Les problèmes locaux associés, le problème homogénéisé et les expressions des coefficients effectifs sont dérivés pour chaque niveau d'organisation en utilisant le principe de correspondance et la transformation de Laplace-Carson. En travaillant avec une fonction stratifiée générale, des composants anisotropes et un contact parfait aux interfaces, les problèmes locaux et les coefficients effectifs sont calculés analytiquement dans l'espace de Laplace-Carson. L'inversion numérique dans l'espace temporel d'origine est également effectuée. Une approche pour la modélisation du derme en tant que matériau composite viscoélastique hiérarchique est proposée et le calcul du module de relaxation effectif est effectué.

Introduction

The performance of mechanical properties as weight, heat resistance, corrosion, among others are optimized thanks to the use of composite materials. Specifically, those with viscoelastic properties are widely applied in the aerospace, aeronautical and automobile industry, in bioengineering, as well as in the design of durable and sustainable structural components. Some of these composites, such as biological and synthetic viscoelastic materials, often possess a hierarchical structure and exhibit a multiscale character. The modeling of composite materials requires the development of micromechanics techniques to predict the general (or effective) properties of the heterogeneous structure from the properties, density, proportion and arrangement of its constituents. An excellent review on these methods can be found in [START_REF] Kalamkarov | Asymptotic Homogenization of Composite Materials and Structures[END_REF] and [START_REF] Sevostianov | Generalization of Maxwell homogenization scheme for elastic material containing inhomogeneities of diverse shape[END_REF]. The existence of a sequence of scales of decreasing order allows to perform a succession of homogenization steps (see [START_REF] Yang | A three-scale homogenization algorithm for coupled conduction-radiation problems in porous materials with multiple configurations[END_REF], [START_REF] Yang | High-order three-scale method for mechanical behavior analysis of composite structures with multiple periodic configurations[END_REF]). In every step, the homogenization scheme requires the solution of a cell problem with data corresponding to the homogenized material properties of the previous step. On the other hand, many complex heterogeneous structures are characterized by more general periodic functions (see [START_REF] Tsalis | Charalambakis Homogenization of structures with generalized periodicity[END_REF]). These functions, called of stratification, describe the microstructure of the composite material and they are related to homogenization problems of shell-type structures of widely applied technological interest (nano-hulls, fiber-reinforced polymers (FRP), civil engineering structures repair, modeling of human heart tissue, among others). Actually, the investigation of the effective properties of non-ageing viscoelastic composites are mainly based on the correspondence principle and the Laplace transform (see [START_REF] Hashin | Viscoelastic behavior of heterogeneous media[END_REF]). The procedure, see for example [START_REF] To | Overall viscoelastic properties of 2D and twophase periodic composites constituted of elliptical and rectangular heterogeneities[END_REF], consists into the change of the convolution constitutive law, which describes the nonageing viscoelastic behavior, to a fictitious elastic one in the Laplace domain. Then, the inversion of Laplace transform is considered to derive the effective behavior in the time domain. The present work deals with three-scale Asymptotic Homogenization Method (AHM) which is proposed to modeling a non-ageing linear viscoelastic composite material with generalized periodicity and two hierarchical levels of organization, characterized by the small parameters ε 1 and ε 2 (see Fig. 1). The theoretical aspects of the method are rigorously developed by [START_REF] Allaire | Multiscale convergence and reiterated homogenisation[END_REF], [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF], [START_REF] Trucu | Three-scale convergence for processes in heterogeneous media Applicable Analysis[END_REF]. The approach aims to generalize the results for elastic composite materials obtained in [START_REF] Ramírez-Torres | Three scales asymptotic homogenization and its application to layered hierarchical hard tissues[END_REF] and to extend them to non-ageing viscoelastic ones with generalized periodicity at each scale of different length. The method is finally applied to investigate the effective viscoelastic properties of the dermis.

The linear viscoelastic problem for hierarchical structures

Let us denote by Ω ∈ R 3 a linear viscoelastic heterogeneous material possessing two hierarchical levels of organization and a generalized periodicity (see more details of the concept in [START_REF] Tsalis | Charalambakis Homogenization of structures with generalized periodicity[END_REF], [START_REF] Tsalis | Homogenization of elastoplastic composites with generalized periodicity in the microstructure[END_REF]). Three different scales, namely d 1 , d 2 and L, are considered. They characterize the different sizes of the structural levels and are assumed to be well-separated, i.e.

ε 1 = d 1 L << 1 and ε 2 = d 2 L << ε 1 . (Eq. 1)
At the ε 1 -structural level, the domain Ω is occupied by a two-phases quasi-periodic composite such that

Ω = Ω ε 1 1 ∪ Ω ε 1 2 , Ω ε 1 1 ∩ Ω ε 1 2 = / 0. We assume that Ω ε 1 1 = ∪ N 1 α=1 α Ω ε 1
1 and the interface between Ω ε 1 1

and Ω ε 1 2 is denoted by Γ ε 1 . Here, Y denotes the unit cell. At the ε 2 -structural level, we consider that for each α = 1, ..., N 1 , α Ω ε 1 1 is a two-phase quasi-periodic composite material. Hence, we define, α Ω

ε 1 1 = Ω ε 2 1 ∪ Ω ε 2 2 , Ω ε 2 1 ∩ Ω ε 2 2 = / 0, Ω ε 2 1 = ∪ N 2 β=1 β Ω ε 2
1 and the interface between Ω ε 2 1 and Ω ε 2 2 is denoted by Γ ε 2 . In addition, Z is the corresponding unit cell. Besides, x x x (x i ) represents the global spatial coordinate and two local and independent variables are introduced

y y y = ρ ρ ρ (y) (x x x) ε 1 and z z z = ρ ρ ρ (z) (x x x) ε 2 , (Eq. 2)
The functions ρ ρ ρ (i) : R 3 -→ R 3 are the so-called stratified functions. The parametric equations ρ ρ ρ (i) (x x x) = constant with i = y, z describe the ε 1 -and ε 2 -structural levels, respectively. Besides, the stratified functions satisfy ρ ρ ρ (i) ∈ C ∞ (Ω) (see [START_REF] Tsalis | Charalambakis Homogenization of structures with generalized periodicity[END_REF]).

Considering that the constitutive response of all the constituents of the composite body is linear viscoelastic, and ignoring inertial terms, the problem in Ω consists in, where σ σ σ (σ i j ) represents the second-order stress tensor, ξ ξ ξ (ξ i j ) denotes the second-order strain tensor,

∇ •σ σ σ(x x x,t) + f f f (x x x,t) = 0 0 0 in (Ω\(Γ ε 1 ∪ Γ ε 2 )) × R (Eq. 3) u u u(x x x,t) = ū u u on ∂Ω d × R σ σ σ(x x x,t) •n n n = S S S on ∂Ω t × R u u u(x x x,t) = 0 0 0 in Ω × {0}
u u u (u i ) is the viscoelastic displacement, f f f ( f i ) denotes the action of external volume forces that satisfies f f f (x x x,t) ∈ L 2 (Ω × R).
Moreover, ū u u ( ūi ) and S S S ( Si ) are the prescribed displacement and traction on the boundary ∂Ω = ∂Ω d ∪ ∂Ω n , with ∂Ω d ∩ ∂Ω n = / 0 and n n n (n i ) is the outward unit vector normal to the surface ∂Ω. Additionally, the scale-dependent constitutive law which relates the second-order stress and strain tensors is written as follow, (see [START_REF] Christensen | Theory of Viscoelasticity-2nd Edition An Introduction[END_REF])

σ σ σ(x x x,t) = t 0 R R R (x x x,y y y,z z z,t -τ): ∂ξ ξ ξ(u u u(x x x, τ)) ∂τ dτ. (Eq. 4)
Here, the relaxation modulus is denoted by R R R (R i jkl ) and it fulfills the following symmetry properties

R i jkl = R jikl = R i jlk = R kli j . Also, we assume R R R ∈ L ∞ (Ω × R)
and that it's positively defined. Continuity conditions for displacement and traction are imposed on both Γ ε 1 and Γ ε 2 , i.e.

[[u u u(x x x,t)]] = 0 0 0, σ σ σ(x x x,t) •n n n ( j) = 0 0 0, ( j = y, z) (Eq. 5)

wherein the outward unit vectors to the surfaces Γ ε 1 and Γ ε 2 are represented by n n n (y) = (n (y)

1 , n (y) 2 , n (y)
3 ) and n n n (z) = (n

(z) 1 , n (z) 2 , n (z)
3 ), respectively (see Fig. 1). The operator [[•]] denotes the jump across the interface between the two constituents. In the context of small displacements, the components of the second-order strain tensor ξ ξ ξ are given by the relation

ξ kl (u u u (x x x,t)) = 1 2 ∂u k (x x x,t) ∂x l + ∂u l (x x x,t) ∂x k . (Eq. 6)
The statement of the scale-dependent constitutive law (Eq. 4) corresponds to the special form of nonageing linear viscoelastic materials (see [START_REF] Maghous | Periodic homogenization in thermoviscoelasticity : case of multilayered media with ageing[END_REF]). Therefore, the problem can be transformed into an elastic one using the Laplace-Carson transform. The aforementioned is known as the correspondence principle. Applying the Laplace-Carson transform to (Eq. 3)-(Eq. 5) and taking into account the summation convention, the mathematical statement for the linear quasi-static viscoelastic heterogeneous problems in the Laplace-Carson space is written,

∂ ∂x j R i jkl (x x x,y y y,z z z, p) ξ kl (u u u (x x x, p)) + f i (x x x, p) = 0, in (Ω\(Γ ε 1 ∪ Γ ε 2 )) × [0, +∞).
(Eq. 7)

The corresponding boundary conditions associated to (Eq. 7) are

u i (x x x, p) = u i , on ∂Ω d × [0, +∞),
(Eq. 8) R i jkl (x x x,y y y,x x x, p) ξ kl (u u u (x x x, p)) n j = S i , on ∂Ω n × [0, +∞), (Eq. 9)

and the initial condition is

u i (x x x, p) = 0, in Ω × {0} .
(Eq. 10)

The interface contact conditions become,

[[u i (x x x, p)]] = 0, R i jkl (x x x,y y y,z z z, p) ξ kl (u u u (x x x, p)) n (m) j = 0 (m = y, z) . (Eq. 11)
In order to avoid problems with the notation and make the work clearer, from now on, the homogenization process is developed in the Laplace-Carson space and the functions that depend on the parameter p (e.g. Φ(p)) are indicating that we are working on that space.

3.

The three-scale asymptotic homogenization method to solve the heterogeneous problem

The aim of this section is to solve the heterogeneous problem (Eq. 7)-(Eq. 11) by using AHM. As the material property is regular in x x x, and periodic in y y y and z z z then, according to the chain rule, the derivative in relation to the global coordinate applied to any function in the form Φ Φ Φ (x x x,y y y,z z z) yields the transformation

∂Φ i ∂x j → ∂Φ i ∂x j + 1 ε 1 ∂ρ (y) l ∂x j ∂Φ i ∂y l + 1 ε 2 ∂ρ (z) m ∂x j ∂Φ i ∂z m , (Eq. 12)
and using a similar idea, the equation (Eq. 6) becomes,

ξ kl Φ Φ Φ x x x, ρ ρ ρ (y) (x x x) ε 1 , ρ ρ ρ (z) (x x x) ε 2 = ξ kl (Φ Φ Φ (x x x,y y y,z z z)) + ε -1 1 ξ y kl (Φ Φ Φ (x x x,y y y,z z z)) + ε -1 2 ξ z kl (Φ Φ Φ (x x
ξ (z) kl (Φ Φ Φ (x x x,y y y,z z z)) = 1 2 ∂ρ (z) r (x x x) ∂x l ∂Φ k (x x x,y y y,z z z) ∂z r + ∂ρ (z) s (x x x) ∂x k ∂Φ l (x x x,y y y,z z z) ∂z s .
(Eq. 15)

Homogenization procedure

The solution to the problem (Eq. 7)-(Eq. 11) is proposed as follows,

u u u ε (x x x, p) = ũ u u (0) (x x x,y y y,z z z, p) + ∞ ∑ i=1 ũ u u (i) (x x x,y y y,z z z, p)ε i 2 , (Eq. 16)
where ũ u u (0) is defined as ũ u u (0) (x x x,y y y,z z z, p) = u u u (0) (x x x,y y y,z z z, p) + ∞ ∑ i=1 u u u (i) (x x x,y y y,z z z, p)ε i 1 .

(Eq. 17)

Here and subsequently (unless necessary), the variable dependence is dropped out for convenience. Then, replacing expansion (Eq. 16) into (Eq. 7) and (Eq. 11), grouping in powers of ε 2 and considering the solvability condition reported in [START_REF] Bakhvalov | Homogenisation : Averaging processes in periodic media[END_REF], the sequence of problems depicted below can be solved in recursive form. A summary for each problem is proposed.

Problem for ε -2 2 ∂ρ (z) m ∂x j ∂ ∂z m R i jkl ξ (z) kl ũ u u (0) = 0, in (Z\Γ z ) × [0, +∞), (Eq. 18) ũ(0) i = 0, on Γ z × [0, +∞), (Eq. 19) R i jkl ξ (z) kl ũ u u (0) n (z) j = 0, on Γ z × [0, +∞).
(Eq. 20)

Since the right hand side of (Eq. 18) is zero, the solvability conditions are satisfied and ũ u u (0) is a solution of (Eq. 18) if and only if it is constant in relation to the variable z z z, i.e. ũ u u (0) = ũ u u (0) (x x x,y y y, p) (Eq. 21)

and taking into account (Eq. 17), it is obtain u u u (0) = u u u (0) (x x x,y y y, p) , u u u (i) = u u u (i) (x x x,y y y, p) .

Problem for ε -1 2 Using (Eq. 21) and the fact that ξ (z) kl ũ u u (0) = 0, the problem leads to

∂ρ (z) m ∂x j ∂ ∂z m R i jkl ξ (z) kl ũ u u (1) + ξ kl ũ u u (0) + ε -1 1 ξ (y) kl ũ u u (0) = 0 (Eq. 22) in Z\Γ z × [0, +∞) ũ(1) i = 0, (Eq. 23) R i jkl ξ (z) kl ũ u u (1) + ξ kl ũ u u (0) + ε -1 1 ξ (y) kl ũ u u (0) n (z) j = 0, (Eq. 24) on Γ z × [0, +∞).
The divergence theorem and the z-periodicity condition of R i jkl ensure the fulfillment of the solvability conditions in (Eq. 22). Thus, the existence and uniqueness of a solution for the problem (Eq. 22)-(Eq. 24) is guaranteed (see [START_REF] Bakhvalov | Homogenisation : Averaging processes in periodic media[END_REF]). Applying separation of variables, a general solution for (Eq. 22)-(Eq. 24) is given by ũ(1) m (x x x,y y y,z z z, p) = χklm (x x x,y y y,z z z, p) Ũ(0) kl (x x x,y y y, p), (Eq. 25)

where Ũ(0) kl (x x x,y y y, p) = ξ kl ũ u u (0) (x x x,y y y, p) + ε -1 1 ξ

(y) kl ũ u u (0) (x x x,y y y, p) .

(Eq. 26)

The third rank tensor χklm is y-and z-periodic. Then, substituting (Eq. 25) into (Eq. 22)-(Eq. 24) and after some simplifications the following auxiliary problem, which we call the ε 2 -cell problem, is obtained

∂ρ (z) m ∂x j ∂ ∂z m R i jkl + R i jpq ξ (z) pq ( χ χ χ kl ) = 0, in (Z\Γ z ) × [0, +∞), (Eq. 27) [ [ χklm ] ] = 0, on Γ z × [0, +∞),
(Eq. 28)

R i jkl + R i jpq ξ (z) pq ( χ χ χ kl ) n (z) j = 0, on Γ z × [0, +∞).
(Eq. 29)

Before proceeding, let's introduce first the following notation. Since the quantities involved vary on the y and z scales, the following cell average operators are defined

• y = 1 |Y| Y (•)dy and • z = 1 |Z| Z (•)dz, (Eq. 30)
where |Y| and |Z| represent the periodic cell volumes.

Problem for ε 0 2

The expressions are grouped in powers of ε 0 2 . Applying the average operator • z , taking into account the divergence theorem and the z-periodicity of the involved functions, we obtain

∂ ∂x j + ε -1 1 ∂ρ (y) n ∂x j ∂ ∂y n Ři jkl Ũ(0) kl + f i = 0 in Ω ε 1 1 × [0, +∞), (Eq. 31)
where

Ři jkl = R i jkl + R i jpq ξ (z) pq ( χ χ χ kl ) z (Eq. 32)
is the effective coefficient at the ε 1 -structural level. Note that Ři jkl = Ři jkl (x x x,y y y, p).

We have carried out the analysis for the ε 1 -structural level. In order to find the effective behaviour of the hierarchical composite, a similar procedure is performed. Then, using relations (Eq. 17) and (Eq. 26) into equation (Eq. 31) and grouping in powers of ε 1 , the following theoretical results are obtained, Problem for ε -2

1
The solvability condition is also satisfied and the following result is reached

u u u (0) = u u u (0) (x x x, p) . (Eq. 33)
Problem for ε -1

1

A general solution for the problem is given by u

m (x x x,y y y, p) = χ klm (x x x,y y y, p) ξ kl u u u (0) (x x x, p) (Eq. 34)

where the third order tensor χ klm is y-periodic and is the solution of the following problem, referred to as the ε 1 -cell problem,

∂ρ (y) m ∂x j ∂ ∂y m Ři jkl + Ři jpq ξ (y) pq (χ χ χ kl ) = 0 in Y \Γ Y × [0, +∞), (Eq. 35) [ [χ klm ] ] = 0 on Γ Y × [0, +∞), (Eq. 36) Ři jkl + Ři jpq ξ (y) pq (χ χ χ kl ) n (y) j = 0 on Γ Y × [0, +∞).
(Eq. 37)

Problem for ε 0 1

The homogenized problem is obtained and it can be written in the form

∂ ∂x j Ri jkl ξ kl u u u (0) + f i = 0 in Ω × [0, +∞) (Eq. 38)
where Ri jkl = Ři jkl + Ři jpq ξ is the expressions for the effective relaxation modulus of the hierarchical composite material. Finally, from (Eq. 8)-(Eq. 10) and applying the cell average operators over Z and Y, the boundary conditions for (Eq. 38) are

Ri jkl ξ kl u u u (0) n j = S i on ∂Ω n × [0, +∞), (Eq. 40)

u (0) i = u i on ∂Ω d × [0, +∞), (Eq. 41)
and the initial condition u

(0) i = 0 in Ω × {0} .
(Eq. 42)

4.

Effective coefficients for hierarchical laminated composites with generalized periodicity Laminated composites materials can be described by using stratified function in the form ρ : R n → R with n = 2, 3 (see [START_REF] Tsalis | Charalambakis Homogenization of structures with generalized periodicity[END_REF]). Now, we consider the general case when the stratified functions are ρ (i) : R 3 → R with i = y, z. Hence, following the same procedure carried out by [START_REF] Guinovart-Sanjuán | Analysis of effective elastic properties for shell with complex geometrical shapes[END_REF], the local problems (Eq. 27)-(Eq. 29) and (Eq. 35)-(Eq. 37) and the effective coefficients (Eq. 32) and (Eq. 39) are transformed using the Voigt's notation and can be solved analytically in recursive form.

In the process, we integrate each equation in relation to the local variables and determining the constants of integration by using the corresponding cell average operator. It is really worthy to remark that the components of the local problems are raised for the anisotropic case which is an advantage of this approach.

5.

An approach for modeling the mechanical properties of the dermis Nowadays, many researches are focus on the study of the components and the structure of the skin.

A better understanding of this soft biological tissue has a real impact on biomedical applications and also inspires modern technology such as flexible electronics, soft robotics and prosthetics (see [START_REF] Sherman | Structural characterization and viscoelastic constitutive modeling of skin[END_REF]). Skin has three main layers : the epidermis, the dermis and the hypodermis. The most important mechanical and thermal unit of the skin is the dermis ; it represents the 90% of the thickness of the skin (see [START_REF] Benítez | The mechanical behavior of skin : Structures and models for the finite element analysis[END_REF]). Also, the dermis is considered a multilayer collagen-reinforced structure. Two of their principal components are the collagen and the ground substance (supporting matrix). The collagen contributes to 75% of the fat-free dry mass and 18%-30% of the volume of dermis and can be considered like a linear elastic material. 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 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 [START_REF] Xu | Introduction to Skin Bio-Thermo Mechanics and Thermal Pain[END_REF]) In this section, an approach for modeling the dermis as an hierarchical viscoelastic composite material is proposed. We follow the structural ideas depicted in the Fig 3 of [START_REF] Sherman | Structural characterization and viscoelastic constitutive modeling of skin[END_REF]. It is stated that the stiffness of collagenous materials, which are within the ground substance [START_REF] Joodaki | Skin mechanical properties and modeling : A review[END_REF] and form the dermis, is a consequence of the properties, arrangement, and geometric distribution of collagen fibrils. These structures are long strands with wavy effects and are arranged forming wavy parallel fibers, which are flattened in the plane of the dermis and determine its properties. We present the problematic from the point of view of laminated materials (see Fig 2). Mechanical properties of materials can be found in Tab. 1 (see [START_REF] Panchal | Vibration analysis of healthy skin : toward a noninvasive skin diagnosis methodology[END_REF]). Besides, the viscoelastic constituent (ground substance) can be modeled using normalized Prony series,

G(t) = G ∞ + N ∑ n=1
G n e (-t/τ n ) , (Eq. 43)

where τ n are the relaxation times, G n are the modulus coefficients, = G ∞ is the long-term modulus (Tab. 2). For sake of simplicity in the model, only two term in the Prony series (see [START_REF] Joodaki | Skin mechanical properties and modeling : A review[END_REF]) are considered.

Young modulus (GPa) The average operator is calculated

f = V 1 f (1) +V 2 f (2) ,
(Eq. 44)

where the subscripts (1), ( 2

Fig. 1 .

 1 Fig. 1. (a) Macroscale : viscoelastic heterogeneous material. (b) ε 1 -structural level. (c) Mesoscale : quasi-periodic cell. (d) ε 2 -structural level. (e) Microscale : quasi-periodic cell. The inclusions do not intersect the boundaries.

m

  x,y y y,z z z)) , (Eq. 13) where ξ (y) kl (Φ Φ Φ (x x x,y y y,z z z)) (x x x) ∂x l ∂Φ k (x x x,y y y,z z z) ∂y m + ∂ρ (y) n (x x x) ∂x k ∂Φ l (x x x,y y y,z z z) ∂y n , (Eq. 14)

  (y) pq (χ χ χ kl ) y .(Eq. 39)

  ) are indicating the corresponding material and V i represents the volume fractions of each constituent. The outcomes in the calculation of the effective relaxation modulus are displayed in Fig 3. The methodology allows to estimate the effective behavior for the dermis.

Fig. 2 .Fig. 3 .

 23 Fig. 2. The hierarchical structure of the dermis.

  Tab. 1. Mechanical properties for the constituents of the dermis. Tab. 2. Coefficients of the Prony series.

					Poisson ratio
	Collagen	1			0.48
	Ground of substance	0.0102		0.48
		G 1	G 2	τ 1	τ 1
	Ground of substance 0.22 0.28 0.31 6.96

Acknowledgements

OLCG kindly acknowledges 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 for its support. RRR thanks to Departamento de Matematicas y Mecanica, IIMAS and PREI-DGAPA at UNAM, for its support to his research project. ART kindly acknowledges the Dipartimento di Scienze Matematiche (DISMA) "G.L. Lagrange" of the Politecnico di Torino, "Dipartimento di Eccellenza 2018-2022" ("Department of Excellence 2018-2022").