GANGHOFFER, to whom this work owes much, and to whom I owe more humanely. I would like to salute here his great scientific competence, his listening and his rigor. At his contact I learned the virtues of doubt, the requirement of precision, the overcoming of certainties too quickly acquired. If this work could be carried out, it is thanks to his precious help and his constant encouragement. Let every student have a leader like him. I would also thank Mr. Hassan Lakiss for his confidence in me. I appreciated the great freedom of autonomy that he gave me during this work. I would then like to thank the members of my jury for having done me the honor of evaluating my work. I would also like to thank the people whom I call "resources" in my thesis and who have allowed me to better understand the current and past functioning of this field. During these three years, I had the pleasure of exchanging scientific work with Yosra Rahali, Khaled El Nady, Ibrahim GODA, Kamal berkache. Thank you to all of you. I would also like to thank all the staff and members of the LEMTA laboratory for the welcome, the privileged working conditions that have been offered to me and the good atmosphere they maintain within the group. Special mention and special thank to Mr Rachid RAHOUADJ, for your kindness, your support, your help and your repeated encouragement during the thesis. My thanks go also to my friends who, with this recurring question, "when do you support this thesis? ", Although agonizing in frequent period of doubts, allowed me never to deviate from my final objective. I would like to thank also the Islamic center association for their supports and their financial aids during the thesis, I'm sure they helped a lot of other students to continue their studies. I would also like to thank Mr. Hassan Lakiss, my supervisor at the master, for allowing me to meet Mr. Jean-François GANGHOFFER. I convey my gratitude through these lines. These thanks cannot be ended, without a thought to my parents, my fiancé and my brothers, their encouragement in good and bad days, are for me the founding pillars of what I am and what I do. I love you. And again: Thank you all! ...

Résumé étendu :

Les métamatériaux sont des structures artificielles avec une microstructure optimisée qui ne sont pas trouvés dans la nature, conçus pour avoir des propriétés spéciales qui répondent aux besoins dans l'industrie, en particulier l'aérospatiale, le génie maritime et l'industrie automobile. L'explosion de l'intérêt pour les métamatériaux, tels que les mousses, les matériaux treillis en treillis et le matériau de grille, est due à l'augmentation spectaculaire de la capacité de manipulation par rapport aux propriétés et à la résistance mécaniques, aux capacités d'absorption des ondes sonores, aux propriétés dispersives des ondes mobiles pour l'isolation thermique et acoustique Propriétés.

Le concept des métamatériaux a d'abord été proposé par Veselago pour le champ d'ondes électromagnétiques. Il a prédit que chaque milieu ayant une permittivité négative et une perméabilité négative doit avoir un indice de réfraction négatif. Il est sans doute intéressant de pouvoir concevoir des métamatériaux pour contrôler, diriger et manipuler les ondes sonores, les soi-disant métamatériaux acoustiques. Le contrôle de la propagation des ondes sonores s'effectue principalement par des valeurs négatives des paramètres constitutifs: module de masse β, densité de masse ρ, rapport de poison et également par chiralité. La densité et le module massique sont les analogues des paramètres électromagnétiques, respectivement la permittivité et la perméabilité dans les matériaux à indice négatif.

L'une des caractéristiques les plus intéressantes dans les métamatériaux acoustiques est la génération de phénomènes d'atténuation des ondes: dans certains intervalles de fréquence, les ondes ne peuvent pas se propager à travers les périodiques mais sont atténuées rapidement et complètement au sein des cellules périodiques. La gamme de fréquences correspondante constitue l'intervalle de bande de fréquence.

Par conséquent, de nombreuses activités de recherche ont favorisé la création artificiellement d'écarts de bande en construisant des métamatériaux acoustiques avec des valeurs négatives des paramètres constitutifs, par exemple en ajoutant des résonateurs locaux dans les structures (périodiques): résonateurs à anneaux fendus pour μ négativité, résonateurs Helmholtz Pour un module négatif et un résonateur à membrane pour une masse dynamique négative. L'objectif général de la thèse est de développer des modèles et des outils pour la prédiction du comportement acoustique des structures périodiques et des métamatériaux. Les questions scientifiques examinées dans ce travail comprennent l'impact de la dissipation et de l'amortissement interne, l'analyse des effets d'échelle dus à la microstructure existante et l'impact de grandes déformations sur la propagation des ondes. Le but ultime sera de concevoir et tester des matériaux d'architecture et bio-inspirés avec des propriétés dynamiques et acoustiques spécifiques aux petites longueurs d'onde, en s'appuyant sur des modèles homogènes de continuum.

L'analyse de la propagation des ondes dans les matériaux périodiques 2D et 3D a été réalisée en tenant compte des effets d'échelle microstructurale et du comportement viscoélastique des poutres, en s'appuyant sur une méthode d'homogénéisation spécifique des structures discrètes pour construire un milieu continu de substitution. Deux stratégies différentes pour l'enrichissement du continuum effectif ont été envisagées, soit un milieu micropolaire et un milieu du second gradient. L'analyse dynamique de réseaux fibreux aléatoires a également été faite en 2D. La propagation des ondes linéaires et non linéaires dans le réseau représenté par milieu continu a été analysée en deux parties successives dans la thèse.

Notre contribution aborde deux questions importantes: l'impact de l'enrichissement du continuum efficace sur les courbes de dispersion et la propagation des ondes a été analysé pour les matériaux de réseau 2D et 3D, en considérant le milieu micropolair et le milieu du second gradient. Cette méthodologie a été spécialement appliquée aux méta matériaux montrant un comportement auxétique, qui montrent un gap de bande partielle amélioré à basse fréquence. Nous analysons dans l'impact de l'amortissement des ondes sur les caractéristiques de dispersion des réseaux périodiques, qui sont modélisés sous forme des structures périodiques. La bande de fréquence et le taux d'amortissement sont calculés pour différents réseaux répétitifs 2D et 3D, basés sur la réponse homogène de l'architecture de réseau initialement discrète, modélisée sous la forme des poutres viscoélastiques de Kelvin-Voigt. Trois de ces réseaux (hexagonal réentrant, diamant chiral, hexachiral) sont des métamatériaux auxétiques, puisqu'ils présentent un coefficient de poisson négatif. Les comportement effectives élastique et visqueuse des structures périodiques sont calculées sur la base de la technique d'homogénéisation discrète. Dans le cas 3D, une application sur le textile est présenter en calculant aussi le coefficient de transmission et en étudiant l'influence de temps de retardation sur les courbes de dispersions.

Nous analysons, dans une deuxième étape, la dispersion des ondes élastiques dans des réseaux de faisceaux périodiques basés sur des modèles de second gradient obtenus par l 'homogénéisation du réseau initialement discret, en s'appuyant sur la méthode asymptotique discrète étendue jusqu'au deuxième gradient du déplacement. Les faisceaux de réseau ont un comportement viscoélastique décrit par le modèle de Kelvin-Voigt et le modèle de viscoélasticité du second gradient homogénéisé reflète à la fois la topologie du réseau initial, l'anisotropie et les caractéristiques microstructurales en termes de ses paramètres géométriques et micromécaniques.

Les modèles de continuum enrichis avec les gradients d'ordre supérieur du déplacement et de la vitesse introduisent des paramètres caractéristiques de longueur qui tiennent compte des effets de microstructure au niveau mésoscopique. Une étude comparative des relations de dispersion et des évolutions du rapport d'amortissement des ondes longitudinales et de cisaillement a été réalisée pour quatre réseaux (le réseau de diamant chiral, les réseaux classiques et réentrants et le pantographe). Le modèle développé permet d'analyser à la fois les effets de l'amortissement et l'échelle de longueur interne à travers les gradients de déplacement du second ordre sur les caractéristiques de propagation des ondes. On observe une augmentation importante de la fréquence naturelle due aux effets de second ordre. Pour le réseau pantographique, la vitesse de phase pour les modes longitudinal et de cisaillement est identique et n'est pas influencée par la direction de propagation de l'onde. Les résultats obtenus montrent globalement que le réseau du pantographe présente les meilleures caractéristiques acoustiques.

L'analyse dynamique des milieux fibreux aléatoires a été présenter dans cette contribution. Afin de contourner la complexité de l 'exécution de calculs dynamiques à l' échelle microscopique du réseau fibreux aléatoire, nous développons et identifions des modèles de "Couple stress" et du milieu du seconde gradient comme des milieux continus efficace au niveau mésoscopique des fenêtres d 'analyse afin d' analyser les effets de taille de ces réseaux sur leurs propriétés dynamiques. Les propriétés mécaniques statiques qui sont à la base de l'analyse dynamique sont calculées grâce à des simulations de EF effectuées sur des fenêtres d'analyse soumises à des conditions de limites mixtes permettant de capturer les modules efficaces classiques et non classiques. Les propriétés acoustiques sont capturées par les diagrammes de dispersion et les tracés des phases; Nous analysons l'influence sur les propriétés dynamiques de trois grandes quantités d'intérêt, à savoir la longueur de flexion de la fibre, la taille de la fenêtre d'analyse et la densité des fibres.

La deuxième question importante est l'analyse de la propagation d'ondes non linéaires dans des réseaux répétitifs soumis à des changements de configuration sous grandes charges cinématiques appliquées et reflétant des non-linéarités géométriques dans certains cas et dans des autres la non linéarité matériel du au fragilité des matériaux. Dans la seconde partie du travail, un schéma incrémental pour la mise à jour de la fréquence et de la vitesse de phase en fonction de la déformation appliquer au structure a été développé, en considérant successivement des matériaux de réseau répétitif 1D, 2D et 3D basant sur la méthode de Linstedt-Poincaré. Il incorpore une mise à jour de la fréquence et de la vitesse de phase des ondes de propagation en fonction de la densité effective et de l'état de déformation finie du continuum efficace utilisé comme support de substitution pour le réseau répétitif initial. La déformation appliquée est montrée avoir des effets significatifs sur la fréquence d'onde et la vitesse de phase. L'influence de la densité effective sur la relation de dispersion et les diagrammes de bandes sous l'application d'une déformation incrémentielle sur la cellule d'unité de réseau est représentée. Un aspect original préconisé dans le présent travail est la dérivation de la méthode de perturbation dans un contexte 2D pour des structures périodiques non linéaires couvrant des non-linéarités géométriques et matérielles; En basse fréquence, la méthode est basée sur les propriétés effectives du milieu continu de substitution, en utilisant la méthode d'homogénéisation. En plus haute fréquence, une analyse sur la structure périodique est menée en utilisant la méthode de perturbation. Une comparaison entre les deux méthodes permet de valider la méthode d'homogénéisation en basse fréquence.

Nous analysons la propagation d'ondes non linéaires dans des structures périodiques discrètes, en considérant successivement des situations 1D, 2D et 3D. L'analyse des ondes est réalisée sur la base de la construction de la densité d'énergie de déformation effective des structures périodiques dans le régime non linéaire, dont la cinématique intègre la déformation de Green-Lagrange. Le milieu continu non linéaire obtenu est du second gradient et il présente deux modes de propagation : un mode subsonique évanescent qui disparaît au-delà d'un certain nombre d'onde critique et un mode supersonique caractérisé par une augmentation de la fréquence avec le nombre d'onde. Dans le cas de faible non-linéarité, un mode supersonique se produit et les courbes de dispersion se situent au-dessus de la courbe de dispersion linéaire. Pour une nonlinéarité plus élevée, l'onde passe d'un mode supersonique à un mode subsonique évanescent et l'onde ne passe plus pour certaines valeurs du nombre d'onde. Une diminution importante de la fréquence se produit pour les modes subsonique et supersonique lorsque la structure devient auxétique, les modes longitudinaux et transversaux étant très proches. L'influence des paramètres géométriques du réseau sur les relations de dispersion est analysée.

Nous avons également construit des milieux hyperélastiques par des tests virtuels reposant sur la méthode d'homogénéisation à la base de cette thèse et nous avons obtenu différents types d'équations de propagation d'onde (équations de Burgers et de Boussinesq) et de solutions (ondes de choc, ondes solaires et ondes harmoniques planaires) dont les propriétés acoustiques ont été étudiées.

Nous analysons les propriétés acoustiques des faisceaux microstructurés comprenant un matériau de réseau répétitif subissant des changements de configuration conduisant à des non-linéarités géométriques. La loi constitutive effective est évaluée successivement comme un continuum non linéaire de grade 1D efficace de premier et second ordre, basé sur un schéma incrémentiel entraîné par la contrainte écrit sur la cellule de l'unité de référence, en tenant compte des changements de la géométrie du réseau. Les équations dynamiques du mouvement sont ensuite écrites, conduisant à des relations de dispersion spécifiques. L'équation de Burgers non visqueux est obtenue comme équation de propagation d'onde spécifique pour le continuum effectif du premier ordre, alors que l'équation de Boussinesq est obtenue en considérant un continuum efficace de gradient de déformation. La présence d'un terme d'ordre du second gradient dans l'équation non linéaire du mouvement conduit à la présence de deux modes différents: un mode subsonique évanescent pour une non-linéarité élevée qui disparaît au-delà de certaines valeurs de nombre d'onde et un mode supersonique pour une non-linéarité faible. Cette méthodologie est appliquée pour analyser la propagation des ondes au sein de différentes microstructures, y compris les hexagones réguliers et réentrants et le motif textile à armure simple.

La recherche présentée dans cette thèse apporte les contributions originales suivantes :

• Une méthodologie générale pour analyser les effets de dispersion dans les réseaux de poutres périodiques comportant des effets de microstructure, basée sur la construction d'un milieu homogénéisé du second gradient ;

• Analyse de la propagation des ondes dans les milieux architecturés dissipatifs, prenant en compte également les effets de microstructure;

• Calcul des caractéristiques de propagation des ondes dans des milieux fibreux aléatoires basés sur des milieux effectifs du second gradient et des milieux effectifs micropolaires ;

• Analyse de la propagation des ondes des milieux continus non linéaire construit par l'homogénéisation de la réponse de réseaux répétitifs non linéaires;

• Outils d'analyse numérique efficaces pour évaluer les effets de dispersion non linéaire dans les réseaux de poutres périodiques ;

• Méthode générale de perturbation pour les non-linéarités géométriques pour l'analyse de la propagation des ondes à des fréquences plus élevées dans des réseaux non linéaires. Metamaterials are artificial structures with an optimized microstructure that are not found in nature, designed to have special properties that meets the needs in industry, especially aerospace, marine engineering and automotive industry. The explosion of the interest in metamaterials, such as foams, lattice truss materials and grid material, is due to the dramatically increased manipulation ability over mechanical properties and strength, absorption capacities of sound waves, the dispersive properties of traveling wave for thermal and acoustics insulation properties.

List of Figures

The concept of metamaterials was first proposed by Veselago [1] for the electromagnetic wave field. He predicted that each medium with negative permittivity and negative permeability shall have a negative refractive index. It is undoubtedly of interest to be able to design metamaterials to control, direct, and manipulate sound waves, the so called acoustic metamaterials. Control of the propagation of sound waves is mostly accomplished through negative values of constitutive parameters: bulk modulus , mass density ρ, poison ratio and also via chirality. The density and bulk modulus are the analogs of electromagnetic parameters, respectively permittivity and permeability in negative index materials. [2] presents a thorough history of acoustic metamaterials. In recent years, the synthesis of novel acoustic metamaterials with negative Poisson's ratio (defined as the negative of the ratio between transverse and longitudinal strains in uniaxial elastic loading) has been developed. In contrast to conventional materials, these so-called "auxetic" metamaterials contract in the transverse directions when compressed uniaxially. The auxetic behavior of the structure has significant effect on the absorption energy and also on the band gap.

It turns out that the constitutive parameters like effective density and bulk modulus can be adjusted by modifying the microstructure in man-made materials in order to get partial or full band gaps; there also in nature periodic systems that have great ability to adjust and enhance their performances and specially the frequency spectrum. The microstructures inspired from nature but which are then modified by man to produce repetitive networks has a great potential in acoustics.

Examples of such natural and artificial periodic networks include fibrous reinforcements, repetitive structures used in civil engineering, polymeric foams, trabecular bone, and the filamentary network of biological membranes, as illustrated on Fig. 1.1. A specific naturallyoccurring periodic system is the hexagonal honeycomb structure constructed by various bee species. Honeycomb is one of the most important structures in periodic composites for its ability to form complete band gaps [3]. At larger length scales, periodic structures that respond to entirely different frequency ranges appear in many engineering structures (railroads, aircraft structures, buildings, bridges). Only within the last forty years have researchers began to explore the unique system dynamics of these materials in depth [4]. It turns out that the ability to completely reflect incident waves at some frequencies (band gaps) while allowing others to propagate (pass bands) is a universal feature among periodic materials Fig. 1.2.

Motivation of the work

The general objective of the thesis is to develop models and tools for the prediction of the acoustic behavior of periodic structures and metamaterials. The scientific issues considered in this work include the impact of dissipation and internal damping, analysis of scale effects due to the existing microstructure, and the impact of large deformations on wave propagation. The ultimate goal will be to design and test architectured and bio-inspired materials with specific dynamic and acoustic properties at small wavelengths, relying on homogenized continuum models.

Effect of damping

The existing analyses in the literature of the propagation of elastic waves within periodic structures and metamaterials mostly deal with non-damped systems, and few works treat damped structures [5][6][7]. The analysis of the dispersion relations for dissipative periodic lattices having the attributes of metamaterials is an actual research topic raising the interest of researchers in the field of wave propagation phenomena. The effect of damping on the dispersion features is very important, as shown in recent studies devoted to composite materials [8][9][10][11].

Studies of the propagation of elastic waves in damped media started more than 40 years ago with an analysis of damping in infinite periodic structures by Mead [12]. From this early period, studies considering partial aspects of the effect of viscous damping on the band structure appeared in the specialized literature; see [8] and references therein, without however providing a detailed analysis of the effects of damping on the dispersion band structure.

In a second stage, the impact of damping on the band structure of infinite periodic crystals was studied, considering either free wave propagation (for impulsive loading), or time harmonic wave propagation (in case of a forced harmonic loading); an overview of the works considering these two situations is presented in the recent contribution of Wang et al. [13]. For free wave propagation, the wavenumber is a real number indicative of propagating models without taking into consideration the geometrical attenuation, while the frequency is a complex number, in which the imaginary part measures the attenuation in time due to the presence of damping.

In [13], the authors analyze the propagation of elastic waves in acoustic metamaterials, and show the dispersive and dissipative effect of viscosity on the complex band structure and transmission spectra.

We analyze as a novel aspect the impact of damping on the dispersion features of periodic dissipative lattices, taking into addition the impact of an auxetic behavior on wave propagation.

The band diagram structure and damping ratio shall be evaluated for different repetitive lattices, based on the homogenized continuum response of the initially discrete lattice architecture, modeled as Kelvin-Voigt viscoelastic beam elements.

The impact of the microstructure on the dispersion relation will be accounted for by enhancing Cauchy continuum in two different ways. A first strategy is the consideration of additional rotational degrees of freedom in the context of an effective micropolar continuum obtained by homogenization of viscoelastic periodic networks, a second strategy by incorporating additional higher-order gradients of the displacements in the context of second gradient continuum.

Microstructural effects adopting higher order media

Scale effects due to the presence of a microstructure play an important role in the design and acoustic properties of metamaterials. Many advanced theories and models have been proposed to study wave propagation problems accounting for non-locality and microstructural effects in materials. One category of approaches is multiscale homogenization techniques, which compute asymptotically the solutions of the wave equations involving multiple spatial and temporal scales, to capture the long-term response of the homogenized response. A lot of works in this direction include contributions of Mindlin [14][15][16], Aifantis [17][18][19][20] in connection with the higher-order strain theories, which enhance the classical continuum theories by considering higher order gradients of either the displacement or strain fields, accompanied by internal length scales [21,22]. Different homogenization techniques have been proposed to build continuum descriptions of discrete materials [23,24]; an extended list of references and an overview on gradient elasticity theories applied in dynamics can be found in recent works [25][26][27][28][29][30]. Following another strategy, micropolar theories incorporate additional rotational degrees of freedom, and have been developed by Eringen and co-workers [31,32], a subclass of models being the couplestress theories [33][34][35][36]. No classical theories have proven successful in dynamic problems of wave propagation in beams and half-spaces [37][38] by removing singular behaviors (like discontinuities or singularities of the fields) inherent to classical theories and showing their ability to capture the expected size and wave dispersion effects specific to the simultaneous existence of several lengths scales.

Microstructural effects adopting higher order displacement gradients

The presence of homogenized models based on classical Cauchy-type elasticity theory is able to provide realistic predictions of many effects arising from small scales, amongst of them wave dispersion. Classical theories based upon the sole first order displacement gradient lack indeed internal length parameters, characteristic of the underlying microstructure. This explains the success of gradient-enriched theories in capturing microstructural effect on the macroscopic behavior of materials, by including high-order gradients associated to internal lengths representative of the microstructure.

Gradient elasticity theories constitute an extension of the classical equations of elasticity by incorporating additional higher-order gradients of the displacements. Many theories based on an enrichment of the classical elasticity framework by higher order gradients have been proposed in the past to overcome deficiencies of classical elasticity [39][40] and plasticity theories [41] in describing both static and dynamic phenomena, including size effects, strain and stress fields in the neighborhood of singularities [42][43] and wave dispersion in dynamics [44][45].

Much research has been devoted to model static and dynamic problems arising from the various developed gradient theories; all obtained results showed an enhanced stiffness concomitant to the increase of the gradient coefficients, see e.g. [46] and references therein. Note that different homogenization techniques have been proposed to build gradient elasticity continuum models of discrete materials [47].

Generalized continuum theories have been shown to offer an attractive alternative for capturing dynamic behaviors overlooked by classical elasticity, especially dispersion relations.

Applications of gradient elasticity in dynamics have fostered extensive research [48][49][50]. In [51],

the dynamic behavior of periodic lattice materials is investigated using an equivalent higher-order continuum model obtained by the homogenization of the equations of motion. Considering dynamic aspects and especially wave propagation phenomena, the impact of the microstructure in heterogeneous materials on the dispersive propagation of elastic waves was first recognized and the analyzed in seminal paper [52], who proposed several non-local continuum models to capture the dispersion relation of planar waves. It has however been proven by experiments that most waves are dispersive, that is, each wavenumber travels with a different phase velocity. This explains the success of gradient-enriched theories in capturing dynamic behaviors overlooked by classical elasticity.

A detailed comparison between the dispersive characteristics of various simplified models of gradient elasticity can be found in [53].

The analysis of wave propagation in periodic beam networks made of viscoelastic beams of type Kelvin-Voigt enriched by second order gradient models, obtained by the homogenization of the initially discrete network shall be exposed, relying on the discrete asymptotic method extended up to the second gradient of the displacement.

Wave propagation in random fibrous media

As mentioned in very recent contributions [54], the dynamics of periodic fiber networks has recently raised a lot of interest, whereas the study of wave propagation in random fibrous networks has not deserved yet the attention of researchers. Although fibers networks are often subjected to dynamical loadings like vibrations, most of the works have indeed been devoted to the analysis of their static behavior. We will use generalized continuum theories at an intermediate mesoscopic level in order to address the issue of size effects related to the random fibrous microstructure; those theories have been extensively used to explain size effects for a wide class of materials, but not for random fibrous networks to our knowledge. This constitutes one original aspect advocated in the present work. Real materials such as biological membranes and tissues often exhibit a number of important length scales, which must be included in any realistic model.

Effect of the structure nonlinearity on the dispersion curve

The propagation of elastic waves in a linear framework has deserved a lot of works in the literature, whereas only a few authors analyzed so far wave propagation in nonlinear media. The incorporation of nonlinear aspects of wave propagation in structures is necessary whenever large deformations occur [55][56][57][58], but it remains a considerable challenge nowadays. Two types of nonlinearities may be present in a broad sense, which can be classified as material nonlinearities and geometrical nonlinearities [59]. The propagation of nonlinear waves in periodic structures is accompanied by a number of new phenomena that are different and can never be observed for linear media. The presence of a nonlinearity in periodic structures results in the dependency of the wave propagation, phase and group velocities upon the amplitude of the wave; this phenomenon deserves the name amplitude dependent dispersion relation. It opens new possibilities for a passive tuning of the dispersion band structure through an amplitudedependency of propagating waves, thereby going beyond a mere control of the dynamic and acoustic properties of repetitive structures by the design [60] or by application of an external stimulus.

This entails that solutions of the wave propagation equations are much more complex compared to harmonic plane solutions of the linear wave equations: nonlinear periodic structures support a variety of wave solutions depending on wave amplitude, waves interactions, and the type of nonlinearity; one can for instance mention solitary wave solutions for Boussinesq type equation, and shallow water waves for Burger's equation [START_REF] Daraio | Tunability of solitary wave properties in onedimensional strongly nonlinear phononic crystals[END_REF][START_REF] Sen | Solitary waves in the granular chain[END_REF][START_REF] Vakakis | Forced localization in a periodic chain of non-linear oscillators[END_REF][START_REF] Vakakis | Nonlinear wave transmission in a monocoupled elastic periodic system[END_REF][START_REF] Vakakis | Non-linear normal modes (nnms) and their applications in vibration theory: an overview[END_REF][START_REF] Vakakis | Resonant oscillations of a weakly coupled, nonlinear layered system[END_REF][START_REF] Starosvetsky | Traveling waves and localized modes in one-dimensional homogeneous granular chains with no pre-compression[END_REF].

The effect of a pre-stress or pre-strain on wave propagation in homogeneous anisotropic media has raised the interest of many authors but it is still an important challenge to study the effect of preexisting finite elastic deformations on wave propagation. The initial deformation must be large enough to change the geometry of the medium, since an infinitesimal initial deformation would not affect the properties of the material based on the superposition principle valid for small deformations. The incremental effective properties of pre-stressed homogeneous media undergoing large deformation have been analyzed in [START_REF] Bhatnagar | Nonlinear Waves in One-Dimensional Dispersive Systems[END_REF][START_REF] Andrianov | Dynamic homogenization and wave propagation in a nonlinear 1D composite material[END_REF], wherein the authors put some restrictions on strain energy sufficiently enough to allow elastic waves to propagate within the material.

We analyze in this thesis the propagation of nonlinear waves in homogenized periodic nonlinear networks, considering successively 1D, 2D and 3D situations based on a homogenized continuum of the initial repetitive network. The obtained nonlinear effective constitutive law of Cauchy type shall be enriched by second order gradient terms to account for microstructural effects. This will be the basis for analyzing the impact of large deformations on wave propagation. An incremental scheme for the update of the frequency and phase velocity of the homogenized medium has been developed to compute modifications of the existing band gaps under applied deformation, considering successively 1D, 2D and 3D repetitive network materials.

Organization of the thesis

The thesis is divided into two parts, the first part dealing with wave propagation aspects in a small strains context, which is then extended to the nonlinear setting. This document is organized as follows.

In Chapter 2, we compute the dispersion relations for dissipative 2D periodic lattices having the attributes of metamaterials modeled as viscoelastic beam-lattices taking into account the microstructural effects by higher order medium. In Chapter 3, we extend the previous analyses to 3D viscoelastic periodic structures, and study the impact of wave damping on the dispersion features. Chapter 4 deals with the computation of the dispersion of elastic waves in periodic beam networks based on second order gradient models obtained by the homogenization of the microstructure. The dynamical analysis of random fibrous networks modeled on the basis of couple stress and second gradient effective continua is exposed in Chapter 5.

The second part of the thesis deals with the nonlinear effects on wave propagation.

The influence of large pre-strain on the propagation of acoustic waves is detailed in Chapter 6, whereby an incremental scheme for the update of frequency and phase velocity of the computed homogenized medium is developed for 1D, 2D and 3D structures successively. In chapter The dynamical equations of motion are obtained for an equivalent homogenized micropolar continuum evaluated based on the homogenized properties, and the dispersion relation and damping ratio are obtained by inserting an harmonic plane waves Ansatz into these equations.

The comparison of the acoustic properties obtained in the low frequency range for the four considered lattices shows that auxetic lattices attenuate waves at lower frequencies compared to the classical hexagonal lattice. The diamond chiral lattice shows the best attenuation properties of harmonic waves over the entire Brillouin zone, and the hexachiral lattice presents better acoustic properties than the reentrant hexagonal lattice. The range of validity of the effective continuum obtained by the discrete homogenization has been assessed by comparing the frequency band structure of this continuum with that obtained by a Floquet-Bloch analysis.

Overview

Architectured materials, and especially repetitive network materials made of structural elements like beams, constitute a wide class of structures having the capacity to filter waves in certain directions and frequency range. The dynamic response and wave propagation properties of periodic lattices and structures have raised numerous studies especially in aeronautics, for the objective of reducing or absorbing vibrations, shock and sound in structural components [1,2].

The mechanical response of such networks has fostered a lot of research activity in the literature, but the evaluation of their dynamical and acoustic properties, especially in the high frequency domain, remains a scientific challenge [3][4][5][6][7][8][9].Materials [10,11] structures [1,2], and devices [12] exploiting spatial periodicity are involved in a growing number of areas, such as ultralight architectured materials [10,11,13], phononic crystals [14][15][16][17][18][19][20][21][22], [8,5] or acoustic metamaterials [3][4][5][6][7] and [23][24][25]. These structures raised in the recent years a considerable interest [26], most of the research works concentrating on the conception of materials and structures having a periodic microstructure in order to get complete sound propagation in a certain frequency range, called the spectral band gap (a frequency band gap in which waves are forbidden irrespective of propagation within the structure). Many applications of the propagation of elastic waves deal with non-damped systems and structures [27,28] but up to now few works treat damped structures [29][30][31]. The computation of the dispersion relations for dissipative periodic lattices having the attributes of metamaterials is an actual research topic raising the interest of researchers in the field of acoustics and wave propagation phenomena. Damping has an important influence on the dispersion features, as shown in recent studies devoted to phononic composite materials [32][33][34][35]. The nature of damping (proportional or not in the case of viscous damping) and its intensity can lead to specific phenomena, such as branch inversion (high frequency branches take over low frequency branches), branch switching on or off translate the fact that certain dispersive branches do not cover the integrality of the first Brillouin zone, and completely or partially forbidden branches in certain domains of variation of the wave number [32]. An extension of the method initially proposed by Rayleigh [36,37] has been done in [32] for the analysis of Bloch dispersive damped waves in phononic materials in the presence of viscous damping, coined by the same previous authors Bloch-Rayleigh perturbation method. One of the advantages of the method lies in the fact it relies on the standard non damped dispersion relations, without increasing the order of the eigenvalues problem. Rayleigh method has been used for a wide class of non-viscous models [37].

We analyze as a novel aspect in this chapter the impact of wave damping on the dispersion features of periodic dissipative lattices, considering especially the impact of an auxetic behavior on wave propagation. The band diagram structure and damping ratio shall be computed for different repetitive lattices, based on the homogenized continuum response of the initially discrete lattice architecture, modeled as Kelvin-Voigt viscoelastic beams. We shall show that auxetic lattices attenuate waves at lower frequencies compared to the classical non auxetic lattices. The impact of the microstructure inherent to the unit cell is assessed through the consideration of additional rotational degrees of freedom of the effective micropolar substitution continuum.

Scale effects due to the presence of a microstructure play an important role in the design and acoustic properties of metamaterials [38]. The presence of a microstructure gives rise to specific mechanical properties and behaviors that cannot be reached with homogeneous materials, like negative Poisson's ratio, associated to a class of materials called auxetic, due to the occurrence of one or several negative Poisson's ratio [39][40][41][42][43][44], their structure being characterized by a periodical network of four or six elastic ligaments. Lakes was the first in [45] to propose such a microstructure, in reference to the case of planar isotropic hexachiral networks, consisting of a periodic array of rigid rings connected by softer elastic ligaments. Later on, [45] have shown that this material has a Poisson coefficient equal to -1, under the assumption that one can neglect the axial deformation of the ligaments. Starting from the work of Lakes devoted to the hexachiral lattice, a series of contributions from the literature was further devoted to the computation of the effective mechanical properties of the hexachiral lattice and variants of it (like the tetrachiral lattice) [38] and [45][46][47]. The effect of the localized resonant units in these two specific lattices can be associated to localized rotational modes at the level of their unit cell. The discrete rotations within the lattice lead to rotational degrees of freedom at the scale of the homogenized continuum, corresponding to an effective micropolar model. The computation of the homogenized effective static properties of such period lattices is the first step towards the characterization of their band structure.

These metamaterials can lead to a nearly complete sound attenuation over a certain frequency range, coined the acoustic wave spectral gap or band gap [38]. Acoustic metamaterials may exhibit band gaps at low frequencies associated to resonant elements; in [48], the authors derive the dispersion relations and the damping factor for a 1D damped mass-chain. This goes in line with the development of sonic crystals presenting such spectral gaps, based on the idea of internally resonant structural units leading to negative effective elastic constants [47,49]. Sonic crystals exhibiting stopband have been developed in [26,7,47,49]. More recently, [23] proposed a periodic metamaterial with resonant internal units, allowing adjustable band gaps at low frequencies; this effect is associated to rotational localized modes of deformation obtained by the presence of a chiral like microstructure of the periodic unit cell. In [50], the authors recently obtained a band gap structure in periodic tetrachiral material.

This chapter is organized as follows: we compute in section 2 as a first step the homogenized viscous behavior of auxetic lattices giving rise to a micropolar viscoelastic continuum model, based on an extension of the discrete homogenization technique for beam-lattices, developed in the pure elastic case in [51]. The extension to a viscous behavior of the beam elements of the dedicated homogenization method proposed in the present work constitutes a novel aspect in comparison to the existing literature. We next introduce the effective constitutive laws into the dynamical planar equilibrium equations, and formulate an eigenvalue problem on view of the computation of the dispersion relation, relying on the generalized plane wave propagation equations (section 3). The analysis of free wave motion in dissipative micropolar models at low frequencies is performed in terms of the dispersion relations and damping factor in section 4. The evolution of the phase velocity with frequency together with the evolution of the modulus of the phase velocity versus the wavenumber for the effective viscous medium is analyzed in section 5.

The deformation of the hexagonal lattice unit cell versus time is determined in terms of the Eigen waves in section 5. We assess the range of validity of the homogenized properties by comparing the dispersion curves to those obtained by Bloch's theorem (section 7). Finally, we conclude in section 8 by a summary of the main results and a few perspectives for future developments.

Homogenized viscoelastic micropolar behavior of periodic beam lattices

A summary of the discrete asymptotic homogenization method has been presented in [51]; it proceeds in two steps. In a first step, the equilibrium equations for an Euler-Bernoulli beam are written at the nodes; a second homogenization step leads to the continuous effective model. The discrete homogenization method consists in assuming asymptotic series expansions of node displacements, tension, moments and external forces versus a small parameter labeled ε, defined as the ratio of a characteristic length of the basic cell to a characteristic length of the lattice structure. Those expansions are then inserted into the equilibrium equation, conveniently expressed in weak form. The balance equation of the nodes, the force-displacement relation and the moment-rotation relations of the beams are developed by inserting those series expansions and using Taylor's expansion of finite differences. The discrete sums are finally converted in the limit of a continuous density of beams into Riemann integrals, thereby highlighting continuous stress, deformation and velocity of deformation measures.

The kinematic and static variables for any beam in the lattice are represented in the local coordinate system associated to the viscoelastic beam with Bernoulli-Euler kinematics (Fig. 2.1).

Fig. 2.1 Kinematic and static beam variables.

2.2.1Expressions of forces

A few words regarding notations are in order. Vectors and second order tensors are denoted with boldface symbols; the tensor product of two vectors is denoted and the inner product of two tensors , AB is denoted by the dot product .

AB.

The forces at the extremity of the beam can be related to the kinematic variables in the global coordinate system as follows

             p s s F K T U C T U (2.1)
Where  K is the stiffness matrix,   C the viscosity matrix (evaluated in the Appendix A) and

 

T the local to global transformation matrix, consisting of two 3x3 rotation matrices R.

Therefore, based on the framework of viscoelastic Bernoulli beam theory, the efforts and moments can be expressed as  is the extensional viscosity. Here and in the sequel, the dot over a variable denotes a time derivative.

        A E O E O E E L L A e N      b b u u u u e e (2.2)       .. 3 3 12 L 12 L 2 2 E E E E O O EO O e EI I T L L                               b b uu e u e u (2.3)         .. 2 2 2 2 2 3 2 3 O OE z O E O E O E EI I M L L L L e                         b b e u u e u u (2.4)         .. 2 2 2 2 2 3 2 3 E EO E O E O z E O EI e I M L L L L                            b b e u u
The material coefficient E is the Young modulus of the beam material; parameters A and L are the cross-sectional area and length of the beam respectively; I is the quadratic moment of the beam. For simplicity reasons, a rectangular section of the beam is considered, with a constant unit thickness e =1 and width t. Hence, the cross-sectional area is equal to 2 At  and the quadratic moment of the beam is evaluated in bending as

3 12 t I 
, and vectors , bb  ee are respectively the unit director and transverse unit vector for each beam.

Main steps of the discrete homogenization method

The main steps of the discrete homogenization method involved to compute the effective anisotropic micropolar 2D viscoelastic continuum are summarized in the sequel in algorithmic format.

1. Asymptotic expansion of geometrical and kinematic variables for each beam b, in curvilinear coordinates denoted  in the sequel:

-The beam length l εb =εL b -The beam width t εb =εt b -The nodal displacement:

        2 0 1 2    O O O          u u u u (2.6)         2 0 1 2 0 0 1 2 2 () ( ) ( ) ( )             E ib ib ib ib E E E E i                           u u u u u u u u (2.7)
This leads to the relative displacement

              ( ) ( ) ( ) ( ) b 2 0 O O 1 1 2 1 2 2 U                     uu U U b b b b b b ib E E b b E O i                  u u u u u (2.8)
ib  is the shift factor, which belongs to the set   1, 0,1  ; it traduces the fact that the end node belongs to the next neighboring cell (Fig. 2.2).

-The relative velocity can be obtained similarly as,

              ( ) ( ) ( ) ( ) b2 2 1 0 OO 1 1 2 2 U                 u u u u u u u b b b b b EE b ib EO i b b                  U U (2.9)
-The expansion of the nodal microrotation writes (limited to the first order in )

( ) ( ) ( ) ( ) 01 0 0 1 bb b b OO b O E E ib i                   
(2.10)

The rotation velocity can be obtained similarly.

Write the asymptotic expansion of forces ,

bb NT  and moments

( ) ( ) , O b E b z z M M  (exposed in Appendix B).
3. Write the equilibrium of forces and moments in virtual power form

      3 . ( ) ( ) 0     vv i R b b Β O b E b     F (2.11)       3 ( ) ( ) ( ) ( ) 0      M w M w i R O b E b b Β O b E b      (2.12)
with  v therein a virtual velocity field chosen to vanish on the edges and with w  the virtual rotation rate. The force vector F b decomposes into a normal and a transverse contribution as

      i ( ) ( )     v v v ib O b E b      (2.14) 0 ( ( )) ( )  w w Ob     0 () 0 ( ) ( )      w ww E b i ib i         (2.
      i i vZ F v F x dx     i . 0      v S i d    (2.16)   i . 0      μ i d    w (2.17)
with i S the stress vector therein decomposed into (2.18) and with i μ the couple stress vector expressing at two orders (2.19) 6. Constitutive law: the stress and couple stress tensors express versus the kinematic variables as

    1 2 1 1 2 2 sec i i i b b b b ib b b b b ib b B b B R R first power of ond power of N T N T                   S S S e e e e
    2 2 1 2 1 1 2 2 1 1 2 2 R i i i E O ib E O ib z z z z z z b B b B R M M M M                             μ μ μ e e
                    i i i i i 1 2 1 2 i i i i i i 2 i i 2 i 1 2 1 2 i i i i 1 1 1 1 1 1 1 1                                                                       R R R R S S S S S R R R R m μ μ μ μ μ e v e v e v e v EE BB BB K K g g g g g g g g          (2.20)
Previous relations can therefore be written as

                              e e v v viscous part elastic part E B E B                                  e e v v viscous part elastic part B K B K m  (2.21)
With g the Jacobian of the transformation from Cartesian to curvilinear coordinates and R the position vector for any point.

                                                                                                     (2.22)
This constitutive law can be further simplified basing on symmetry properties of the studied lattices: it has indeed been shown that for centro-symmetrical lattices the pseudo-tensors [B e ] and

[C e ] vanish (similarly, the tensors [B f ] and [C f ] vanish) [51]. A solid which is isotropic with respect to coordinate rotations but not with respect to inversions is called noncentro-symmetric (it does not have a center of symmetry). The centrosymmetric nature of the lattice leads to an important simplification of the stress and couple stress vectors, which write

  2 2 2 1 2 i i E O ib z z z bB R MM             μμ e (2.24) with 1 1 , b b N Tand, 2 b M
, respectively, the first order normal and transverse effort and the second order moment.

In this study, four lattices are selected for their expected interesting acoustic properties: the hexagonal lattice, the re-entrant lattice (in fact the hexagonal with a negative angle), the diamond chiral lattice and the hexachiral lattice, pictured in Fig. 2.3. The mechanical and geometrical parameters for the four studied lattices are given in Table 1. Hexachiral and tetrachiral beam networks presenting local resonators at the nodes have been studied in [38], relying on micropolar homogenization. The authors compare the results of micropolar and second order gradient homogenization analyses, in order to evaluate the domain of validity of the beam network model. These structures with a periodic network give rise to an auxetic chiral behavior and present band gaps, as shown in [47,49]. The dynamical equations are those of an effective Cosserat continuum and the elastic moduli and inertia terms have been obtained for hexachiral and tetrachiral lattices. The effective constitutive equations of the network are then used to study the propagation of planar waves [38] along the directions of the lines connecting the centers of the network. The static and dynamic properties (frequency band gaps) of anti-tetrachiral cellular solids are analyzed in [41]. The homogenization of periodic hexachiral and tetrachiral honeycombs is performed based on both micropolar homogenization and a second gradient homogenization scheme in [42]. An analytical model for the computation of the effective classical properties of 3D reentrant auxetic materials is developed in [49]. The static behavior and wave propagation analysis of rectangular chiral lattices is analyzed in [50],

based on a 2D orthotropic chiral micropolar theory. Other studies devoted to the homogenization of auxetic hexachiral materials have been done in [46]. In [47], the authors proposed a metacomposite incorporating a 2D hexachiral network with resonant elastic inclusions, in order to obtain band gaps at low frequencies. This meta-composite has been analyzed via a numerical model in which the ligaments have been modeled by beam like elements and the inclusions by a 2D finite element model.

Dynamical equilibrium and characteristic equation

u v m m I                (2.25)
The effective density therein is given in general by , with J 1 the inertia of the cell [38]; these quantities are next computed specifically for all four lattices.

In the specific case of the hexagonal lattice parameterized by the angle  (for re-entrant lattices, it holds  <0) (Fig. 2.3 a, b), one obtains

    2 cell A 2L cos (1 sin )    (2.26)
together with the following inertia parameters:

        * * 3 , 2Lcos
(1 sin ) 2cos (1 sin )

t tL I            (2.27)
For the hexachiral lattice, one obtains

            (2.31) The curvatures 1 , 1 2 ,2 , к к 
   (2.32)
define two additional kinematic quantities that have to be considered together with the constitutive equations.

We recall accounting for the assumption of centrosymmetric unit cells that the general constitutive equations of lattices made of viscoelastic Kelvin-Voigt beams is written as

               
                                                                                          (2.33)
The matrices , , ,

e e v v
A D A D are dependent on the specific considered lattice.

For a harmonic wave propagating along an axis in an infinite planar micro polar medium, the generalized displacement field with components , , UV at a point r takes the following form [58]: The eigenvalue problem for Eq. (2.35) yields a characteristic equation of the form:

6 5 4 3 2 0. a b c d e f              (2.36)
The roots of Eq. ( 36) may be expressed as

               (2.38)
Relying on these expressions, we can plot the dispersion curve (a plot of frequency against wave number) for the dissipated frequency and the damping ratio versus the wave vector on the edges of the irreducible part of the first Brillouin zone (O, A, B), which is defined as the smallest area allowing a full representation of   k  , or by changing the direction of wave propagation and the wavenumber. The first Brillouin zone for the hexagonal, re-entrant and hexachiral lattices is described in details in [52,53], and in Appendix D for the diamond lattice. Note that for a viscous damped system, the natural frequency ( ) ns k  is equal to the un-damped natural frequency.

Dispersion relations and damping ratio for the selected lattices

The frequency and damping ratio dispersion curves can be obtained for any lattice based on the previous methodology after calculating the homogenized micropolar continuum. In this section, we shall illustrate this aspect for the classical and re-entrant lattices, the hexachiral and diamond chiral lattices, and especially illustrate their higher capacity to absorb and dissipate energy, which could be useful for packaging materials, to absorb shocks or to improve the acoustic damping

properties. An extensive review of the performance and applications of auxetic materials can be found in [47]. A preliminary work devoted to the computation of the effective static mechanical properties of auxetic lattices has been done in [54]. The abbreviations H, R, D and H-C stand for the wordings hexagonal, Re-entrant, Diamond Chiral and Hexachiral respectively.

We first evaluate the frequency and damping ratio dispersion curves of the hexagonal and reentrant lattices. The first set of results is shown in Fig. 2.4: the snapshots therein display the frequency band structure and damping ratio for each of the modes (RO, L, S) of the hexagonal and re-entrant lattice without damping ( e = 0 ) and with damping ( e = 500Mpa.sec and e = 800Mpa.sec ). The results in Fig. 2.4 show shifts in the frequency band diagrams due to the presence of damping, these shift being more pronounced as the viscosity coefficient e increases. We also observe that the damping ratio values increase when the viscosity coefficient e increases across the entire Brillouin zone, for all three modes of propagation. These results are in very good agreement with those presented in [48,55] as to the damping ratio and frequency band structure (the same result can be observed for the other lattices). It can be seen in Fig. 2.4 that the Rotary branch is more sensitive to the viscosity coefficient (damping) and presents a higher value of damping, in comparison to the L and S modes (for the 2 lattices). Let recall that the secant slopes of the line connecting the origin point O to the point of interest on the curve gives the phase velocity. Thus, the two lines starting at point O and with slopes C l and C t constitute the best approximation to the dispersion curves in the frequency diagram for n  (without damping), for each lattice [13]. In the case of a damped medium, a decrease of the slope occurs for each lattice (corresponding to a decrease in the frequency), since the phase velocities in a damped medium are lower in comparison to the non-damped medium. We present in Fig. 2.5 a comparison between the dispersion relation and damping ratio for the classical and reentrant hexagon lattices in a damped medium with an extensional viscosity e = 500Mpa.sec. larger than the band gap in the hexagonal lattices and the opposite situation occurs between L and S modes.

The comparison between both lattices highlights that the reentrant lattice is able to dissipate more energy than the hexagonal lattice (Fig. 2.5), since the reentrant lattice shows damping frequencies lower compared to the hexagonal lattice for each mode.

We represent in Fig. 2.6 a comparison of the damping ratio evolutions and damping frequency band structures for the 3 modes of propagations, for the 4 considered lattices (hexagonal, reentrant hexagonal, hexachiral and diamond chiral) in the damped medium for an extensional viscosity e = 500Mpa.sec . The diamond chiral lattice presents the lowest damped frequency in the 3 modes; the hexagonal lattice presents the higher damping ratio for the rotary and longitudinal mode. We also observe that the re-entrant effect for the hexagonal lattice affects the shear mode [52] and this lattice presents a high damping ratio and frequency band structure, close to those of the diamond chiral lattice.

ice. We analyze in Fig. 2.7 the influence of the slenderness ratio  on the damping ratio for the four studied lattices.

Phase velocity in a damped medium

The phase velocity in a damped medium can be expressed as follows:

   C u p k  (2.39)
where k is the wavenumber (the modulus of wave vector), the complex frequency and u is a unit vector indicating the direction of wave vector.

The modulus of the phase velocity is plotted in Fig. It can be observed from Fig2.8 that for waves travelling in a damped medium, shifts in the phase velocity occur for the three modes, these shifts being more pronounced as the wavenumber is increased; they vanish for a specific value of k for which no wave propagation occurs (for the hexagonal lattice, the RO mode vanishes for k=7). These shifts in the phase velocity are more pronounced for the re-entrant and hexagon lattices; this can be explained due to higher values of the damping ratio ζ for the three modes of propagation for these lattices (Fig. 2.5), whereas lower values are obtained for the two other lattices. For an elastic medium (without damping, µe=0), the phase velocity for the 3 modes does not depend on the wavenumber for all studied lattices (as we find in Christoffel equation for the undamped medium); it only depends on the direction of propagation.

Eigenvectors for the classical hexagon lattice

Inserting the expression of the frequency, (Eq.2.37), in the general form of the expression of plane wave (Eq.2.34), we can conclude that the eigenvectors (the wave amplitude) depend on time as follows (the expression of the amplitude with respect to V and  are same as for U)

 

    

; after an elapsed time t =5τ, the lattice has returned to equilibrium (it is clear that the lower the characteristic time, the faster the unit cell reaches equilibrium).

We show in Fig. 2.9 the eigenvector of the hexagonal lattice at initial time and at the subsequent At point O, the unit cell exhibits only a rigid body translation for the 3 branches; at points A and B, the unit cell exhibits both translation and rotation. The first column describes the transverse wave motion; the second gives the second branch of the dispersion curve, representing a longitudinal wave. In the rotational mode, the Eigenwave exhibits a combined transverse and longitudinal motion.

Bloch theorem

The characteristics of wave propagation within a structure are finally evaluated based on the generalized Bloch theorem, in order to investigate the range of validity of the homogenized equations in terms of the frequency range. Bloch theorem states that the amplitude of wave propagation does not depend on the location of unit cell in the structure and accordingly, the characteristics of the periodic structure can be identified by an analysis performed over the reference unit cell only. The wavevector k is defined in a reciprocal lattice (also periodical) in the wavevector space whose basis vectors * j e are given by:

* . i j ij   Ye
with the symbols ij  therein denoting the Kronecker delta function.

According to Bloch theorem, the displacement of a point P in another cell is obtained by translating by a factor n 1 the period in direction 1 Y , and n 2 times in direction 2 Y , thus it can be written in term of the displacement of the reference unit cell as:

1 1 2 2 ( ) ( ) ( ) ,
k n k n q p q r e   (2.41)

Construction of stiffness, mass and viscosity matrices for 2D Bernoulli beams

The continuum counterpart of the kinematic variables interpolating the discrete kinematics of a typical beam element is approximated by

    6 1 , ( ) nn n u x t a x q t        6 1 , ( ) nn n v x t b x q t   
The functions   n ax ,   n bx therein are the shape functions that satisfy the boundary conditions at the end of the beam. We next as a first step evaluate the stiffness and damping matrices based on the computation of the elastic and dissipative energies.

The elastic energy stored in a beam [56,57] is associated to uniform extensional strain and flexion; it can be evaluated separately for each deformation mode. Under extension, the elastic strain energy writes

    6 n 2 2 x n 1 0 0 ax 1 1 EA(ε ) dx EA[ q t ] dx 2 2 x L L e n x x U           (2.42)
The stiffness coefficient (for indices 1 and 4 corresponding to the horizontal displacement components) are identified to

        2 6 n ij n 1 i j 0 0 ax 1 1 K = EA q t EA q q 2 x 2 L L j i n x x da x da x dx dx dx dx                   (2.43)
Under bending, the elastic stored energy in a Euler-Bernoulli beam is the strain energy expressing as (shear deformations are neglected, thus the following kinematic constraint holds,

2 xx 2 v u y θ v = θ, ε = = - = -y x x x x                   2 n 2 zn 1 00 2 6 2 bx 11 EI( ) dx EI q t dx 22 x LL f n xx U             (2.44)
The associated stiffness coefficients (for indices 2, 3, 5 and 6) are

        2 2 2 2 6 j n i ij n 2 2 2 1 ij 0 0 bx b x b x 1 1 K = EI q t EI dx q q 2 x 2 x x L L n x x dx                       (2.45)
in which E is the Young modulus, A the beam cross section area and I the moment of inertia. The dissipated energy due to the damping effects can be written in a form similar to the elastic energy [57,58] and may be evaluated separately for extension and bending effects. Under extension, the dissipated energy writes    

        2 6 j i ij 1 ii 0 0 ax ax 1 1 C= eA eA dx 2 2 x x L L n n n x x a x q t dx x q t q                        (2.47)
Similarly, under bending, we obtain the dissipated energy

    2 2 6 z f 1 0 0 1 1 D= dx e I e t I 2 2 n L L n n x x x q t d x b x t                      (2.48)
This leads to the viscosity coefficients (for indices 2, 3, 5 and 6)

        2 2 2 2 6 j n n i ij 2 2 2 1 ij 0 0 bx b x b x b x 1 1 C= e I e I dx q q 2 x t 2 x x L L n x x dx                         (2.49)
The equations of motion are obtained by applying the Euler-Lagrangian equations without external forces:

, T U D t q q q               (2.50)
This equation of motion can be written for each beam; the equation of motion for the unit cell can then be derived in the following form:

0,    q q q M K C (2.51)
where M, K and C are respectively the assembled global mass, global stiffness and global dissipated matrices, and ,, qqq are the acceleration, velocity and nodal displacement vectors. The resulting dispersion relations can be conveniently represented by letting k vary along the contour of the first Brillouin zone, as described into details in [13]. They are used to quantify the error introduced by the homogenized method for both damped and undamped cases, and for the regular hexagonal and diamond lattice, considered here as reference examples.

Wave propagation in hexagonal and hexachiral lattices

Inserting the generalized equation of plane harmonic waves (2.34) into the equation of motion leads to the following equation , ,

k k k k q e q q e q q e q     
After introducing these relationships into the equation of motion for each of these two lattices, we obtain the final form of the equation of wave propagation within the structure, depending upon the wavenumber and complex frequency.

Hexagonal lattice Hexachiral lattice

1 12 2 ( , , ) 0 q N k k q      1 2 3 1 1 2 4 5 6
( , , ) 0

q q q N k k q q q                    
The eigenvalue problem for N and The main purpose for presenting the generalized Bloch theorem is to delineate the validity of the homogenized method at low frequency for a damped periodic structure. We show in Fig. by applying Bloch's theorem and the homogenization method, at the edge of the first Brillouin zone for both hexagonal and hexachiral lattices for a damped medium with µe = 500 Mpa.sec.

Conclusion

We have computed the band diagram structure and damping ratio for different repetitive lattices, based on the homogenized continuum response of the initially discrete lattice architecture, modeled as Kelvin-Voigt viscoelastic beams. Three of these lattices (reentrant hexagonal, chiral diamond, hexachiral lattice) are auxetic metamaterials since they show negative Poisson's ratio.

The methodology developed in this chapter is general and is applicable to any lattice in the low frequency range: the effective viscoelastic anisotropic continuum behavior is first computed in terms of the homogenized rigidity and viscosity matrices. In a second step, the dynamical equations of motion are written, based on the homogenized properties, and the dispersion relation and damping ratio are obtained by inserting harmonic plane wave expressions into these equations. The comparison of the acoustic properties of the four lattices has been done considering that all lattices share the same beams, with identical microstructural parameters (geometry and mechanical moduli). We have represented the eigenvectors corresponding to the different Eigen modes for the hexagonal lattice, and the application of Bloch theorem in dissipative medium.

We have shown that the auxetic lattice attenuates waves at lower frequencies compared to the classical hexagonal lattice. The diamond chiral lattice shows the best attenuation properties of harmonic waves over the entire Brillouin zone. The hexachiral lattice shows better acoustic properties than the reentrant hexagonal lattice.

The range of validity of the effective continuum obtained by the discrete homogenization has been delimited by comparing the frequency band structure of this continuum with that obtained by Bloch theorem.

The slenderness ratio of the beams is shown to have an important impact on the attenuation properties; therefore, the consideration of Timoshenko thick beams constitutes an important perspective of this work. The versatility of the developed method enables to explore a wide range of 2D and 3D auxetic lattices and more generally metamaterials exhibiting interesting acoustic properties in future contributions.

Appendix A: expression of the stiffness and damping matrices

For homogeneous isotropic viscoelastic Bernoulli-Euler beams obeying a Kelvin-Voigt behavior, the stiffness and damping matrices are given by K=

        εb b εb εb b b β b b b β b b β b 1 2 1 2 1 2 εb εb AE A e N = ε Δ + ε Δ + ε Δ + ε Δ = εN + ε N ll e U U e U U (B.1)           Rb R(b) R(b) R(b) R(b) R(b) R(b) R(b) b εb εb E O E O εb b β b ib 0 1 2 0 0 1 1 3 i εb εb εb O E O E b 2 b ib 0 1 2 0 0 1 1 3 i εb θ 12E I l T= ε Δ + ε Δ - θ + θ + ε θ + θ + 2 l θ 12 eI l + ε Δ + ε Δ - θ + θ + ε θ + θ + = 2 l b b                                         e U U e U U b 2 b 12 εT + ε T (B.2)               Rb R(b) R(b) R(b) R(b) R(b) R(b) R(b) b εb b εb E O E O O(b)ε b β b ib 0 z 1 2 0 0 1 1 2 εb i εb εb εb O E O E b 2 b 0 1 2 0 0 1 1 2 εb i εb θ 6E I E I M = ε Δ + ε Δ + 4θ + βθ + ε 4θ + βθ + β l l θ 6 eI eI + ε Δ + ε Δ + 4θ + βθ + ε 4θ + βθ + β l l b b b b                           e U U e e U U e ib 2 O 3 O 1z 2z = ε M + ε M                       (B.3)               Rb R(b) R(b) R(b) R(b) R(b) R(b) R(b) b εb b εb E O E O E(b)ε b β b ib s s 0 z 1 2 0 0 1 1 2 εb i εb εb εb O E O E b 2 b 0 1 2 0 0 1 1 2 εb εb 6E I E I θ M= ε Δ + ε Δ + βθ + 4θ + ε βθ + 4θ + 4 l l θ 6 eI eI + ε Δ + ε Δ + βθ + 4θ + ε βθ + 4θ + 4 l l b b b b                           e U U e e U U e ib i 2 E 3 E 1z 2z = ε M + ε M                       (B.4)
Since each beam is self-equilibrated, forces and moments can be recovered at the origin node, from the writing the equilibrium equations. Those expressions of the resultant forces and moments are next involved in the homogenization of the initially discrete lattice towards an equivalent continuum.

Appendix C: expression of the first and second order stress vectors

The first and second order stress vectors 1 2 , i i S Sand the first and second order couple stress vectors express as the following sums over all beams of the reference unit cell

            Rb Rb Rb Rb b O E b b 3 b b 1 1 0 0 1 b O E b b 3 b b e 1 e 1 0 0 L E Δ + E Δ - θ + θ + 2 L µ Δ + µ Δ - θ + θ 2 i b Br S                                e U U e e U U e (C.1)         Rb Rb Rb Rb b O E b b 3 b b ib 0 2 2 1 1 i i 2 b O E b b 3 b b ib 0 e 2 e 2 1 1 i θ L E Δ + E Δ - θ + + θ 2 S= θ L +µ Δ +µ Δ - θ + + θ 2 b Br                                         e U U e e U U e (C.2) Rb Rb Rb Rb b 2 b 2 E O E O i 3 3 1 0 0 e 0 0 (L ) (L ) = E (θ -θ ) + µ (θ -θ ) 12 12 b Br          (C.3)     Rb Rb Rb Rb 2 2 b b E O E O 3 ib 3 ib 0 0 2 1 1 e 1 1 i i L L θ θ E θ -θ + +µ θ -θ + 12 1β i b Br                            (C.4)
1-Select any lattice point in the reciprocal lattice as the origin and connect it to neighboring points.

2-Construct the perpendicular bisectors of these lines; the region of intersection is the first Brillouin zone. The Table 2 shows the basis vector of direct and reciprocal diamond chiral lattice and the coordinate of the edge of IBC.

Y = 2 2Li Y = 2 2L j * 1 * 2 1 e = i 2 2L 1 e = j 2 2L   1 1 1 O 0, 0 ; A 0, ; B , 2 2 2             a) b)

Chapter 3: Wave propagation in 3D viscoelastic auxetic and textile materials by homogenized continuum micropolar models

Summary

We extend the work presented in the previous chapter and analyze the impact of wave damping on the dispersion features of 3D viscoelastic periodic structures modeled as Kelvin-Voigt viscoelastic beam-lattices. The band diagram structure and damping ratio of these lattices are computed for different repetitive structures, based on the homogenized micropolar viscoelastic response of the effective medium obtained by the homogenization of the initially discrete lattice architecture. The employed methodology is herewith exemplified for the cases of the 3D hexagon structure, which shows negative Poisson's ration for reentrant configurations, and the 2D plain weave textile structure. The effective viscoelastic anisotropic continuum behavior of both structures is first determined in terms of the homogenized stiffness and viscosity matrices, based on the discrete homogenization method. The dynamical equations of motion are obtained for an equivalent homogenized micropolar continuum evaluated based on the homogenized properties, and the dispersion relations and damping ratio are obtained by inserting an harmonic plane waves Ansatz into these equations. The 3D reentrant hexagon is an auxetic metamaterial showing excellent wave absorption capability, as reflected by the high incidence of partial band gaps. The textile plain weave structure shows a complete band gap for low frequencies. Comparing the frequency band structure of the effective continuum with that obtained by Floquet-Bloch analysis extended to a 3D context shows that the wave propagation characteristics can be directly obtained from the homogenized continuum at low frequencies.

Overview

A large diversity of novel cellular solids and lattice materials have been developed for a wide range of engineering applications, due to the fact that such lattice-like materials show an enhanced static and dynamic response in comparison to their bulk counterpart. They have accordingly raised the interest of many researchers regarding both their static and dynamic behavior. Engineering applications using lattice materials include light-weight structures, vibration control devices or energy absorption systems. The enhancement of acoustic properties is dependent on the bulk material, the relative lattice density and the internal lattice architecture [1]. Since it is usually too expensive from a computational viewpoint to analyze the dynamics of complete repetitive networks due to the huge number of structural elements, a convenient strategy is to represent them as an equivalent homogeneous material, based on a representative unit cell (RUC) of the initial network with effective mechanical properties. Suitable homogenization schemes are then developed to provide the RUC with an effective behavior including information about the microstructure at a continuum level of description.

A well-known drawback of the classical theory of elasticity lies however in its inability to predict in a correct manner the dispersion of waves at very short wavelengths, since it predicts that frequency and wave number are linearly related, so that the phase and group velocities of propagating waves are independent of the frequency; it thus fails to model the impact of microstructural features on the wave dispersion characteristics. This restriction is due especially to the lack of internal length scales which are required when the characteristic wavelength of the exciting incident field is comparable to some internal characteristic intrinsic length of the medium. In such situations, averaging of the individual response of the carriers of the microstructure does not make sense, and the axiom of locality becomes invalid.

Many advanced theories and models have been proposed to study wave propagation problems accounting for non-locality and microstructural effects in materials. One category of approaches is multiscale homogenization techniques [2], which compute asymptotically the solutions of the wave equations involving multiple spatial and temporal scales, to capture the long-term response of the homogenized response. Following another route, the mechanics of generalized continua [3][4][5] accounts for the non-locality of the elastic fields due to microstructural effects in a macroscopic manner, by introducing higher-order strain gradients or additional degrees of freedom, like in micropolar and couple stress theories. Pioneering works in this direction include contributions of Mindlin [6][7][8], Aifantis [9][10][11][12]in connection with the higher-order strain theories, which enhance the classical continuum theories by considering higher order gradients of either the displacement or stain fields, accompanied by internal length scales [13,14]. Different homogenization techniques have been proposed to build continuum descriptions of discrete materials [15,16]; an extended list of references and an overview on gradient elasticity theories applied in dynamics can be found in recent works [17][18][19][20][21][22].

Following another strategy, micropolar theories incorporate additional rotational degrees of freedom, and have been developed by Eringen and co-workers [23,24], a subclass of models being the couple-stress theories [25][26][27][28].No classical theories have proven successful in dynamic problems of wave propagation in beams and half-spaces [29][30][31] by removing singular behaviors (like discontinuities or singularities of the fields) inherent to classical theories and showing their ability to capture the expected size and wave dispersion effects specific to the simultaneous existence of several lengths scales. Analytical studies of the dispersive character of such high order continua have been initially conducted in depth by Mindlin [6], Eringen and Suhubi [23],

and more recently by Engelbrecht et al. [32,33], Papargyri-Beskou et al. [34], Fafalis et al. [35] and Berezovski et al. [36] for linear elastic solids.

In this chapter, we analyze the dynamics of periodic network materials relying on dedicated homogenization techniques developed to substitute to the initial discrete periodic lattice an effective micropolar viscoelastic continuum. The essential objective is to analyze the impact of the viscous dissipation and of the microstructural degrees of freedom of the effective medium on the wave propagation characteristics. We have selected in this chapter two network materials for their expected specific wave propagation characteristics: the 3D hexagonal lattice, having a reentrant configuration associated to an auxetic behavior and showing an increase number of partial band gaps, and a 2D plain weave textile structure for which two modes of wave propagation completely vanish.

The main originality advocated in the present chapter is the derivation of the Bloch theorem in a 3D context for a dissipative viscous periodic lattice; it allows a validation of the dispersion relation obtained for the effective continuum in the low frequency range. We study the influence of the effective retardation time of the homogenized continuum on the dispersion relations and band diagrams, based on the homogenization of the initially discrete repetitive networks. We show a complete attenuation of the first rotary wave for the two studied lattices. The expressions of phase and group velocities in damped medium are obtained.

Construction of 3D viscoelastic micropolar effective continuum

The design of periodic elastic materials such as phononic crystals for wave propagation and acoustic applications has pushed the technological barriers and has opened a new realm of applications in which damping results in the temporal attenuation of the propagating waves. In many situations of interest, at least one of the constituents show a dissipative behavior, for instance the viscoelastic matrix phase in phononic crystal composites.

Studies of the propagation of elastic waves in damped media started more than 40 years ago with an analysis of damping in infinite periodic structures by Mead [37]. From this early period, studies considering partial aspects of the effect of viscous damping on the band structure appeared in the specialized literature; see [38] and references therein, without however providing a detailed analysis of the effects of damping on the dispersion band structure.

The overview of the existing literature on wave propagation in phononic crystals and metamaterials in presence of viscous damping reveals that works in this field started by considering wave propagation with damping in finite structures, see [39][40][41].The static and dynamic properties (frequency band gaps) of anti-tetrachiral cellular solids are analyzed in [42]. The homogenization of periodic hexachiral and tetrachiral honeycombs is performed based on both micropolar homogenization and a second gradient homogenization scheme in [43].

In a second stage, the impact of damping on the band structure of infinite periodic crystals was studied, considering either free wave propagation (for impulsive loading), or time harmonic wave propagation (in case of forced harmonic loading); an overview of the works considering these two situations is presented in the recent contribution of Wang et al. [44]. For free wave propagation, the wavenumber is a real number indicative of propagating models, while the frequency is a complex number, in which the imaginary part measures the time damping.

In [44], the authors analyze the propagation of elastic waves in acoustic metamaterials based on locally resonant viscoelastic phononic crystals, and show the dispersive and dissipative effect of viscosity on the complex band structure and transmission spectra.

Expressions of forces

We consider a homogeneous isotropic viscoelastic 3D beam of Kelvin-Voigt type, slender enough so that the shear deformation can be neglected (Euler-Bernoulli beam). Normal and transverse efforts as well as moments are exerted on the two beam extremities.

The forces at the extremity of the beam can be related to the kinematic variables in the global coordinate system as follows

             p s s F K T U C T U (3.1)
Where  K is the stiffness matrix [45], 

C the viscosity matrix (evaluated in the Appendix A)

and  

T the local to global transformation matrix, consisting of four 3x3 rotation matrices R.

Therefore, based on the framework of viscoelastic Bernoulli beam theory, the efforts and moments can be expressed as The coefficients s E and s G are the tensile and shear modulus of the beam material respectively.

        . .     e u u e u u E s x x E O x E O AE Ae F L L  (3.2)             3 3 12 12 . . . . 2 2                     e u u e e u u e E sz z y y E O z O E x E O z O E EI eI L L F L L      (3.3)             3 3 12 12 . . . . 2 2                     e u u e e u u e s y y E z z E O y O E z E O y O E E I eI L L F L L      (3.4)         ..     ee O s x x O E x O E GJ J M LL      (3.5)         ..     ee E s x x E O x E O GJ J M LL      (3.6)                 2 2 6 6 . . 4 2 . . 4 2         e u u e e u u e s y s y y y O y z E O y O E z E O y O E E I E I eI eI M L L L L       (3.7)                 2 2 6 6 . . 2 4 . . 2 4         e u u e e u u e s y s y y y E y z E O y O E z E O y O E E I E I eI eI M L L L L       (3.8)                 2 2 6 6 . . 4 2 . . 4 2           e u u e e u u e O s z s z z z z y E O z O E y E O z O E E I E I eI eI M L L L L       (3.9)                 2 2 6 6 . . 2 4 . . 2 4           e u u e e u u e E s z s z z z z y E O z O E y E O z O E E I E I eI eI M L L L L       (3.
Parameters A and L are the cross-sectional area and length of the beam respectively; I y , I z are the quadratic moments of the beam and J is the torsional constant.

For simplicity reasons, a circular section of the beam is considered (with constant radius r).

Hence, the cross-sectional area is equal to

Homogenization steps

We will use the same methodology for the discrete asymptotic homogenization method presented in the chapter 2for 3D structures with elastic beams (for details see [45,46]). This method requires the development of all geometrical variables (length, thickness, width) and kinematic variables (displacements, velocity and rotations at the lattices nodes) as Taylor series expansions versus a small parameter labeled ε, defined as the ratio of unit cell size to a macroscopic length characteristic of the entire lattice. These expansions are thereafter inserted into the equilibrium equation of forces and moments, expressed in weak form. After resolution of the unknown displacements and rotations in the localization problem posed over the identified reference unit cell, the stress and couple-stress tensors are constructed versus their conjugated kinematic variables, respectively the strain and curvature tensors, thereby defining the homogenized constitutive law [47][48].

In the local coordinate systems attached to each beam, the asymptotic expansion of forces b

x F  , b y F  , b z F  and moments ( ) ( ) ( ) , , , O b E b E b x x x M M M  ( ) ( ) ( ) ,, O b E b O b y y z M M M    and () Eb z M  are exposed in Appendix B.
Finally the homogenization continuum self-equilibrium can be obtained from the discrete equilibrium equations using the form of Riemann integrals

  3 3 0 lim ( )       i i vZ F v F x dx     i . 0      v S i d    (3.11)   i . 0      μ i d    w (3.12)
with i S the stress vector therein decomposed into

    1 2 1 1 1 2 2 2 sec             S S S i i i b b b ib b b b ib x x y y z z x x y y z z b B b B R R first power of ond power of F F F F F F      
e e e e e e (3.13) and with i μ the couple stress vector decomposed into

        2 1 2 1 1 1 1 2 2 2 2 2 2 1 1 2 2 1 1 2 2                            μ μ μ R i i i E O E O ib y y y z z z bB E E O E O ib x x y y y z z z bB R M M M M M M M M M      
e e e e e (3.14) The Constitutive law: stress and couple stress tensors versus the kinematic variables,

                    i i i i i 1 2 1 2 i i i i i i 2 i i 2 i 1 2 1 2 i i i i 1 1 1 1 1 1 1 1                                                                       R R R R S S S S S R R R R m μ μ μ μ μ e v e v e v e v EE B B BB K K g g g g g g g g          (3.15)
and can therefore be written as

                              e e v v viscous part elastic part E B E B                                  e e v v viscous part elastic part B K B K m  (3.16)
The constitutive law can further be simplified basing on symmetry properties of the studied lattices: it has indeed been shown [49] that for centro-symmetrical lattices the elastic pseudotensor e B   vanish (and similarly its viscous counterpart v B   also vanishes).

Dynamical equilibrium and characteristic equation

In order to obtain the equations of the group and phase velocities as well as the dispersion relation and the damping ratio in a damped micropolar medium, we write the dynamical equilibrium equations which are based on the constitutive equations of the obtained continuum

micropolar medium * * * * * *                                                         
                    (3.17)
Here, , uv and w are the components of the acceleration vector and ,, In order to obtain the displacement formulation of the equations of motion, the compatibility conditions involving the macro strain components are written as follows , , ,

             x x yy zz z z xy x y z yx y x z zy z y x yz y z x xz xx x z y zx z y xy y u u u u u u u u u      (3.18) , , , , , , , , , , , , , ,  
and the curvature as

                , , , , , ,,  
                      x x y y z z xy y yx x yz z zx x xz z x y z x y yzx                 (3.19)
In order to obtain the velocity and the gradient of the linear and angular velocities, one has to derive the expressions (3.18) and (3.19) versus time.

In order to compute the acoustic characteristics (i.e., the group and phase velocities, the dispersion relation and damping ratio), one adopts for the solution of the dynamical equilibrium equation the form of a generalized plane wave.

For a harmonic wave propagating along an axis in an infinite planar micropolar medium, the generalized displacement field with components ( , , , , , )

x y z q U V W    
at a point r is assumed in the following form (see [50]):

  , . ti q Q e    kr (3.20)
In ( 32), ( , , , , , )

x y z Q U V W     is the amplitude of the wave motion,   1 2 3 k ,k ,k  k
the wave vector and a complex frequency (function of the wave vector) that permits wave attenuation in time. The wave vector k is a complex number: its real part represents the attenuation in a

Cartesian coordinate system in 3D, and its imaginary part represents the phase constants.

For a plane wave without attenuation in the Cartesian coordinate system in 3-D, the propagation constants along the x, y and z directions are 11 sin( ) cos( )

k k i    , 22 sin( )sin( ) ki k    and 3 3 cos( ) k k i    .
In the following the units of wavenumber is m -1 .

Substituting the generalized plane wave (Eq. (3.20)) in the dynamical equilibrium equations, delivers the following algebraic equation:

  1 2 3 k , k , k , 0               x y z U V W D    (3.21) A relation between 1 2 3
,, k k k and  is obtained by solving the eigenvalues problem in (3.21) which represents plane waves propagating in dissipative medium at the frequency λ. The roots of the characteristic equation can be expressed in the following form:

2 ( ) ( ). ( ) . ( ) 1     s s ns ns s k k k i k      (3.
      2 2 2 ( ) ( ) , ( ) 1 , ( )       ns s s ds ns s s s ns k real imag k k real          (3.23)
These expressions allow plotting the dispersion curve for the dissipated frequency and the damping ratio versus the wave vector k, as shall be done later.

Each eigenvalue of (Eq.3.21) has a corresponding eigenvector that describes the displacement of the particles caused by each mode of propagation.

Due to the presence of damping (viscosity), the eigenvectors can be expressed versus time as follows

    n s t n U t U e   (3.24)
with n U the n th eigenvector at initial time, for which a pure elastic response is obtained.

We introduce the characteristic time of the lattice, 1 ( ).

( )

s ns kk    ; after an elapsed time 5 t  
, the lattice has returned to its equilibrium state. The lower the characteristic time, the faster the unit cell reaches equilibrium.

Results: dispersion relations and damping ratio evolutions

The 3D hexagonal lattice, parameterized by the angle  (re-entrant case for 0   ), and a plain weave fabric are considered to exemplify the wave propagation phenomena, as illustrated in Fig.

3.2.

They are modeled as periodic network of structural elements modeled as Bernoulli beams.

For the two lattices, the beams with lengths h, L, Lv are given three distinct retardation times,

quantities 1 2 3 ,,    respectively.

Effective properties of the considered lattices

We evaluate in Fig. In the specific case of 3D hexagonal lattices parameterized by the angle , the area of the periodic cell can be expressed as

    cell A 2 cos θ ( sin θ ) L Lv h L   (3.25)
and the inertia parameters are expressed as

                          3 2 3 2 2 ** 3 3 3 2 3 2 2 3 2 2 ** 2 . 33 22 , , 2 cos θ sin θ β . cos θ sin θ 2 2 2 2 cos ( ) sin ( ) . 3 3 3 3 3 , , 2 . cos θ sin θ β . cos θ sin θ z yx h r L Lv h r h L Lv I L Lv h L L Lv h L Lv h Lv r L r L Lv h II L Lv h L L Lv h L                                     (3.26)
in which  is the bulk density of the beam.

For the plain weave fabric, the area of the periodic cell can be expressed as:

  2 2 4 cos θ A Lv  (3.27)
and the inertia parameters are

              22 ** 2 2 22 2 2 3 3 2 3 ** 2 2 3 1 12 2 3 8 , , cos θ 6 cos θ 8 cos θ 4 λ 8 , cos θ        z yx r L Lv r L Lv I Lv r Lv L Lv L LvL Lv II Lv       
(3.28)

Dispersion diagrams

Based on Eq. (3.23), we represent in Fig. 3.4 the dispersion curves and damping ratio for the obtained five types of branches in a general damping situation, for the three considered lattices. that the L branch is more sensitive to the viscosity coefficient and presents a higher value of damping, in comparison to the SH and SV modes (for the hexagonal and re-entrant lattices). This can be explained based on the homogenized viscosity matrix: we find that the first Lamé coefficient is greater than the second one, resulting in a lower damping coefficient for the shear modes. This explains the cut-off branches in the longitudinal modes appearing earlier than the shear modes (corresponding to a damping coefficient equal to unity). Comparing the three examples, we obtain that the cut-off branch for longitudinal and A2 arbitrary waves (named also rotary waves in the literature [51]) occurs earlier for the hexagonal lattice than for the re-entrant lattice in a damped medium ( e  =200); for the A1 modes, the cut-off branch appears much earlier for the textile structure (for a damped medium with a viscosity e  =100).

It is known that [48,52] the plain weave textile shows very good in-plane mechanical properties but it suffers from some weakness in the out-of-plane direction; this behavior can be clearly observed from the disappearance of the longitudinal and vertical shear modes in Fig. 3.4(c) and 4(d). In order to highlight this behavior, the dispersion relation for the textile is pictured in Fig.

3.5 with a direction of propagation θ=π/6, =π/β (in-plane propagation). We shall note that for the three investigated examples, the global assembled mass matrix for the unit cell plays an important role (as important as the stiffness and damping matrices) on the dynamical characteristics. The total mass of the unit cell increases with the beam density, resulting in shifts of the band structure and therefore shifts of the frequency vibrations of the structure (this result is not shown here).

Phase and group velocities

The phase velocity is the speed at which the single plane wave propagates within the medium.

The phase velocity vector, in a direction of wave vector, is defined by:

   C u p k  (3.29)
where k is the wavenumber (the modulus of wave vector) and u is a unit vector indicating the direction of wave vector.

The group velocity defines the direction of energy flow within the structure; it represents the velocity at which the envelope of a wave packet propagates.

The group velocity vector in the k-space can be given by:

1 2 3 ,,              T k g c k k k    (3.30)
The group velocity can be transformed in the Cartesian space by a straightforward manipulation of vectors: In (3.31), E is a matrix containing the components of the direct lattice basis with respect to the Cartesian reference system. The group velocity is therefore the modulus of the group velocity vector, which has its own propagation angle.

The phase and group velocities in the direction of the wave vector can be written in complex form as follows: The dissipative phase and group velocities are given as follows:

c k k c k i c k c k k c k i c k             (3.
2 2 ( ) ( ) () ( ) ( ) 1 ( ) () ( ) ( ) 1 ( )       p p p ds ds ns s g g g ds ds i ns i s i wk c k c k k wk c k c k k   (3.33)
Let recall that the secant slope of the line connecting the origin point to the point of interest on the curve in the dispersion diagram gives the phase velocity, while the tangent to the dispersion curve at any point gives the group velocity [53].

We note that the group velocity is often lower than the phase velocity; however, in reality the group velocity depends on the properties of the damped medium in which the wave propagates.

In the specific case ( ) ( )

p g s s c k c k 
, the media is called non dispersive; a media in which It can also be deduced from Fig 3 .8 that for waves travelling in a damped medium, shifts in the phase velocity occur for the five modes, and these shifts are more pronounced as the wavenumber is increased, and they vanish for a certain value of k.

We observe that, in the case of an elastic medium (case of no damping, µe=0), the phase velocity in the 5 modes does not depend on the wavenumber k (as we find in Christoffel equation for the undamped medium): it only depends on the direction of wavevector. One can also observe from Fig. 3.8, a, c that for a wavenumber k =15, when the cut of branch occurs in some modes, the phase velocity tends to zero, indicating no wave propagation.

For the SH modes, for the three investigated examples, the phase velocity is the same for both the damped and undamped medium (no shift in the phase velocity in the two directions when we have damping medium). This can be explained as follows: the horizontal shear wave depends on the material rigidity and viscosity in the z-direction for the damped medium (they can be decoupled from the other modes of propagation), and the vertical component of the effective viscosity is small enough so that it does no more influence the phase velocity.

We shall note that we obtain the components of the phase velocity vector by the relations: cos( )sin( ), sin( )sin( ), cos( )

p p p p p p x y z C C C C C C        
with p C therein the magnitude of the phase velocity vector (presented in Fig. 3.8). To highlight the influence of the viscosity coefficient on the group velocity, we evaluate in Figs.3.9, 3.10, and 3.11each component of the group velocity vector for the hexagonal lattice. For each mode of propagation, we have a set of group velocities. As we mentioned before, the group velocity depends on the phase velocity for the wave propagation, so that for each mode, we have a set of group's velocities as shown in Figs.3.9, 3.10 and 3.11. We observe that the horizontal shear wave presents the higher values of the group velocity for each component. As for the phase velocity, there are shifts in the group velocity for a damped medium ( e 200  ) and these shifts are more pronounced when the wavenumber increases (the same behavior is obtained for the other lattices).

1 2 3 ,, k k k                in
One can conclude from the previous three figures that for an elastic medium, the group velocity does not depend on the wavenumber, but it only depends on the direction of the wave vector and on the geometrical and mechanical characteristics of the medium. This feature can also be found in the expression of the group velocity presented in Lane [54].

Transmission loss properties of the homogenized medium

The transmission loss properties of the 3D hexagonal lattice and the 2D plain weave fabric are next evaluated. An incident plane propagating in the direction θ=π/6, =π/6with unit amplitude is sent at the boundary for each lattice. The transmission loss is defined by the following formula:

  10log  ti TL u u (3.34)
where i u and t u are the amplitudes of the incident and transmitted waves respectively. This coefficient represents the fraction of incident energy which is transmitted by the structure.

The transmission loss for both periodic lattices without viscosity ( e0  ) and with viscosity ( e 200  ) is represented in Fig. 3.12 for the two rotary modes. It is seen that the transmission loss in the elastic medium increases as the wavenumber is increased; this increase in transmission is mostly due to the increase of the frequency in the modes. For the viscoelastic medium, the transmission loss decreases as expected, since the presence of viscosity implies attenuation of the wave propagation for each mode, and then attenuation in the transmitted waves. Beyond a specific value of the wave number, we observe that the transmission loss is quite small for the first rotary wave (for the 2 lattices), due to the cutoff branch which occurs at this value of the wave number.

We shall note here that the transmission wave is clearly attenuated inside the complete and partial band gaps.

Bloch theorem

Wave propagation in three-dimensional dissipated periodic lattices is studied through the application of the generalized Bloch theorem, which constitutes an original aspect in this contribution. There is to our knowledge no 3D extension of Bloch theorem including viscous dissipation [55].

Any direct lattice structure in a three-dimensional space can be obtained by translating the unit cell along three linear basis vectors i e .The displacement () j qw for a point in the reference unit cell (with position j w ) can be expressed as a planar wave with complex frequency  as follow:

( . . ) ( ) , j t i k w jj q w q e    (3.35)
in which j q is the amplitude of the wave and k the wavevector. The frequency  shall be a complex parameter to take into account the temporal attenuation of the medium due to the presence of damping (viscosity). The displacement of a point P in another cell, obtained by translating by a factor n 1 the period in the direction 1 e , n 2 in the direction 2 e and n 3 in the direction 3 e , can be written in term of the displacement of the reference unit cell according to Bloch theorem as:

1 1 2 2 3 3 ( ) ( ) ( ) , k n k n k n j q p q r e    (3.36)
In (3.36) i k represents the complex component of the wavevector k; its real part represents the attenuation of a wave when traveling from one unit cell to the neighboring one and its imaginary part represents the phase constant. Bloch theorem states that the amplitude of wave propagation does not depend on the unit cell location in the structure and accordingly, the characteristics of the periodic structure can be identified by an analysis performed over the unit cell only. Note that the wavevector k is defined in a reciprocal lattice (also periodical) in the wavevector space whose basis vectors j b are given by:

.

j i ij be  
with the symbols ij  therein denoting the Kronecker delta function.

Evaluation of the mass, stiffness and viscosity matrices

To obtain the eigenvalue problem and then the dispersion relation, we shall combine the dynamical equilibrium equation of the unit cell and the generalized Bloch theorem.

Each beam of the unit cell is considered as a Euler-Bernoulli beam (neglecting the shear deformation) with 6 degrees of freedom at each nodes, that is 3 translations and 3 rotations   , , , , ,

i i i i xi yi zi q u v w     .
The continuum counterpart of these kinematic variables interpolating the discrete kinematics of a typical beam element is approximated by:

    12 1 , . ( ), nn n u x t a x q t        12 1 , . ( ), nn n v x t b x q t        12 1 , . ( ), nn n w x t c x q t    (3.37)       , , n n n
a x b x c x are the shape function.

The kinetic energy per unit thickness of an Euler-Bernoulli beam (without effect of inertia and neglecting the torsion) is given by:

2 2 2 0 0 0 12 12 .. 11 0 1 1 1 2 2 2 1 ( ) , 2 L L L L r s r s r s rs rs u v w T S dx S dx S dx t t t S q q a a b b c c dx                                         (3.38)
The potential energy is given by:

2 2 2 2 2 2 2 0 0 0 2 2 2 2 12 12 2 2 2 2 11 0 1 1 1 , 2 2 2 1 ( ) , 2 L L L z y L s s s r r r r s z y rs u v w U ES dx EI dx EI dx x x x a b c a b c q q Es EI EI dx x x x x x x                                              (3.39)
The dissipated energy is given by:

                                                             (3.40)
where L, E, Iy, Iz, S, e  are respectively the length, Young modulus, second moment of area with respect to y and z axis, the cross sectional area and the extensional viscosity of the Kelvin-Voigt viscoelastic beam.

The equations of motion are obtained by applying the Euler-Lagrangian equations without external forces:

, T U D t q q q               (3.41)
This equation of motion can be written for each beam; then the assembled equation of motion for the unit cell is derived as:

0,    q q q M K C (3.42)
where M, K and C are respectively the assembled global mass, global stiffness and global dissipated matrices and ,, qqq are the acceleration, velocity and nodal displacement vectors.

Analysis of wave motion

After formulating the global equation of motion for the unit cell, we introduce the generalized Bloch theorem to investigate the wave propagation inside the dissipated periodic lattice. The equation of motion is represented as follows, using a condensed notation:

2 ( ) 0, M K C q N q      (3.43)
According to Bloch theorem, the relationships between the nodal displacements write as follows (see Fig. 3.13): q e q q e q q e q q e q Using these relationships, we can write the following transformation:

k k k k q I q I q qq eI T q qq eI q eI q eI                            (3.45)
For any i th node, the nodal displacement vector is given by   , , , , ,

i i i i xi yi zi q u v w     .
Finally, the governing equation of motion takes the form, with matrix N defined in (55):

1 1 1 2 3 2 2 ( , , , ) 0, h q q T NT N k k k q q                (3.46)
where T h is the Hermitian transpose of the transformation T.

The eigenvalue problem for equation (3.46) gives a relation between 1 2 3 ( , , , )

k k k  that represents
a plane wave propagating at the complex frequency . We shall note that we can restrict the values of the wavevector to the edges of the irreducible part of the first Brillouin zone.

We represent in Fig. 3.14 the frequency band structure for the 3D hexagon on the edge of the first Brillouin zone with damping ( e 400  ) and without damping ( e0  ). As previously mentioned, shifts in the frequency occur for a damped medium, which increase with the viscosity coefficient.

Fig. 3.14 Frequency band structure for a damped and undamped 3D hexagonal lattice.

We observe in Fig. 3.14 that the entire lattice presents a complete and large bandgap between the seventh and eighth branches of the dispersion curves. We note a partial band gap situated between the other branches. It is easy to observe from Fig. 3.14 that g c0  when the wave frequency approaches the full band gap (the tangent to the dispersion relation is horizontal), which means that the energy distributed in the system becomes stationary.

The results show the presence of other interesting wave phenomena, such as partial and complete band gaps, and the veering of modes. Veering consists in two branches coming close to each to other without crossing or touching, and it can occur when the coupling coefficient between different modes is very weak.

The main purpose for presenting the generalized Bloch theorem is to delineate the validity of the homogenized method at low frequency and the band gap presented by the lattice. We show in One can deduce from Fig. 3.15 that the homogenized method is in good agreement with the Bloch theorem, since the error is negligible in the low frequency region (below 0.75 rd/s). The same results can also be found in the case of the damped medium which shows a shift in the dispersion curve (as in Fig. 3.14).One concludes that the wave propagation in dissipative periodic lattices can in the range of low frequencies be directly investigated from the dynamics of the effective continuum constructed by homogenization.

Conclusion

In this contribution, we have analyzed the impact of wave damping on the dispersion features of 3D viscoelastic periodic lattices. We analyzed wave propagation in a 3D hexagonal lattice presenting a reentrant configuration leading to an auxetic behavior, and in a textile monolayer structure selected for its expected excellent wave absorption properties. The band diagram structure and damping ratio have been computed for the two repetitive lattices, based on the homogenized micropolar viscoelastic response of the effective medium obtained by a homogenization method dedicated to networks materials. In a first step, the effective viscoelastic anisotropic continuum behavior of the selected lattices has been determined in terms of the homogenized stiffness and viscosity matrices. The dispersion relation and damping ratio are obtained by inserting a harmonic plane waves Ansatz into the formulated dynamical equations of motion.

From the obtained results and discussions, we can draw the following conclusions:

(1) A non-complete band gap (Fig. 3.6) occurs for the hexagonal lattice between the different branches; the excess in partial band gap exhibited by the re-entrant hexagonal lattice shows special acoustics properties for this auxetic metamaterial. The complete band gap occurring in the textile structure for a direction of propagation θ=π/6, =π/6 highlights the possibility to generate new fabric materials with significantly better acoustic properties.

This band gap disappears for in-plane wave propagation due to the vanishing of modes at specific values of the wave length.

(2) The location of the branches in the frequency domain drops by increasing the retardation time of the homogenized medium and a cut-off branch appears at certain values of wavenumber.

(3) The analytical expressions of the group and phase velocities in a damped medium are evaluated. As for the dispersion relation, a shift in the 2 velocities occurs. For an elastic medium, the group velocity does not depend on the wavenumber, but it only depends on the direction of the wave vector and on the characteristic of the medium.

(4) The transmission loss decreases because the presence of viscosity entails an attenuation of wave propagation for each mode, and thus an attenuation of the transmitted waves.

Beyond a specific value of the wavenumber, we observe that the transmission loss becomes very small for the first rotary wave for the two lattices, due to the cut-off branch which occurs at this value of the wavenumber. The transmitted wave is clearly attenuated inside the complete and partial band gaps.

(5) The entire 3D hexagonal lattice presents a complete and large band gap between the seventh and eight branches of the dispersion curves, and a partial band gap between the other branches. These phenomena occur over the whole range of frequencies.

(6) The wave propagation characteristics can be directly obtained from the homogenized continuum at low frequencies, as shown by comparing the frequency band structure of the effective continuum with that obtained by Floquet-Bloch analysis.

APPENDIX A: Determination of the viscosity matrix for a Kelvin-Voigt beam type

The viscosity matrix is given by:  

L C A e A L                                 3 2 3 2 3 2 3 2
                                                                                            
Eq. (A.1)

APPENDIX B: Expressions of forces and moments for a viscoelastic beam

        b 2 b 2 1 2 1 2 2 12
.

.

bb b b b b s x x x b b bb xx AE Ae F l l FF                     e U U e U U Eq. (B.1)           ( ) ( ) ( ) ( ) ( ) ( ) ( ) b 2 b 0 1 2 0 0 1 1 3 2 0 1 2 0 0 1 1 3 12 . . 2 12 . . 2 Rb R b R b R b R b R b R b R b bb b E O E O b ib sz y y z i b b b O E O E b b ib z y z i b EI l F l eI l l                                                                        e U U e e U U e           2 1 2 b b y y F F                              Eq. (B.2)           ( ) ( ) ( ) ( ) ( ) ( ) ( ) b 2 b 0 1 2 0 0 1 1 3 2 0 1 2 0 0 1 1 3 12 . . 2 12 . . 2 Rb R b R b R b R b R b R b R b bb b E O E O sy b ib z z y i b b b O E O E y b b ib z y i b EI l F l eI l l                                                                        e U U e e U U e           2 1 2 b b z z FF                             Eq. (B.3) ( ) ( ) ( ) ( ) () 0 0 1 1 1 1 2 3 1 2 . . R b R b R b R b bb b O E O E O b ib ib s x x x b i b i O O x x GJ J M l l MM                                                                   e e       Eq. (B.4) ( ) ( ) ( ) ( ) () 0 0 1 1 1 1 23 12 . . R b R b R b R b bb b E O E O E b ib ib s x x x b i b i EE xx GJ J M l l MM                                             e e       Eq. (B.5)               ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) b 2 b 0 1 2 0 0 1 1 2 2 1 2 0 0 1 1 2 6 . . 4 2 4 2 2 6 . . 4 2 4 2 2 Rb R b R b R b R b R b R b R b b b b b E O E O s y s y O b ib y z y b i b b b O E O E y y b b z y b b E I E I M l l eI eI l l                                                           e U U e e U U e          0 2 3 1 2 ib i O O y y M M                               Eq. (B.6)               ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) b 2 b 0 1 2 0 0 1 1 2 2 1 2 0 0 1 1 2 6 . . 2 4 2 4 4 6 . . 2 4 2 4 4 Rb R b R b R b R b R b R b R b b b b b E O E O s y s y E b ib y z y b i b b b O E O E y y b b z y b b E I E I M l l eI eI l l                                                           e U U e e U U e          0 2 3 1 2 ib i E E y y M M                               Eq. (B.7)               ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) b 2 b 0 1 2 0 0 1 1 2 2 1 2 0 0 1 1 2 6 . . 4 2 4 2 2 6 . . 4 2 4 2 Rb R b R b R b R b R b R b R b b b b b E O E O O b ib s z s z z y z b i b b b O E O E b b z z y z b b E I E I M l l eI eI l l                                                            e U U e e U U e          2 3 0 1 2 2 ib O O z z i M M                                Eq. (B.8)               ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) b 2 b 0 1 2 0 0 1 1 2 2 1 2 0 0 1 1 2 6 . . 2 4 2 4 4 6 . . 2 4 2 4 4 Rb R b R b R b R b R b R b R b b b b b E O E O E b ib s z s z z y z b i b b b O E O E b b z z y z b b E I E I M l l eI eI l l                                                           e U U e e U U e          0 2 3 1 2 ib i E E z z M M                                Eq. (B.9)

Chapter 4: Analysis of dispersive waves in repetitive lattices based on homogenized second-gradient continuum models

Summary

We analyze the dispersion of elastic waves in periodic beam networks based on second order gradient models obtained by the homogenization of the initially discrete network, relying on the discrete asymptotic method extended up to the second gradient of the displacement. The lattice beams have a viscoelastic behavior described by Kelvin-Voigt model and the homogenized second gradient viscoelasticity model reflects both the initial lattice topology, anisotropy and microstructural features in terms of its geometrical and micromechanical parameters. The continuum models enriched with the higher-order gradients of the displacement and velocity introduce characteristic lengths parameters which account for microstructural effects at the mesoscopic level. A comparative study of the dispersion relations and damping ratio evolutions for the longitudinal and shear waves has been done for four lattices (the chiral diamond lattice, the classical and reentrant lattices, and the pantograph). The developed model allows analyzing both the effects of damping and internal length scale through the second order displacement gradients on the wave propagation characteristics. An important increase of the natural frequency due to second order effects is observed. For the pantograph lattice, the phase velocity for the longitudinal and shear modes is identical and is not influenced by the direction of wave propagation. The obtained results show overall that the pantograph lattice present the best acoustic characteristics.

Overview

New classes of cellular solids and lattice materials have over the last decade found a wide range of engineering applications, such as light-weight structures, vibration control devices, systems for energy absorption, relying on the fact that such lattice-like materials enhance the static and dynamic responses in comparison to their solid counterpart. This improvement of the properties depends on the bulk material of the lattice, its relative density, and the internal geometrical structure [1] Periodic lattices can be considered as prototype models of many systems whose description can be simplified as assemblies of beam elements rigidly connected or joined by hinges. The dynamic behavior of such periodic network has raised the interest of many researchers, especially due to the use of non-destructive techniques for accessing mechanical properties of the investigated material. Due to the prohibitive cost of computing the dynamical response of periodic networks including many elements (thus a huge number of d.o.f.), it proves more economical to represent the network materials on the macroscale as an equivalent homogeneous material obtained from the homogenization over a suitable unit cell consisting of a rigid joint network of beams.

The homogenized moduli contain information about the microstructure, although in an average sense. Homogenized models based on classical Cauchy-type elasticity theory are however not able to provide realistic predictions of many effects arising from small scale, amongst of them wave dispersion. Classical theories based upon the sole first order displacement gradient lack indeed internal length parameters, characteristic of the underlying microstructure. This explains the success of gradient-enriched theories in capturing microstructural effect on the macroscopic behavior of materials, by including high-order gradients associated to internal lengths representative of the microstructure.

Gradient elasticity theories constitute an extension of the classical equations of elasticity by incorporating additional higher-order gradients of the displacements. Many theories based on an enrichment of the classical elasticity framework by higher order gradients have been proposed in the past to overcome deficiencies of classical elasticity [2][3][4] and plasticity [5][6][7][8] theories in describing both static and dynamic phenomena, including size effects, strain and stress fields in the neighborhood of singularities [9][10][11][12][13] and wave dispersion in dynamics [14][15][16][17][18].

Since the early works in the 1960s by [19,22] there has been an abundant literature devoted to the topic, see the review article of [4], which includes an historical perspective.

The initial research in this field was devoted to the elaboration of gradient theories without paying so much attention to materials or structures that could obey such theories and with little consideration of sound physical interpretation of the new phenomena exhibited by these theories.

Much research has been devoted to model static and dynamic problems arising from the various developed gradient theories; all obtained results showed an enhanced stiffness concomitant to the increase of the gradient coefficients, see e.g. [23][24][25][26] and references therein. Note that different homogenization techniques have been proposed to build gradient elasticity continuum models of discrete materials [27][28][29].

Generalized continuum theories have been shown to offer an attractive alternative for capturing dynamic behaviors overlooked by classical elasticity, especially dispersion relations, [30] and references therein. Applications of gradient elasticity in dynamics have fostered extensive research [31][32][33]. In [30] the dynamic behavior of periodic lattice materials is investigated using an equivalent higher-order continuum model obtained by the homogenization of the equations of motion. Considering dynamic aspects and especially wave propagation phenomena, the impact of the microstructure in heterogeneous materials on the dispersive propagation of elastic waves was first recognized and analyzed in seminal paper of [20], who proposed several non-local continuum models to capture the dispersion relation of planar waves.

Gradient elasticity models are useful to predict the wave dispersion characteristics in heterogeneous or discrete systems, [20,17,34,35]. A thorough analysis of the effects brought about by gradient elasticity models can be found e.g. in [13,32], especially the impact of the higher-order inertial terms [4,5]. A detailed comparison between the dispersive characteristics of various simplified models of gradient elasticity can be found in [6].

The descriptions of microstructural effects analyzed in wave propagation are in most of the time phenomenological in nature, following either the line of micromorphic [36] or second gradient continuum models [13] formulated at the macroscale. The phenomenological nature of these models entails that the material parameters involved in the wave propagation formulations need to be calibrated to experiments, and consequently such models are not predictive. The large number of intrinsic parameters to be identified requires a somewhat difficult measurement protocol; this issue is addressed especially in [37][38][39].

In order to circumvent this drawback, multiscale methods have been developed to link the dispersive aspects of wave propagation to the microstructure of the material. Amongst these, homogenization techniques derive the dynamical macroscopic behavior by upscaling the microscopic one, so obtaining a more smeared description of the dispersive wave propagation features, but still reflecting in a predictive manner the impact of the underlying material microstructure. To circumvent the difficulty tied to the large number of intrinsic parameters, a multitude of reduced models with a tractable number of length-scale parameters have been proposed, such as static theories of gradient elasticity involving a single length-scale parameter, used e.g. to analyze stress singularities in the vicinity of the crack tip [2,14,15,40].

One of the biggest challenges of gradient elasticity theories is a physical interpretation of the involved length-scale parameters in terms of given microstructure [4]. An overview of the existing formulations for both static and dynamic problems can be found in [4].

In many studies, gradient elasticity theories have been derived from the continualization of the response of a discrete lattice consisting of discrete masses and springs, see for instance an early work of [41] or the more recent contribution and overview of [4] including many references.

Amongst the employed methods, asymptotic homogenization is quite appealing and has developed importantly in the field of wave propagation over the last decade, [42][43][44]. The long wave propagation in periodic media based on homogenization has been proposed by [45],

extending the approach developed by [46][47][48]. Homogenization in the time domain has also been proposed by several authors, see e.g. [10,49] and references therein, considering the fast and a slow time scales.

Although a wide body of research has been devoted to gradient-enriched theories for both elasticity and phenomena described by internal variables, gradient viscoelasticity theories have deserved much less interest in the literature. One of the very few models of viscoelastic materials with account on length scale effects through strain gradients is by [51], followed by [50], without however a clear relation between the introduced length scale and the material's microstructure in this last work. A viscoelasticity theory with a micro-inertia gradient has been advanced [52] to analyze the effect of both a gradient enhancement and a viscous behavior on wave dispersion in periodic composites.

First strain gradient materials addressed by [20,21] have further interest due to possibility of associating to them higher order inertia terms, such that the kinetic energy depends on the velocity gradient. These effects allow to better describe strain singularities at sharp crack tips and to capture size effects within the dynamic material behavior (dispersion relations) manifested by real materials such as porous materials, polymer foams [4,20,32]. These higher order inertia effects will however not be investigated in the present contribution, which focuses on the impact of second order gradient effects on the dispersion relations obtained for periodic lattices.

The outline of the present chapter is as follows: the homogenized viscoelastic behavior of repetitive planar lattices consisting of viscoelastic Kelvin-Voigt type beams is determined in section 2, based on an equivalence between the writing of the principle of virtual work for the lattice and the posited second-gradient continuum. The constitutive relation for general repetitive lattices exhibiting arbitrary anisotropy is also expanded in matrix format in section 2, based on the introduction of stress and hyperstress vectors reflecting the lattice topology and microstructural parameters. The effective constitutive laws are next introduced into the dynamical planar equilibrium equations (section 3). The dispersion relations and damping ratio evolutions versus the wave number are evaluated for different lattices in section 4, and the phase velocity for both longitudinal and shear waves in section 5. In section 6 we present the effect of the internal length to the dispersion relation. We conclude in section 7 by a summary of the main results and mentioning a few perspectives for future work

Homogenized viscoelastic second gradient behavior of periodic beam lattices

Expressions of forces

We consider a homogeneous isotropic viscoelastic 2D beam with a Kelvin-Voigt type behavior, slender enough so that shear deformation can be neglected. We will describe the deformation of the beam as a function of the end displacements, with these kinematic variables changing in time.

The vector of kinematic degrees of freedom is written as

  1 1 2 2 , , ,  q u v u v
            2 1 2 2 2 3 1 3 3 . . 22 3. b b b b b b b b E e EI M e e I LL                         ee u u u u (4.3)
In these expressions, nodes The kinematic and static variables for any beam in the lattice are represented in the local coordinate system associated to the viscoelastic beam of the so-called simplified Bernoulli model, in the sense that nodal rotations are not considered.

Asymptotic development of the kinematic variables

The beam length L εb can be expressed as the following asymptotic expansion L εb =εL b , and a similar expansion holds for the beam width t εb . The small parameter  is the ratio of the unit cell size to a characteristic size of the entire lattice.

For simplicity reasons, a rectangular section of the beam is considered, with a constant unit thickness e =1. Hence, the cross-sectional area is equal to the width, . The asymptotic expansion of the nodal displacement and velocity, successively quantities  u and

ε d dt   u u
, are written up to the second order versus the small parameter  , in curvilinear coordinates denoted β in the sequel (it is a vector with two components in 2D), as:

        2 0 1 2               u u u u (4.4)         2 0 1 2               u u u u (4.5)
Thus, the displacement difference between the extremity and origin node of each beam is expressed by a Taylor series development versus  , as

          2 0 1 2 1 O O O                   u u u u u (4.6)           2 0 1 2 0 0 1 0 1 2 2 2 22 2 2 2 2 ( ) ( ) () ( ) .. ( ) .. ( ) .. 2 E i ib i ib i ib i ib ii E E E E E b b i i i i i L L L L L L                                                               u u u u u u u u u u u (4.7)
This leads to 22 00 1 11

1 22 2 2 2 2 ( ) ( ) () () 2 i ib E O i ib i E E O E O ii ib b i b L LL UU                                          uu u uu uu uu (4.8)
The relative velocity can be obtained similarly as the relative displacement (or taking directly the time derivative of (4.8)),

0 0 1 1 1 2 2 22 1 2 2 2 2 ( ) ( ) () ( ) 2 i ib E O i ib i E EO ib b EO i i b i L L L U U                                           u u uu u u u uu (4.9)
with i L the periodic length, the index   The normal and transverse efforts as well as the moment exerted on a given beam of Kelvin-Voigt element, in (4.1) to (4.3), can be expressed then under

    2 2 1 2 2 1 2 1 2 . . b b b b b b E b b E E e U U E N U U N N                                e e . . (4.10)
    3 2 2 1 2 2 1 3 2 1 2 . . b b b b b b b E E E b e E T U U U U T T                                 e e . . (4.11)
   

2 3 2 3 1 2 1 2 2 1 3 3 3 2 1 1 . . 2 2 b b b b b b b E b E b b E e M U E U U U M L L M                                   e e . .
(4.12)

Lattice equilibrium equations in virtual power form

The virtual power over the whole lattice express as a sum of the products of internal forces by their conjugated kinematic quantities at each node, as   ()

i i i i bb P T V N U     4.(4.13)
with i U and i V therein the longitudinal and the transverse components of the virtual velocity field, chosen to vanish on the domain edges. The virtual velocity field has the same shape as the displacement field computed in the previous subsection; in the sequel, we select a field of simple form reducing to the first order term in.

The discrete equilibrium of the moments is similarly expressed in virtual power form, such that after development it is self-equilibrated; thus, equilibrium is then automatically satisfied, expressing at the center of each beams as the following equation 117 One can decompose the sum in (4.13) as a double summation on the cells of the lattice and on the nodes of a reference unit cell:

i i i i e c n B c P T V N U P (4.16)
with the set of lattice cells, the set of nodes of the elementary cell, and e P the virtual power of the boundary nodes in an elementary cell, which can be decomposed as

  e O O O O E E E E b P T V N U T V N U      (4.17)
The cell is equilibrated at each of his node, and the efforts are periodical with a period of length equal to the cell width; this means each boundary node n also exists on an opposite edge of the unit cell. One can thus write the following self-equilibrium equations:

O E n n T T    (4.18) O E n n N N    (4.19)
One can thus simplify (4.17) to

      e E E O E E O b P T V V N U U      (4.20)
Let next develop the expressions of the relative longitudinal and transverse velocities  

 E O V V,   EO UU  in a Taylor series expansion:         2 2 2 2 2 ... ( ) ( ) ( ) .. 2 OO O i ib O E O i ib O i ib i i VV V L V V V L V L                              
This entails the relative transverse velocity

2 2 2 2 2 ( ) ( ) () 2 O i ib O E O i ib ii V L V V V L                (4.21)
One obtains in a similar manner the relative longitudinal velocity

2 2 2 2 2 ( ) ( ) () 2 O i ib O E O i ib i i U L U U U L                 (4.22)
Finally, after straightforward computations, expression (4.20) becomes

2 1 1 2 2 2 12 2 3 2 2 2 12 2 2 42 ( ) ( ) ( ) ( ) 2
( ) ( ) 2 OO E i ib E i ib ii i ib O O E E i ib ii e i ib O O E E i ib ii i ib E VU T L N L L V V T T L P L U U N N L L T                                                                                              2 2 2 2 2 2 22 ( ) ( ) 22 b O i ib O E ii V L U N                                                              (4.23)
The continuous formulation of the virtual power is obtained by passing to the limit in the previous discrete sum, assimilating the elementary cell to a small surface element, thus delivering the Riemann integral form of the virtual power when 0   : , denoting therein g the determinant of the Jacobian matrix of the transformation:

119 12 12 [1] [1] det [2] [2] RR g RR              (4.25)
Vector R in (4.25) designates the position vector of any point within the lattice when changing from Cartesian to curvilinear coordinates, expressed under the form

    1 1 1 2 2 2 1 1 1 2 2 2 i,j [1] [1], [2] [2] L Y L Y L Y L Y        R (4.26)   1 1 i,j θ θ cos sin     1 Y   2 2 i,j θ θ cos sin     2 Y
are the periodicity vectors by which one can generate the whole lattice. Thus, the continuum form of the virtual power of internal forces [29] is written as

1 1 2 2 2 1 2 2 2 2 2 0 1 2 2 22 ( ) ( ) ( ) ( ) 2 1 lim ( ) ( ) 2 O O E i ib E i ib i i i ib O O E E i ib i i i ib O O E E i ib i i i E VU T L N L L V V T T L P g L U U N N L L T                                                                                                  2 2 2 2 2 2 2 2 2 ( ) ( ) 22 b ib O i ib O E i i dV V L U N                                                                                               (4.27)
One restricts in the sequel to periodic uniform structures with central symmetry, for which expression (4.27) becomes,

1 1 2 2 2 2 2 2 0 2 2 2 2 2 ( ) ( ) 1 lim ( ) ( ) 2 2 O O E i ib E i ib i i b i ib O i ib O E E i i V U T L N L P dV g L V L U T N                                                                                               (4.28) After developing the quantities, 2 2 2 2 ( ) ( ) ( ) , , 2
O O i ib O i ib i ib i i i V U L V L L                  and 2 2 2 2 () 2 i ib O i LU     
into the Cartesian basis, expression (4.28) then becomes

O O O O E O O O O E E V V V V T L L x y x y U U U U N L L x y x y L P T g                                                                                          2 2 2 2 1
O O O O O O O O O E V V V x y x y V V V L x y x y U U U L x y x y N L                                                                  2 2 2 2 2 2 2 2 2 2 2 2 cos θ sin θ βsinθ cosθ 2 b O O O dV U U U x y x y                                                                                                                                                       (4.29)

Equivalence with a viscoelastic second order grade continuum

The previous formulation of the virtual power of internal forces is next linked with the continuum form of the constitutive equations of a linear viscoelastic second order grade continuum, in order to express the stress and hyperstress tensors together with the continuum equilibrium equations they satisfy.

Second order gradient viscoelastic constitutive law

Let first recall the constitutive law for a homogeneous anisotropic viscoelastic second order grade continuum, written in index format:

                    e e v v e e v v

viscous part elastic part viscous part elastic part

A B A B S B D B D                                             
with , , , , , therein are respectively the first and second order elasticity and viscosity coefficients, the coupling moduli, which all depend on the specific considered lattices.

ij
From a general viewpoint, the density of the virtual power of internal forces writes [53]:

. .

: i P                   D S D    (4.30)
involving the second order Cauchy stress , the third order hyperstress tensor S with the

symmetry ijk ikj SS 
, and the virtual rate of deformation . D , with the gradient operator considered as the vector

i i x     e .
Considering virtual velocity fields that vanish on the domain boundary  and in the absence of body forces and couples, we obtain the following equilibrium equation of the second order grade continuum [53] in both tensor and index format:

  2 0 0, 1..2 kij kp p i j symmetrical S k x x x                          S  (4.31)
This equation can further be developed as the system of two differential equations

2 2 2 2 k1 k 2 k11 k12 k 21 k 22 1 2 1 1 1 2 2 1 2 2 S S S S 0, k 1, 2 x x x x x x x x x x                               (4.32)
We next determine the expression of the stress and hyperstress tensors based on the virtual power of internal forces and highlighting the analogy with the expression of the virtual power obtained for the initially discrete lattice (4.29):

    . i P dV         SD  (4.33)
After expansion, previous expression can be rewritten under the form 2 .

.

i q pq q p q P dV x x x                               DD FH (4.34) With the pair of indices           , 1,1 , 2,2 , 1,2 pq 
here and in the sequel, and introducing the following force and hyper force vectors therein: (4.38) Introducing therein the following vectors

q q ij i j q iq i e e
1 1 2 2 1 1 2 2 , , and
b b b b E E E E E E E E T T N N      T e T e N e N e
The components of these vectors include both elastic and viscous contributions, as can be seen from (4.10). One can reformulate expression (4.29) under the form

2 0 .. lim q pq b q p q P dV x x x                          F H uu (4.39) such that     11 1 1 1 1 2 2 2 θ θ 1, 2 , cos + cos EE q L L g           F TN (4.40)   1 1 2 1 1 1 2 2 2 θ θ sin + sin E E L L g          F TN (4.41)
The scalar g is the determinant of the Jacobian matrix. There not four independent combinations, but only three, since one sums on the combination of the crossed derivatives (1, 2) and (2, 1), , see [55,56] for more details.

thus 2 2 2 2 2 2 2 2 22 11 1 1 1 2 θ θ cos cos + 22 EE L L g            TN H (4.42) 2 2 2 2 2 2 2 2 22 1 1 1 2 2 2 θ θ sin sin + 2 2 E E L L g             TN H (4.43)   θ θ θ cos sin + sin cos E E L L g          TN H (4.
The viscous characteristic lengths are elaborated in a similar manner, using the dissipated energy instead of the elastic energy (we shall not repeat the previous writing), of the form:

ij ij ijk ijk 1 W W , = σ +S 2 vv d d ij ijk
, with , ijk jk i

Dynamical equilibrium and characteristic equation

For any homogenized 2D viscoelastic lattice, the equations of motion for a second gradient medium, write in components form as the two following differential equations along the x and y directions of a Cartesian coordinates system,

                             (4.45) 2 2 2 2 21 22 211 212 221 222 1 2 1 1 1 2 2 2 * 1 2 S S S S x x x x x x v x x x x                              (4.46)
Here, u and v are the horizontal and vertical components of the acceleration vector. The effective density therein is given in general by

* 1 cell M A  
, with M 1 the mass of the set of lattice beams, A cell being the area of the periodic cell; these quantities are next computed specifically for all four lattices analyzed in this work.

In order to obtain the displacement formulation of the equations of motion, the compatibility equations involving the macro strain components are written as follows The gradient of deformation is given by its 6 independent components as

      u u u u u u      
In order to obtain the velocity, and the gradient of velocity, one has to derive the expression there above of the displacement and displacement gradient versus time. The deformation, the gradient of deformation, the velocity gradient and the second order velocity gradient are used in the constitutive equation. The wave vector k is a complex number: its real part represents the attenuation in the x-y plane, and its imaginary part is the phase constants. The eigenvalue problem for Eq. (4.49) yields a characteristic equation developed as:

4 3 2 0 a b c d          (4.50)
The roots of Eq. (4.50) may be expressed in the following form:

2 ( ) ( ). ( ) . ( ) 1 s s ns ns s k k k i k          (4.51)
in which s represents the branch type, namely l standing for the longitudinal waves and t for the shear waves. Two pairs of complex conjugates solutions are obtained, corresponding respectively to longitudinal and shear waves.

    2 2 2 () ( ) ( ) , ( ) 1 ,
s ns s s ds ns s s ns real k real imag k k                126
Relying on these expressions, we can plot the dispersion curve for the dissipated frequency and the damping ratio versus the wave vector k.

4.4Application to specific lattices

In this study, four lattices are selected for their expected interesting acoustic properties:the hexagonal lattice, the re-entrant lattice (in fact the hexagonal with a negative angle), the diamond chiral lattice and the pantograph, as pictured in Fig. 4.2. We first compute the homogenized mechanical response of these lattices, as a basis to determine their wave propagation behavior.

The geometrical and mechanical parameters of the four unit cells are given in Table 4.1.

     1400MPa, 0.3, s E    Pantograph 50mm, 1mm, 90 l t      1400MPa, 0.3, s E    3 1000kg/m   Diamond chiral 50 mm, 1mm l t   1400MPa, 0.3, s E    3 1000kg/m   Hexagonal 50 mm, 1mm l t   1400MPa, 0.3, s E    3 1000kg/m  
The effective area and density of these lattices is evaluated in closed form.

For the hexagonal lattice parameterized by the angle θ, the area is given by  

Stress and hyperstress tensors for the considered lattices

The expression of the stress and hyperstress tensors for the investigated lattices made of Bernoulli viscoelastic beam of Kelvin-Voigt type is evaluated based on the general methodology exposed in section 2.

Hexagonal lattice

One obtains the Cauchy stress expressing as -For the elastic part:

  3 1 2
U V U V a K a K a K a K x y y x U V U V KK a K a K a K a K y x x y                                   
-For the viscous part:

128   3 1 2 1 1 2 1 3 3 6 v v v v v v v v l l f f v vv lf v v v v v v v v f f l l a K a K a K a K x y y x KK a K a K a K a K y x x U V U V U V U V y                              . . . . . . . . Where           1 2 3
v v l f l f l f l f l f a K K a K K a K K a K K a K K           and   3 3 v v v l f a K K  
The hyperstress tensor is obtained in its different components as follows -For the elastic part:

 

S K a K a K a S K a K a K a S K a K a K a L S K a K a K a KK S K a K a K a S K a K a             2 2 2 2 2 2 2 3 2 2 2 0 0 9 ee l U x V x U yx V yx U Ka y V y                                         
-For the viscous part:

  KKand taking the velocity instead of the displacement. We accordingly do not write explicitly the viscous stress and hyperstress tensors (which can be obtained based on this substitution).

v v v v v v v f f f v v v v v v v l l l v v v v v v v l l l v v v v v v v vv f f f lf v v v v v v v f f f v v v v v ll xxx yxx xyx yyx xyy yyy S K a K a K a S K a K a K a S K a K a K a L S K a K a K a KK S K a K a K a S K a K a             2 2 2 2 2 2 2 3 2 2 2 0 0 9 vv l U x V x U yx V yx U Ka y V y       

Diamond chiral lattice

The Cauchy stress is expressed, for the elastic part, as

          
l f l l f l U V U V K K K K K K KK y y x x U V U V K K K K K K K K K K x x y y                                         
The hyperstress tensor writes:

          
U x V K K K S x S K K K U S K K L S K K K K K K K S S K K K              
c UV UV c K c K K c K xy yx UV c UV KK c K c K c K K yx xy                           
The hyperstress tensor writes:

  3 1 3 2 1 2 1 2 1 2 1 3 1 3 1 3
K c K c K c S K c K c K c S K c K c K c S L K c K c K c S KK K c K c K c S S K c K c K c                         2 2 2 2 2 2 2 2 2 2 U x V x U yx V yx U y V y                                         where     1 2
; 3 The constitutive equations are then used in the dynamical equilibrium equation of the effective second order grade medium to obtain the dispersion, damping ratio and phase velocity responses versus the wave vector, for all studied lattices, following the methodology described previously in section 3.

U x V x U yx V yx U y V y S K K S K K S L S S K K S K K              

Evaluation of the internal extensional and shear lengths

The calculated elastic and viscous internal lengths are evaluated for the 4 studied lattices based on the general expressions given in previous subsections (Table 4.2). (3 ) 3 4

( 3 )

t L t L t   2 2 2 2 (3 ) 4 ( 3 ) t L t L t   2 L e xy l 2 L 43 L 1 4 L 2 L v xx l 2 L 2 2 2 2 (3 ) 3 4 ( 3 ) t L t L t   2 2 2 2 (3 ) 4 ( 3 ) t L t L t   2 L v xy l 2 L 43 L 1 4 L 2 L
As shown in Table 4.2, the elastic and viscous internal lengths for the 4 lattices do not depend on the beam mechanical properties ( , e E  ), but only on the geometry through parameters (L, t), respectively the beam length and width.

Dispersion relations and damping ratio evolutions

We firstly evaluate the dispersion relation and damping ratio evolutions of the diamond chiral and re-entrant lattices, in order to show the influence of the viscosity coefficient on both the frequency band structure and damping ratio. In the following e  is given in MPa.sec. entrant lattice in the shear mode, as found in [60], while the re-entrant and the hexagonal lattices present the higher damping ratio of all lattices in the 2 modes.

In Fig. 4.7, it can be observed that the damping frequency for the shear mode of all 4 lattices are close to each other, whereas for the longitudinal modes and at low wave number, the diamond chiral and pantograph lattices present the lowest damped frequency. By increasing the wave number, the damping frequency for the hexagonal and re-entrant lattices tends to zero faster than for the diamond chiral and pantograph lattices, due to their higher damping ratio, as shown in In Fig. 4.11, we compare the dispersion relations for the Cauchy and second gradient media, for the diamond chiral and re-entrant lattices, for the 2 modes. We observe that there is an important shift in the natural frequency between Cauchy effective medium and a second order medium; the natural frequency is increased by second order effects. This is due to the contribution of the second order effective modulus of the entire lattice, given from the hyperstress S, which globally increases the effective elastic properties.

Phase velocity

The phase velocity in a direction can be written in complex form as follows: We observe (Fig. 4.15) that the internal length in the shear case do not affect the natural frequency for the diamond chiral lattice and generate a small shift in the 2 modes. For the reentrant lattice, and for the longitudinal mode, it's found that as we increase the internal shear length an important increase of the frequencies occurs. For the shear mode of the same lattice, the 3 branches coincide.

( ) ( ). ( ) . ( ) 1

s s ns ns s c k k c k i c k       (4.

Conclusion

This chapter provides an analysis of the dispersion of elastic waves in periodic beam networks based on second order gradient models obtained by the homogenization of the initially discrete network obtained by the discrete asymptotic method extended up to the second gradient. The We have obtained an important increase of the natural frequency due to second order effects, which overall increases the effective elastic and viscous properties. For the pantograph lattice, it has been obtained that the phase velocity for the longitudinal and shear modes is identical and is not influenced by the direction of wave propagation. Since the phase velocity has the lowest value for the pantograph, this lattice has the best absorption properties. The presented results entail that the pantograph lattice present overall the best acoustic characteristics amongst the studied lattices.

Future work shall analyze the wave propagation in presence of nonlinear deformation of the lattice due to geometrical nonlinearities arising due to the small bending stiffness of the structural beam elements. The developed homogenization technique shall further lead to continuum models enriched with higher-order inertia terms.

Chapter 5: Wave propagation in random fibrous networks based on generalized continuum mechanics

Summary

This work is done in the framework of collaboration with Kamel Berkache, a PhD student from USTHB University in Alger (Algeria).

The dynamic analysis of random fibrous networks is a novel topic, which has many applications since these networks are very often subjected to dynamical loading conditions such as mechanical vibrations. In order to bypass the complexity of performing dynamic computations at the microscopic scale of random fibrous network, we develop and identify couple stress and gradient models as effective continua at the mesoscopic level of windows of analysis, in order to analyze the size effects of such networks on their dynamic properties. The static mechanical properties which are at the basis of the dynamical analysis are computed thanks to FE simulations performed over windows of analysis subjected to mixed boundary conditions allowing to capture the classical and non-classical effective moduli. The acoustic properties are captured by the dispersion diagrams and plots of the phase velocities; we analyze the influence on the dynamic properties of three main quantities of interest, namely the fiber bending length, the size of the window of analysis, and the fibers density. The impact of these parameters is successively assessed for the couple stress and strain gradient substitution continua. A comparison of the acoustic properties of the two effective media is provided as a summary of the present work. The couple stress medium is essentially non dispersive waves in longitudinal mode and dispersive waves in shear mode, whereas the strain gradient medium behaves in a dispersive manner for both modes. Both media show an increase of the frequency with the fiber bending length. Small variations of frequency and partial band gaps occur with fibers density in the affine regime, these effects becoming important in the non-affine regime. The effect of window size on the dispersion diagram and partial band gaps is weak in both affine and non-affine regimes. In affine regime, the influence of the second gradient disappear and no significant effect in the phase velocity. Significant shifts in the partial band gaps occur for the strain gradient effective medium in comparison to the couple stress continuum.

Overview

Nature exhibits a large diversity of materials presenting a stochastic fibrous microstructure, building highly complex and multi-functional parts. Protein for instance, is frequently found in nature in a fibrous form; the most abundant fibrous protein in mammals is collagen, which constitutes the major part of tendons and ligaments, and most of the organic matrix in bone and dentin. It confers mechanical stability, strength and toughness of these tissues [1].

The structural characteristics of random fibrous networks depend upon the properties of the fibers, thus their modeling is necessary in order to understand the mechanism of deformation and failure on system sub-scale, due to the difficulty of measuring the in situ deformation mechanisms of the fibrous microstructure.

The tradeoff between local axial stretching and bending deformations of the fibers has an important impact on the overall mechanical response; especially, the response of the network to imposed deformations is likely to be non-affine, the degree of non-affinity being controlled by the bending length, a scalar quantity which quantifies the relative importance of the bending to the stretching stiffness, as pointed out in [2][3][4] , who evidenced that the network shifts from the nonaffinely deforming structure to an affinely deforming one by increasing the bending length. [5][6] further concluded that the degree of heterogeneity decreases by increasing the network density; one of the most important result as to scale effects is that the heterogeneity leads to a strong dependency of the apparent moduli on the size of the probed network domain; this has been modeled in [3] by evaluating the correlation functions of the tensile modulus versus the window size.

As described in [7], models in the literature developed to simulate the mechanical behavior of fibrous networks fall into two main categories, namely phenomenological models and micromechanical models. Micromechanical models overcome the shortcomings of the phenomenological models, which very often are not able to capture the relation of the fiber properties to the model parameters. The primary focus of micromechanically based constitutive models of non-woven fibrous networks is the elastic behavior. Cox in [8] was one the first to propose a model for the elastic modulus of paper based on the mechanics of the fiber network, with all fibers extending from one end of the mat to the other, and assuming stretching of the fibers as the dominant deformation mechanism. However, since fibers have a relatively low stiffness and are randomly oriented, bending is an important feature, particularly in the absence of a supporting medium [9][10].

As mentioned in very recent contributions [11][12], the dynamics of periodic fiber networks has recently raised a lot of interest, whereas wave propagation in random fibrous networks has not deserved yet the attention of researchers. Although fibers networks are often subjected to dynamical loadings like vibrations, most of the works have indeed been devoted to the analysis of their static behavior. The contribution of [11] seems to be a pioneering work on the wave propagation in random fibrous networks in the small strain range. The authors show via frequency-domain FE computations performed at the microscopic scale of the entire network that the response is non dispersive at long wavelengths, while it becomes dispersive at intermediate and short wavelengths. In this last situation, the Bloch modes are highly non-affine (the degree of non-affinity is increasing when decreasing the wavelength), most of the deformation localizing within the longest fibers of the network.

In this work, we will use generalized continuum theories at an intermediate mesoscopic level in order to address the issue of size effects related to the random fibrous microstructure; those theories have been extensively used to explain size effects for a wide class of materials, but not for random fibrous networks to our knowledge. This constitutes the main originality advocated in the present contribution. Real materials such as biological membranes and tissues often exhibit a number of important length scales, which must be included in any realistic model.

Many advanced theories and models have been proposed to study wave propagation problems accounting for non-locality and microstructural effects in materials. One category of approaches is multiscale homogenization techniques [13], which aims at computing asymptotically the solutions of the wave equations involving multiple spatial and temporal scales, and to capture the long-term response of the homogenized response. The mechanics of generalized continua accounts for the non-locality of the elastic fields due to microstructural effects in a macroscopic manner, by introducing higher-order strain gradients or additional degrees of freedom. An overview of the literature on the advanced theories and models proposed to study wave propagation problems accounting for non-locality and microstructural effects in materials can be found in [12], in which the authors analyze the dynamic properties of periodic textile structures.

5.2Identification of 2D continuum equivalent moduli based on couple stress medium

5.2.1Generation of different RVE's made of random fibrous networks

In this work, we will consider random structures to be a special class of stochastic fibrous networks and classify a random process as one where the events are independent of each other and equally likely, according to the criteria identified in [21].

The generated random network consists of finite fibers distributed with a random orientation, each passing through a point distributed according to a Poisson point process in a plane. Many fibrous networks are very close to being two-dimensional because the dimensions perpendicular to their plane are very small relative to the in-plane dimensions. We consider here systems of two-dimensional networks in which the fibers are of uniform length 0 L and are deposited on squared regions of dimensions L ; a typical network generated under the above mentioned conditions is shown in Fig. 1. Without loss of generality, the thickness in z-direction is set to unity. linear size L occupying a much bigger network. In the sequel, the wording network density shall refer to quantity D.

The cross-links are introduced at all points where fibers intersect, at which the coordination number is 4 z  ; they are here modeled as ``welded'', so that the fibers are loaded both axially and in bending. The angle between intersecting fibers is preserved and leads to a transfer of the bending moments between fibers. Numerical computations with ``Welded'' joints type are more stable than those with ``Pin'' joints type, as discussed in [6]; nevertheless, simulations show that there is no influence of the type of interactions (welded or free) on the effective classical and nonclassical moduli.

Computation of the effective properties of the couple stress substitution continuum

Generalized continuum theories have been extensively used to account for size effects. There are two classes of generalized continuum theories: the strain gradient theory from the class of highergrade continua and the micropolar theory from the class of higher-order continua.

In the former, the gradient of strain is proposed as an additional deformation measure, while in the latter the microrotation gradient is the source of an extra internal energy.

The development of the non-linear theory of elasticity traces back to the seminal work of the Cosserat brothers that did however not get the attention it deserved for a long time. At the beginning of the 1960s, several authors interested in [22] theories, and a special case of the Cosserat continuum theory has been investigated by [17], in which the rotation of the rigid Cosserat triad is defined in terms of the displacement gradients, deserving the name couple stress theory. In the micropolar theory, the deformation is described by the displacement vector u and an independent rotation vector , whereas in the couple stress theory, the rotation vector  is not independent from the displacement vector, since it is given as the antis-symmetric part of the displacement gradient, leading in the present 2D context to the expression of the microrotation around the z axis as

                              (5.2)
As a result of the kinematic coupling (5.1) in the couple-stress theory, the strain tensor ij  is symmetrical with components defined as 1 2

xy yx v u x y               .
Ignoring body forces and body moments, the dynamical equilibrium in translation and rotation, writes as the set of three equations 0, 0, 0

xy yx yy yz xx xz xy yx m m x y x y x y                          (5.3)
The balance equation of internal bending momentum (5.3) implies the equality of both shear stress components xy yx   .

Thus, the constitutive equation can be expressed in the following uncoupled form (for a centrally symmetric unit cell structure) as

AA AA A m D m D                                                      (5.4)
where in the coefficients ij A are the classical Cauchy moduli, while coefficients ij D are the micropolar moduli relating the two independent non nil couple stress components to the corresponding curvatures. The effective Young's moduli can be expressed versus coefficients ij A as: , .

x y

A A E A E A A A    
Furthermore, the effective Poisson ratios are computed as: The main purpose of this section is then to determine the effective constitutive constants of the couple-stress continuum from the response of random fibrous networks within windows of analysis of different sizes (in the sequel, we conveniently use the short cut WOA for window of analysis). We design different boundary conditions for the determination of the independent components of the constitutive (rigidity) constants over a domain  with boundary  . In each case, we force the WOA to bear a set of specific deformation, as described in [18], and compute numerically the total elastic strain energy WOA U stored in the WOA under the corresponding boundary conditions. The numerical procedure used here is similar to the one used in [18] but here restricted to a 2D situation: the total strain energy stored in the WOA is equated to the energy of an equivalent homogeneous couple-stress continuum, thus it holds the identity

2 WOA couple stress ij ijkl kl ij ijkl kl V U U A D             (5.5)
with V the volume of the WOA. The strain energy stored in the effective homogeneous couplestress continuum can be obtained by the prescribed strain/stress fields.

We define the fiber bending length as the ratio between the axial stiffness to the bending stiffness,

parameter b f f l E I E A 
. The evolution of the classical and couple stress moduli versus the fiber bending length b l is plotted in logarithmic axes in Fig. 5.2 and Fig. 5.3, for a constant density of fibers. Low values of b l enhance local rotations of the fibers, which do not follow the imposed deformation over the boundary of the window of analysis, thus the network responds essentially in a non-affine manner; opposite to this, high values of b l lead to a rather affine response, while Fig. 4 shows that the network is nearly isotropic, since the previously defined anisotropic measure  does not exceed 12% for the classical modulus and the deviation from isotropy is 5% for the couple stress moduli. We can thus conclude that the generated fibrous networks are nearly isotropic (up to statistical fluctuations) for all large enough WOA's.

The dynamical response of such fibrous networks can be evaluated either based on computations done at the microscopic level [11], or at the scale of the homogenized continuum. The first approach however involves huge computations due to the large number of fibers (a few thousands), so we shall instead follow the second route in this contribution, and aim at capturing the influence of the microstructure (which controls the dynamic response of the system) by a suitable enhancement of Cauchy elasticity, considering successively higher order and higher grade effective continua. A justification for the need of such an enhancement of Cauchy continuum will be given later on.

Wave propagation analysis based on couple stress theory

In order to set the stage, the dynamical equilibrium equations for the identified effective couple stress medium corresponding to the two independent degrees of freedom , uv in the present planar situation write successively as the two independent equations (using the balance equation of momentum to express the following derivative of the shear stress components

m m u m m v               (5.6)
The effective density therein is given in general by * Inserting the constitutive law, and adopting a plane harmonic wave Ansatz for the solution of the dynamical equilibrium equation, the wave motion equations are further written as:

    1 1 1 1 2 2 2 2 4 * 2 2 2 A k + A k D k k D k -ρ ω k k βA +A D k D k 11 33
U V                              (5.7)
Nontrivial solution of the last equation exists provided the determinant   The phase and group velocities can be expressed as follows 12 ,,

pg cc k k k            (5.8)
In the case of wave propagation in the longitudinal direction (x direction in the Cartesian basis), the dispersive longitudinal and transverse modes have frequencies expressing versus the wave number and effective moduli as

    2
2A +A k k k 4 A k 2 A k 4A k 4 A k k 2A +A k k k 4 A k 2 A k 4A k 4 A k k l t D D D D A A D D D D A A                    
(5.9)

Expressions (5.9) show that the two propagation modes are controlled by the components of the rigidity matrix in traction along x direction, the shear components along the y direction and the couple stress components along the x direction. Nonlinear relations are obtained between frequency and wavenumber, and the frequency reaches a plateau for large values of the wavenumber.

For in-plane wave propagation   12 0, 0 kk  , the presence of the shear strain components in the two tensile stress components , xx yy  introduces a coupling between the shear and tensile effective moduli, which in turn impact wave propagation for the longitudinal and shear modes.

In the sequel, we shall analyze the influence of internal length b l , density D and window size L on wave propagation within the network. We introduce the following non-dimensional parameters: The frequency band structure can be obtained based on the previous methodology, relying on the computed effective moduli of the homogenized couple stress continuum. In all subsequent plots, the continuous line corresponds to wave propagation in the longitudinal mode, while the dashed line corresponds to the shear mode. All frequency band structure is presented for direction of propagation 4  .

Influence of internal length on the dispersion relation and on phase and group velocities

We shall first investigate the effect of the internal bending length b l on the dispersion relation, and evaluate the phase velocity; recall that b l is defined as the ratio of the bending stiffness to the axial fiber stiffness. values of k, based on previous expressions (5.9). The same evolutions of the frequency are obtained in [11], in which the authors perform finite element analysis of the dispersion relations at the fiber level. Works from the literature show that anisotropy of wave propagation changes with frequency, and is not in conflict with the static isotropy of the network. Furthermore, Fig. 5.6 shows that the state of anisotropy is not affected by the internal bending length and wavenumber.

In the longitudinal mode, the medium is non-dispersive as shown in Fig. 5.6 (there is no change in the phase velocity when wavenumber increases), while the shear mode is dispersive (the phase velocity increases with wavenumber). This behavior can be explained by the additional rotational degree of freedom  of the couple stress theory which leads to a second gradient terms in flexion and shear and thus only affects the shear modes, but not the longitudinal mode (which is independent from  ). The dispersive behavior of the shear modes in the couple stress theory is the basic difference in comparison with Cauchy medium which is non-dispersive for both longitudinal and shear modes.

Influence of the fiber density on the dispersion relation and on the phase and group velocities

It is well recognized that the dispersion relation and phase velocity depend on the static properties of the fibers. In order to assess the effect of the density of fibers on the dispersion diagram, we represent in Fig. 5.7 the frequency band structure versus wavenumber for different fiber density values in both affine and non-affine regimes. In the non-affine regime, an increase in frequency and partial band gap width occurs; in contrast to this, for affine deformations, the frequency and the width of the partial band gap remain constant, and all data converge to a horizontal asymptote (stabilization effect), so that the long wavelength speeds of the networks become independent of fiber density. The anisotropic behavior for the longitudinal and shear modes is highlighted by the polar plot of the modulus of the phase velocity for different values of the network density (Fig. 5.8, 5.9), in the affine and non-affine regimes, for two values of the wavenumber. We conclude from these plots that the anisotropic behavior is the same in both regimes, and it is furthermore independent from density. The degree of anisotropy becomes higher for the longitudinal and shear modes when moving from the affine to the non-affine regime, whatever the value of density. In the non-affine regime, the medium remains non-dispersive for the longitudinal modes and dispersive for the shear mode when density increases. In the affine regime, the medium moves from a dispersive to a non-dispersive response as the density increases: this can be explained by the fact that in the affine regime and for increasing density, the macroscopic deformation is close to the microscopic deformation, which entails that the medium behaves as an effective Cauchy medium showing no influence of the microstructure (no dispersion).

Influence of window size on acoustic properties

The evolution of the frequency predicted by the couple stress theory versus window size is pictured in Fig. 5.10 in the affine regime (for large b l ) and non-affine regimes (for small b l ), for the two modes of propagations and for two values of the wave number. The evolution of the network anisotropy and the dispersive character of waves are analyzed from the polar plot of the modulus of the phase velocity in Fig. 5.11, 5.12, in both affine and nonaffine regimes. It shall be pointed out that there is in fact no unique choice for the effective substitution continuum of the initial random fibrous network, since computations performed at the microscopic level of the lattice reveal both important local rotations of the individual fibers, as well as internal strain gradients. Models based on couple stress type theory do however not give realistic predictions of the effective medium properties, such as the dispersion relation (since the longitudinal mode is non-dispersive). It has however been proven by experiments that most waves are dispersive. Accordingly, we explore in the sequel an alternative modeling strategy, considering an enhancement of the Cauchy continuum by higher strain gradients.

Identification of 2D continuum equivalent based on second gradient medium

In classical continuum mechanics, only the first displacement gradient is involved and all the higher order displacement gradients are neglected in measuring the deformations of a body. In this case, the stress tensor at a material point is linked to the strain tensor through the classical elasticity tensor. The second gradient elasticity is a kinematic enhancement of classical elasticity taking into account the second gradient of the displacement field.

The strain gradient theory is next exploited in terms of the strain (first displacement gradient) and second gradient of the displacement field which are, respectively, the second and third order tensors

        ,, , , , 1 1
: 2 2 1 : 2 ij i j j i ijk ij k i jk j ik uu K u u                     ε u u

Ku

(5.10)

The classical infinitesimal strain tensor ε therein is the symmetric part of the displacement gradient, with three independent components in the present 2D situation. The second gradient of the displacement field, tensor K , is symmetric in the last two indices and it has thus six independent components, due to the symmetries ,

ij ji ijk ikj K K    . In 2D, the displacement field is the vector   , T uv  u
, which entails the following strain tensor in vector format 1 ,, 2

T u v u v x y y x                ε (5.11)
Similarly, the strain gradient tensor writes in vector format from the second gradient of the displacement field as

2 2 2 2 2 2 2 2 2 2 , , , , , u v u v u v x y y x x y x y                      K (5.12)
In the strain-gradient theory of linear elasticity, the constitutive law involves the symmetric Cauchy stress tensor and the hyper-stress (or double stress) tensor s . The tensors and s are related to the strain tensor ε and the strain gradient    K ε through the following general constitutive law for a homogeneous anisotropic second order grade continuum: The stress and hyper-stress for the effective 2D second order grade continuum can be defined in vector format as

               s (5.15)
The main purpose of this section is to determine the effective constitutive coefficients of the strain-gradient continuum from the response of random fiber networks; the short cut SG is conveniently used in the sequel to denote the employed effective second gradient continuum.

We next design a set of boundary conditions for the sequential determination of the components of the constitutive (rigidity) constants over various 2D domains  with boundary  , as explained into ore details in [20].

In each case, we force the WOA to bear a set of specific deformation and deformation gradients, as detailed in [20], and compute numerically the total elastic strain energy WOA U stored in the window of analysis under the corresponding boundary conditions. Similar to the identification procedure for the couple stress continuum (section 2), the first and second order gradient elastic moduli are computed by equating the total strain energy stored in the RVE with the energy of an equivalent homogeneous strain-gradient continuum, viz

2 WOA strain gradient ij ijkl kl ijk ijklmn lmn V U U A K D K          (5.16)
The left-hand side in (5.16) is the total elastic strain energy stored in the window of analysis, while the right-hand side is the expression of the energy of the postulated effective strain-gradient continuum.

The characteristic lengths are essential parameters for the second gradient continua; we generalize the definition of these parameters to an anisotropic continuum in terms of the engineering constants. In 3D, the six internal lengths associated to the independent classical moduli A can be identified by the expressions [19]:

1 2 3 1 rr r D l A           (5.17)
The previous equation can be simplified in the present 2D context, resulting in three internal lengths in tension and shear associated to the independent classical and second order elastic moduli as follows: the presently employed second gradient effective medium.

Acoustic properties of the second order gradient substitution continuum

The equations of motion for an effective second gradient continuum are easily obtained; they write in components form as the two following differential equations along the x and y directions of the Cartesian coordinates system

                               (5.
                              
(5.20)

The second order time derivatives u and v therein are the horizontal and vertical components of the acceleration vector. Relying on the plane harmonic wave solutions of the dynamical equations leads to a wave motion equation describing the propagation of longitudinal and shear waves, written in compact form as:

    1 1 2 2 4 4 2 2 * 2 A k + A k D k D k +D k k -ρ ω k k βA +A 11 33 111 122 112 1 2 12 33 1 2 1 2 2 1 0 2 2 1 1 0 2 2 4 4 2 2 * 2 k k 2A +A A k + A k D k D k +D k k -ρω 1 2 12 33 22 33 211 222 212 2 1 1 2 2 1 2 2 U V                       (5.21)
Nontrivial solutions of the last equation exist if the determinant  

12

,, w k k  of the matrix on the left-hand side Eq. (5.21) vanishes, which shall provide the dispersion relations. We shall next investigate the effect of internal bending length b l , fiber density and WOA size on the dispersion relation, and compare them with those obtained for the couple stress effective medium.

Dispersion relations and phase velocity for the second order effective continuum versus internal bending length l b . Comparison with the couple stress effective medium

We evaluate in Fig. 5.14 the frequency for the homogenized second gradient medium versus the internal bending length b l , and compare with the frequency evaluated for the couple stress effective medium. Inspection of Fig. 5.14 shows that the influence of internal bending length is the same for both couple stress and second gradient theories, the partial band gap increasing with b l for both models, while the existing small shifts for the L and S modes increase with b l .

Fig. 5.14 Frequency band structure versus fiber bending length, comparison between second gradient and couple stress theories, for wavenumber k=0.5

The main difference between the couple stress theory and the second order gradient substitution media is clearly visible looking at the polar plot of the modulus of the group and phase velocities (Fig. 5.15). .Orange line: shear mode; 0.5 kL  The dispersive behavior of the second gradient medium in the non-affine regime is evidenced by the modification of the shape of the plot when changing the wavenumber k. In the affine regime, when the macroscopic deformation coincides with the microscopic deformation, the influence of the second gradient disappears and no significant effect in the phase velocity is observed when increasing k, and the medium behaves as a Cauchy continuum.

Effect of density and window size on the dispersion relation, phase and group velocities for the second gradient medium

We next investigate the effect of the density and window size on the frequency band structure for the second gradient medium, and make a comparison with the couple stress theory in both affine and non-affine regimes. The same effect of the density and WOA size on wave dispersion can also be seen on Fig. 5.17; the anisotropy of wave propagation does not change when changing density or WOA size. The shape of the phase velocity plot changes with the wavenumber k, thereby reflecting the dispersive behavior of the second gradient medium, for any values of window size, density and in both affine and non-affine regimes.

Conclusion

The dynamic analysis of random fibrous networks is a novel topic, which has many applications in materials science or in the biomedical area, since these networks are very often subjected to dynamical loading conditions. In order to bypass the complexity of performing dynamic computations at the microscopic scale of the full network, we herewith develop and identify couple stress and second gradient models as effective continua at the mesoscopic level of windows of analysis of different sizes, in order to investigate the size effects of such networks on their dynamic properties. The static mechanical properties which are at the basis of the dynamical analysis are computed thanks to FE simulations performed over windows of analysis subjected to mixed boundary conditions allowing to capture the classical and non-classical effective moduli.

The acoustic properties of these networks are captured by the dispersion diagrams and plots of the phase and group velocities; we analyze the influence on the dynamic properties of three main quantities of interest, namely the fiber bending length, the size of the window of analysis, and the fibers density. The impact of these parameters is successively assessed for the couple stress and strain gradient substitution continua.

The influence of internal bending length is the same in both effective theories, the partial band gap increasing when increasing the internal bending length b l . Density and window size have the same influence on the dispersion relation for both models and in both affine and non-affine regimes. A significant increase in frequency for the longitudinal and shear modes occurs when density or window size increases when moving from couple stress to second gradient medium.

The same effect of density and size of the window of analysis on the dispersion diagram is also obtained. The anisotropic feature of wave propagation does not change when changing the density or the size of window of analysis. The main benefit of using random fibrous networks for wave propagation is the possible control of band gap by different parameters, like fiber density, window size, and fiber length. In the affine regime, when the macroscopic deformation is very close to the microscopic deformation, the influence of the microstructure disappears and no significant effect of the microstructure on the phase velocity occurs, thus Cauchy elasticity is sufficient.

Although the replacement of the initial random network by effective generalized continua is able to capture size effects and is more tractable from a numerical point of view, it remains clear that it loses some of the details of the microstructure response when waves propagate. Accordingly, such mesoscopic models shall be complemented by fully resolved local analyses at the scale of the fibers using Bloch's theorem or suitable extensions of it to better isolate individual phenomena associated to the mechanical response of the fibers themselves.

The influence of large strains developed by random fibrous networks on wave propagation is an important aspect encountered in real situations that shall be investigated in future contributions.

II. Analysis of nonlinear wave propagation in network materials

Chapter 6: Wave propagation in pre-deformed periodic network materials Summary

In this chapter we explore the influence of large deformations on the propagation of acoustic waves in repetitive network materials. The problem of elastic wave propagation in pre-deformed elastic materials and structures is highly interesting in many applications. The development of finite deformations within the network induces a preferred anisotropy associated to privileged directions into the structure, and leads to a continuous modification of the effective mechanical properties, which in turn modifies the wave speed and the width of the existing band gaps. The analysis of incremental motions superimposed on large state of deformation uncovers complex physical mechanisms which require sophisticated methods. Both theoretical and numerical methods are developed in this chapter in order to assess the influence of finite strains developing within repetitive networks on the evolution of their band diagrams. An incremental scheme for the update of frequency and phase velocity of the computed homogenized medium has been developed successively considering 1D, 2Dand 3D structures; it incorporates an update of the frequency and phase velocity of the propagating waves versus the effective density and the state of finite deformation of the effective continuum used as a substitution medium for the initial repetitive network. The applied deformation is shown to have significant effects on the wave frequency and phase velocity. The influence of the effective density on the dispersion relation and band diagrams under the application of an incremental deformation over the lattice unit cell is shown. The effect of the poison ratio on the dispersion relation is well presented and discussed in this chapter for 2D and 3D structures. One originality advocated in the present work is the derivation of the perturbation method in a 2D context for nonlinear periodic structures covering geometrical and material nonlinearities; this extension allows a validation of the dispersion relation obtained for the nonlinear effective continuum in the low frequency range when neglecting the effect of material nonlinearity. The methodologies and trends obtained in this chapter provide insight and guidance on selecting the topology of network materials which may give rise to large amplitude-dependent band gap shifts and phase velocity variations. Amongst these periodic structures, lattice materials are characterized by a regular and periodic microstructure that can be idealized as a network of beams or rods; they offer many advantages from a structural point of view, since they combine low weight, high stiffness and strength, and high energy absorbing capabilities that cannot be realized using homogeneous materials [2][3][4][5].

Furthermore, the topology and microstructural stiffness of lattice materials can be adjusted to fulfill specific requirements like a controlled state of anisotropy and density, or predetermined buckling and collapse modes. The continued interest for developing lattice materials has been accompanied by recent advances in manufacturing techniques like 3D additive printing, allowing the production of lattice materials with perfect control of their microstructure, and at a low cost.

Such technologies allow exploring a quasi-infinite range of possible microstructures, and have raised the need to develop predictive micromechanical models for the analysis of complex components made of lattice materials, to avoid repetitive testing.

The literature devoted to the modeling of such lattice materials most of the time restricts to the geometrically linear regime; in many applications however, the design of bending dominated lattices exploited for morphing structures [6][7][8] requires models of lattice materials in the nonlinear regime. We presently develop the discrete homogenization method for network materials consisting of beam elements advanced in [9], which shall be extended to the nonlinear regime in the present contribution. We assume due to the small bending stiffness of the structural beam elements that the nonlinear lattice response is essentially due to the change of the network configuration (the orientation and length of the beams change with ongoing deformation). The geometrical nonlinear behavior of cellular structures and network materials was extensively studied in [10][11], considering the example of foams, using simplified pin jointed model for which the bending contribution of the skeleton struts was neglected. The authors in [12] proposed another approach where axial, bending and twisting deformations at local level were considered.

More recently, Janus-Michalska [13] extended the linear model developed in [14][15] to construct the stress-strain relation and strain energy function for hyperelastic cellular materials with arbitrary symmetry. An alternative approach was proposed by [16] using a computational homogenization method to derive a nonlinear constitutive model for lattice materials.

Many techniques have been developed to predict the mechanical properties of heterogeneous structures, and especially lattice materials, thereby bypassing the need to resolve the smallest spatial scales. Homogenization methods aim at describing the overall response of heterogeneous structures including composites and periodic structural lattices in terms of effective properties, as presented in the recent contribution [16] and references therein.

The incorporation of nonlinear aspects of wave propagation in structures is necessary whenever large deformations occur [17][18][19][20], but it remains a considerable challenge. Two types of nonlinearities may be present in a broad sense, which can be classified as material nonlinearities [21][22] and geometrical nonlinearities [23]. This last type of nonlinearity is related to the evolution of the microstructure or structure configuration, for instance the change of configuration of a repetitive network, and it can be modeled as a succession of incremental deformations associated to the modification of the structure geometry [23].The presence of a nonlinearity in periodic structures results in amplitude-wave dependency in the dispersion relations.; this opens new possibilities for a passive tuning of the dispersion band structure through an amplitude-dependency of propagating waves, thereby going beyond a mere control of the dynamic and acoustic properties of repetitive structures by the design [24] or by application of an external stimulus.

A nonlinear periodic structure supports a variety of wave solutions depending on wave amplitude, waves interactions, and type of nonlinearity, for example solitary wave solutions for Boussinesq type equation, harmonic plane wave and discrete breathers [25][26][27][28][29][30][31]. In [32] the authors analyze the dispersion relations for axial and flexural elastic wave motion in homogeneous beams subjected to finite strains. The authors in [33] study wave dispersion in a one-dimensional nonlinear elastic metamaterial; the large elastic deformation provide the nonlinearity in the thin rod, whereas the metamaterial behavior is associated with the dynamics of the local resonators. The authors in [34] developed a numerical analysis improved by experimental measures to show the evolution of the locally resonant band gap, under nonlinear pre-deformation. The authors in [35] demonstrate the ability to use deformation to transform phononic band gaps in periodic elastomeric structures. In [24], the authors focus on material which constitutive law contains cubic stress-strain nonlinearity; while the authors in [36] present an analytical and numerical method have been compared to describe the propagation of nonlinear wave within a structure endowed with a square stress-strain relationship.

The effect of pre-stress or pre-strain on wave propagation through homogeneous anisotropic media has raised the interest of many authors; the effect of preexisting finite elastic deformations on wave propagation has been analyzed in [37,18] and the effects of incremental deformations on a homogeneous continuum medium has been studied in [38][39]. The initial deformation must be large enough to change the geometry of the medium, since an infinitesimal initial deformation would not affect the properties of the material based on the superposition principle valid for small deformations. The incremental effective properties of pre-stressed homogeneous media undergoing large deformation have been analyzed in [40,41], wherein the authors put some restrictions on strain energy sufficiently enough to allow elastic waves to propagate within the material.

In [42][43][44][45][46], the authors analyze the propagation of waves in composites consisting of a small number of layers (two or three layers) undergoing sufficiently large deformations.

It is worth mentioning that studies of nonlinear wave propagation in structures essentially deal with one-dimensional (1D) systems, whereas the nonlinear wave dynamics in multi-dimensional (2D) discrete systems has not been thoroughly investigated so far.

In this chapter, we analyze the dynamical properties of periodic network materials subjected to finite strains, relying on dedicated homogenization techniques developed to substitute to the initial discrete periodic lattice an effective Cauchy continuum (relying on Bernoulli beams). An incremental scheme for the update of network geometry, mechanical response and frequency is set up successively in 1D and 2D situations, based on the effective nonlinear medium obtained by homogenization.

Microscopic and mesoscopic nonlinear homogenization problems

In this work, we rely on the discrete asymptotic homogenization method -abbreviated as DH method in the sequel -for beam lattices initially proposed in [9,[47][48], which is utilized to construct Cauchy effective continuum models of repetitive 2D lattices endowed with translation degrees of freedom. The homogenization method relies on the connectivity of the studied network materials and the mechanical equilibrium at each node; it allows the treatment of elementary cells representative of the entire network, including internal nodes.

Generally speaking, the discrete homogenization can be described as a mathematical method to derive the equivalent continuous medium behavior of a repetitive discrete structure made of elementary cells. This technique is inspired by the homogenization of periodic media developed since the early eighties [49][50], with continued works in a more recent period [51,52].More recently, Pradel and Sab in [53] applied the discrete homogenization in combination with the energy method. The interest of the discrete homogenization method is that it delivers the full compliance (or rigidity) matrix of periodic networks, reflecting the sometimes complex anisotropy of the equivalent continuum.

The general idea at the base of the method is the periodic repetition of an elementary cell made of beams connected at nodes to define an infinite lattice; it may be explained as follows [9].

Consider a finite 2D (surface)or 3D structure, parameterized by a small parameter, defined as the ratio between characteristic lengths of the lattice unit cell to a characteristic length of the entire network, scalar quantity L (Fig. 6.1). For a large enough lattices, the ratio of the beam length, the scalar The statics for beams obeying Euler-Bernoulli kinematics (Fig. 6.2) is standard, and leads to the expressions of the normal and transverse forces given in the Chapter 1.

Microscopic incremental problem over the lattice

The analysis is made over the RVE selected as the unit cell, relying on recent results in [54] showing that the RVE size has no influence on the predicted homogenized response in the nonlinear regime, as long as no bifurcation occurs. The developed DH method leads to an algorithm for the computation of the large strains response of the considered networks. The method has been exposed into details in [47,48], thus we shall summarize below the main steps leading to the construction of the effective nonlinear Cauchy continuum.

In a first step, the homogenized constitutive law is evaluated in the linear framework, as previously exposed. As a next step, one set up a kinematic driven scheme through which, for each load increment, the incremental stress tensor is computed versus the imposed mesoscopic loading over the unit cell (transformation gradient).

We write down the non-linear equilibrium problem associated to the large perturbations of the network; the nonlinearity is due to the large displacements, which in turn are responsible for the large changes of beam directors orientation and beam lengths. The kinematic quantities that vary during the incremental scheme are the beam length , both pictured in Fig. 6.27 in the Appendix A. A kinematically driven scheme is written, based on the incremental equilibrium equations in translation and rotation, successively the two following equations 

                                    RR RR RR b k b k b k b k kk 3 b b R , F F . F . x x x x x x bb b k b k b k b k b k b k F F . F . y y y y y y bb b k b k b k b k b k b k F F . F . 0, z z z z z z bb M 3 R ,                                                                     
                                                    R R RR E b k b k E b k b k E b k b k Eb M M . x x y y z z O b k b k O b k b k O b k b k Ob b M M M . x x y y z z Cb bb L. x b E b k E b O b k O b C b bb M . M . l . 0 x bb                                                              
(E (b)) O (b) R R    v v v
the relative virtual velocity of the two extremity nodes of any beam, and w a virtual rotation velocity.

For any beam within the reference unit cell R bB  , the beam vector can be evaluated from the imposed transformation gradient j G and the relative position vector between the two extremity nodes, as

    R R E b O b b j jb     B R R G (6.2.2)
The following notations have been introduced:

o j j     R G
is the imposed mesoscopic transformation gradient over the unit cell, and

    RR E b O b  RR
is the unknown kinematic relative position which is computed incrementally. The variation of the beam length and orientation as well as those of the static forces and moments has been evaluated in [47,48]. The solution for this differenceis obtained by solving iteratively the previously written incremental equilibrium problem posed over the RUC, based on the modified Newton-Raphson method.

This incremental scheme provides an update for the kinematic variables b n1  B at any time step denoted by the integer subscript   n1  , from their values at previous iteration k within the reference unit cell: for a given imposed mesoscopic transformation gradient, the update of the position vector is written as

b b b n n n1     B B B
with from (2.8)

    R R E b O b b j jb n n n     B R R G (6.2.3)
For a given   and for each imposed transformation gradient , 1, 2,3,

o jb j j     

R

one has to determine the relative position vector

    R R E b O b n n  R R
in (2.8) at each time step, which entails the evaluation of quantities ,,,,,,E written at the continuum level of the RUC.

      3 ,

Mesoscopic equilibrium equations over the reference unit cell

The discrete incremental equilibrium takes after homogenization and condensation of the nodal rotations (by using the equilibrium of moments) a form similar to its continuum counterpart, viz

i i R R b b i(k) b b B 0 d0 i                             B F v M W S v (6.2.4) with i(k)
S the incremental stress vector; previous integral formulation constitutes the incremental weak form of the equilibrium (self-equilibrium) posed over the Lagrangian domain  occupied by the reference unit cell.

We next write the incremental constitutive law at the mesoscopic level. Strain-controlled loadings are imposed over the unit cell, so one imposes at each increment the discretized version of the transformation gradient F together with its increment F , specific to each type of loading. In view of setting up an incremental scheme for the resolution of the nonlinear mesoscopic BVP (acronym for boundary value problem), we write the discretized version of the elastic constitutive equation relating the incremental stress

  k n
S to the incremental strain   k Gn E (the subscript k is an iteration counter inside the increment loop)

      k k k S n T, n Gn n 1 n n n n k K : ,             S E S S S S S (6.2.5)
The Green-Lagrange strain G E therein is defined by

  1 . 2  T Ε F F -Ι G .
For a given elementary loading characterized by a continuously varying scalar loading parameter n  (incorporated into the imposed transformation gradient), the incremental Lagrangian strain Gn E is related to the increment of the loading parameter n  and transformation gradient n F as [55]:

            Gn n n n n n n n n n , sym .grad sym .           E F u F F (6.2.6)
together with the elaboration of the incremental transformation gradient

    n n n n n n n I grad grad grad           F u F u u (6.2.7)
The algorithmic material tangent stiffness matrix S T,n K therein results from the assembly of the microscopic tangent stiffness matrices, as exposed in the Appendix A; it is obtained at each increment as the sum of three contributions given in the Appendix A

  R S b b b T,n o,n u,n ,n bB K K K K       (6.2.8)
The updated Kirchhoff stress tensor is next obtained by a push-forward of the incremental Lagrangian stress

  k n S from configuration n  (at increment n) to n1   (at increment n+1), elaborated as             ( k ) k n n kk 1 (k 1) T 1 (k) T 1 T n 1 n n n n n n n n n 1 n 1 n n 1 J . . J . . J . .              F S F F S F F S F (6.2.9)
where   n n J : det  F is the Jacobean, defined as the determinant of the deformation gradient tensor. Kirchhoff stress at increment n1  is accordingly given from its counterpart at previous increment n , tensor   k n1  , based on the initial stress (k ) n S (the first contribution on the right-hand side of (6.2.8)) and the incremental stress

  k n



(the second term on the right-hand side of (6.2.8)). Based on (6.2.9), Cauchy stress is updated over the imposed incremental loading as (the converged values at both times steps are considered)

Effective incremental frequency and phase velocity of a 1D microstructured beam

The dynamical analysis under large strains is first performed in a 1D context, a situation that can be illustrated by a macroscopic beam including a repetitive microstructure.

In order to set the stage, we rewrite the incremental nonlinear constitutive law (Fig. the corrected incremental Young modulus accounting for the nonlinear correction to the linear term nn Ee  in (6.3.1), itself depending on the linearized strain increment. Previous incremental constitutive law describes the beam response to an imposed increment of deformation gradient, itself equal to the small strains increment,

        n n n F x gradu x e x      .
The incremental nonlinear elastic constitutive law written in (6.3.1) is next applied to the analysis of wave propagation through a pre-strained nonlinear microstructured beam. Let consider a beam incorporating many unit cells repeated by periodicity along the longitudinal direction, responding in a nonlinear manner to an imposed incremental strain. All fields in (6.3.1) depend upon the variable x which is the beam axial coordinate; periodicity of the microstructure along x implies that the beam is considered as macro-homogeneous so that the fields experience a smooth variation with x .

Omitting index n, the dynamical incremental equilibrium equation for the continuous displacement of the homogenized continuum writes based on (6.

                      (6.3.4)
Substituting expressions (6.3.4) into the weakly nonlinear wave equation (6.3.3) and ordering versus the successive powers of the small parameter  produces a set of equation as follows

          2 0 2 0 2 0 0 0 11 22 2 1 2 1 2 1 0 0 11 22 2 0 2 0 2 0 2 0 0 1 1 0 1 1 1 1 1 2 2 2
:0 :

2 nn n n n nn n n n t n n n n n n n n n uu OE x uu OE x u u u E x                                             (6.3.5)
The solution of the   0 O  term equation is well-known and is given by planar harmonic waves of frequency 0

1 n   , ( 0 0 1 1 , n n
   is the frequency and the density for increment n respectively, and which will later be denoted ,

nn  ) viz     0 exp cos n u A i kx A kx       (6.3.6)
in which k is the wavenumber and A the amplitude. Subsequent substitution of the expression (6.3.6) into the  

1 O  term results in the equation       2 1 2 1 2 * 2 1 1 2 11 2 2 2 cos t n n n n n n n n n n n uu E E k A kx x                      (6.3.7)
The linear kernel of order  

  2 1 1 2 11 2 0 t n n n n n n Ek           (6.3.8)
Thus, the frequency for the new structure configuration is updated versus the actual tangent modulus and density by the relation

t n n n n n n n Ek              (6.3.9)
The quantities , nn  therein are the effective density and angular frequency of the medium at increment n, according to the  

0 O  equation, 1 n
  is the density for the new configuration and k the wavenumber. The effective incremental phase velocity is given by the following relation versus the updated density and tangent modulus:

1 1 3 22 t n n n n n n n cE c c           (6.3.10)
Parnell [23] considers the simple case of a composite beam, which has experienced some nonlinear deformation; it then behaves incrementally as a linear material; thus, for each increment, a linear dynamical problem is considered and the authors do not adopt an incremental update for the frequency as in the present contribution (reflected by Eq. (6.3.9)). The effective incremental phase velocity of the beam is given by [23]:

* * * p i i E c   (6.3.11) in which * 0 i t in n EE  
 is the sum of all tangent moduli [23]. The extension of this method to 2D and 3D structures is possible from a theoretical point of view, but it requires most of the time lengthy computations; furthermore, the method used by the last author requires calculating the effective density for all increments in order to compute the medium frequency at a given increment. In the following the units of wavenumber is m -1 . The nonlinear phase velocity response of a beam incorporating many hexagonal unit cells (pictured in Fig. 6.5a), repeated periodically along the longitudinal direction and subjected to a nonlinear tensile deformation xx E is computed based on previous incremental scheme (6.3.9), using Parnell method. A comparison of the phase velocity computed by the present incremental method and based on Parnell method is presented in Fig. 6.6, showing very good agreement between both methods, even up to high values of the imposed deformation gradient. One can notice from the values of the phase velocity an important difference between the linear and nonlinear situations under the imposition of a deformation before sending plane waves through the structure. We further observe in Fig. 6.7 that the imposed compression has a significant effect on the incremental phase velocity, which is much more pronounced in traction than in compression, due to the fact that Young modulus varies very little under compression (Fig. 6.7b). The strong variation of Young modulus under traction (Fig. 6.7b) counterbalances the increase of density (Fig. 6.7a) so that the phase velocity in turn increases under a tensile loading (Fig. 6.7c). These results entail overall that the dynamical behavior of the material will be modified by decreasing the wave velocity, especially by imposing a compressive strain before sending the wave into the structure. This feature is especially interesting in situations in which dissipation occurs, since compression tends (as dissipation) to increase wave absorption phenomena.

Fig. 6.7 Density (a), Young modulus (b) and phase velocity (c) variations versus the tensile and compressive strain for a beam incorporating a hexagonal repetitive microstructure.

We next analyze the impact of the unit cell geometry of the beam microstructure on the effective incremental phase velocity versus the imposed deformation xx E , itself built from the imposed transformation gradient xx F .We plot in Fig. 6.8 the evolution of the phase velocity for the hexagonal, re-entrant and Milton lattices shown in Fig. 6.5. We next extend the dynamical nonlinear analysis to planar waves for 2D network materials.

Incremental dispersion relation and phase velocity in a 2D homogenized medium

In a 2D context, the constitutive law governing the evolution of the homogenized network material subjected to large deformation gradients writes relying on the developments at the end of section 2 and as an extension of Eq. ( 6

                           (6.4.1)
where n c is the tangent stiffness matrix of the effective continuum at increment n introduced at the end of section 2 (Eq. 6.2.12 and 6.

ˆt nn n n n n n n div c c div tt                uu ε (6.4.2)
The tangent stiffness matrices , t nn cc therein are considered as constant between increments n and 1 n  , and * 1 n   is as in the 1D case the effective network density for increment 1 n  .

Using the same methodology as in the previous 1D context, and restricting to Centro-symmetrical structures for which the coupling coefficients vanish (the following components of the rigidity matrix vanish, 13 For a plane wave without attenuation in the x-y plane, the propagation wave number constants along the x and y directions are 1 k co k= s(θ) and 2 ksi k= n(θ) , with the angle  indicating the direction of wave propagation with respect to the x axis of the Cartesian basis. We shall note that

l l nn n l n n n a k a k a k a k a k a a k k a a k a a k k a k k a a k k a k a a k k a a k a k a k k                              24
k a k a k a k a k a a k k a a k a a k k a k k a a k k a k a a k k a a k                                         2 4 2 2 2 2 4 33 1 33 1 2 33 2 2 a k a k k a k                  (6.4.3)
We observe in the initial configuration (no applied deformation gradient) a maximum band gap about 14.46 rad/sec, which increases to 85.12 rad/sec when moving from the regular to the reentrant configuration (the hexagonal lattice becomes re-entrant when loaded under a compression strain of 0.5

yy E 

). This effect becomes more pronounced as the compression deformation level is increased. A similar increase is also observed under simple shear, although less pronounced.

The same effects are observed for the hexagonal lattice under an incremental shear load, since the band gap increases from 14.46 rad/sec to 54.2 rad/sec for a shear strain of 0.275 (expressed in percentage of the initial value). The comparison with the compression loading situation shows that the band gap is larger in comparison to shear; this can be assessed by comparing the thickness of the band gap for the hexagonal lattice under a compression strain equal to 0.25 and a shear strain equal to 0.275. The increase of the band gap width from the regular hexagon to the re-entrant configuration is characteristic of the special acoustics properties presented by auxetic metamaterials [57]. These results are in good agreements with [58], in which the authors compare the dynamical properties between the hexagonal and the re-entrant lattices.

In order to highlight the non-isotropic dynamical behavior of re-entrant hexagonal configurations and the impact of anisotropy on the dispersive characteristics, we plot the phase velocity for the hexagonal lattice submitted to a gradient of deformation corresponding to compression (Fig. 6.10) and shear (Fig. 6.11) successively. ). The degree of anisotropy becomes higher as the level of compression increases, due to the modification of the lattice geometry; the anisotropic behavior of the re-entrant lattice is markedly visible from the irregular shape of the phase velocity plot for the longitudinal and shears modes [58].

The study of wave propagation in auxetic materials is an important topic [59], and especially the evolution of the dispersion relations in the large deformations regime is quite interesting. We shall illustrate this aspect in the present work by studying the effective wave propagation within the hexachiral structure, the dynamical behavior of which is computed under incremental uniaxial and biaxial applied loads.

In order to investigate the evolution of the dispersion relation and the bandgap under an incremental deformation gradient, we plot in Fig. 6.12 the dispersion relation for the hexachiral network submitted to uniaxial and biaxial loads, for the longitudinal and shear modes. The remarkable previous effects of the imposed gradient of deformation on the dispersion relation are also observed (Fig. 6.12): an increase of the incremental load results in shifts in the dispersion relation, for both biaxial and uniaxial loads. It is also apparent from Fig. 6.12 that both types of loads have the same influence on the dispersion relation and band gaps for the longitudinal and shear propagation modes. The wave propagation modes for the longitudinal and shear modes are very close, as shown in Fig. 6.12; this can be attributed to the symmetry of the stiffness matrix for the hexachiral lattice (the responses in x and y direction are similar).

Fig. 6.12 Dispersion relation for the hexachiral lattice under biaxial loads (left) and uniaxial loads (right). The solid line corresponds to the longitudinal mode, the dashed line to the shear mode.

The phase velocity plot for the hexachiral lattice submitted to an equibiaxial deformation gradient (Fig. 6.13) highlight its isotropic behavior as well as the non-dispersive nature of wave propagation. The phase velocities for the hexachiral lattice obtained for the longitudinal and shear modes are very close to each other (Fig. 6.13). After introducing the non dimensionless time t   , we obtain the incremental frequency for each mode as follows: We plot in Fig. 6.15, the evolution of the dispersion relation for the hexagonal 3D submitted to an incremental load along the x axis. As expected, as the applied load is increased, a shift in the three modes of propagation occurs, and the width of partial band gap between the modes also increases (Fig. 6.15).

L L l a k a k a k a k a k a k a k a a k k a a k k a a k k a a k k a a k a a k k a k a a k k a a                            22 77
SV SV SV kk a a k a a k k a k a a k k a a k k a k a k a k a a k k a a k k a a k k a k a                                             2 2
k a k a k a k a k a k a a k k a a k k a a k k a a k k a a k a a k k a k a a k k a a k k a a k a a k k a k                   2
SH SH SH a a k k a a k k a k a k a k a a k k a a k k a a k k a k a k a k                                        (6.
The maximum width of the partial band gap occurs between the longitudinal mode and vertical shear mode increases from 28.93 rad/sec at the initial case to 95.02 rad/sec after an applied traction of 0.175 deformation.

In the sequel the red line corresponds to the longitudinal mode, the blue one corresponds to the vertical shear mode and the green one to the horizontal shear mode.

We can be clearly observed from the Fig. 6.16 the anisotropic behavior of the 3D non-regular hexagonal lattices from the variation of the components of phase velocity under the axial incremental deformation. anisotropic is more pronounces as we increase the level of the applied deformation.

Effective dynamical response of 3D auxetics

Fig. 6.17 shows a representative unit cell of the 3D re-entrant honeycomb and pyramid shaped configuration. The evolution of the frequency versus the incremental deformation gradient (traction in x direction) for the pyramid and the re-entrant 3D lattices is illustrated in Fig. 6.18 and Fig. 6.19, for two directions of wave propagation (χ is the angle between the wavevector and z axis while  the angle in the xy plane). Increase in the gradient of deformation entails an increase of the partial band gap between L and SV modes, SV and SH modes, in the direction of propagation χ=π/6ν Ө=π/6 as one can observe in Fig. 6. (18,19). The frequency increases when the applied gradient of deformation is increased.

The pyramid lattice shows very good in-plane mechanical properties but it suffers from some weakness in the out-of-plane direction and that can be clearly observed in the case of in-plane propagation (χ=π/β), how the longitudinal and vertical shear wave are very higher than the horizontal shear mode. However the weakness properties of the 3D re-entrant lattice under shear solicitation results two lower modes of propagation, vertical and horizontal shear modes as shown Fig. 6.19 for in plane propagation (χ=π/β).

It is interesting based on previous results to note that the dynamical behavior of such periodic materials can be improved by imposing an incremental load (gradient of deformation) before sending any wave through the structure, sue to the increase of the band gaps and decrease in the phase velocity, as shown previously on Fig. 6.8.

Effect of Poisson's ratio and density on the frequency band diagram

We investigate the effect of Poisson's ratio and density on the frequency band structure for the 2D hexagon under uniaxial compression and the 3D pyramid under a uniaxial tensile loading.

Note that the variation of Poisson's ratio is presently obtained from the deformation of the structure under the imposed kinematic loading (Fig. 6.22). The effective incremental frequency of such pre-strained structures are applicable only in the low frequency domain, when the wavelength of the propagated wave is much greater than the characteristic length of the representative volume elementary (RVE). At higher frequencies, the effective medium suitable for low frequencies is no more representative, thus one shall have recourse to perturbation method. The extension of the perturbation method to incorporate the effect of a large pre-deformation of the structure is the object of the next section.

Nonlinear dispersion relation based on the perturbation method

Thus far we have only considered the behavior of the material at low frequencies, a domain in which it makes sense to have recourse to the homogenized effective medium. At higher frequencies, the effective homogenized medium is no more representative of wave propagation phenomena within the lattice, thus one has to assess the range of validity of the effective medium via a suitable generalization of Bloch analysis to finite strains.

Bloch analysis has been solved so far in the linear setting, allowing the study of wave propagation in an infinite perfectly periodic linear medium, without the recourse to an effective medium approach, and which is suitable for a comparison of the dispersion analysis obtained by the asymptotic theory.

Each beam of the unit cell is considered as a Bernoulli beam including three degrees of freedom at each node, namely two translations and one rotation, condensed into the vector of kinematic

degrees of freedom   ,, i i zi uv   i q
. The continuum counterpart of these kinematic variables interpolating the discrete kinematics of a typical beam element is approximated by (x denotes the curvilinear abscissa along any beam):

    6 1 , ( ) nn n u x t a x q t        6 1 , ( ) nn n v x t b x q t        6 1 , ( )
z n n n x t c x q t     (6.7.1)
The scalar function (x), (x), ( ) nnn a b c x therein are the shape function.

The kinetic and potential energies per unit thickness of a Bernoulli beam taking into account the large deformation according to the Cauchy-Green tensor relation in the extensional mode only (without effect of inertia and torsion) is given by:

  66 .. 2 2 11 0 0 0 1 1 1 ( ) ( ) ( 2 2 2 L L L r s r s rs rs u v T S dx S dx S q q a a b b dx t t                 (6.7.2)
The potential energy is given by:

2 2 2 2 2 0 0 2 2 66 22 11 0 1 1 1 2 2 2 1 1 2 2 L L z L s s r r r r s z r s t rs Green lagrange relation linear part u u v U ES dx EI dx x x x ab a b a q q ES EI dx q q q ES x x x x                                             666 1 1 1 0 6666 1 1 1 1 0 1 8 L st r s t L s t u r r s t u r s t u nonlinear part nonlinear part aa dx x x x a a a a q q q q ES dx x x x x                                     (6.7.3)
in which , , , , L E Iz S are successively the length, Young modulus, second moment of area with respect to the z-axis and the cross-sectional area of a generic beam within the network unit cell.

Once the element shape functions are constructed, it is straightforward to calculate the linear and nonlinear stiffness and mass matrices after application of the Euler-Lagrange operator to write the nonlinear dynamical equations of motion.

The Euler-Lagrange equations of motion are obtained without external forces:

0, TU t q q           (6.7.4)
This equation of motion can be written for each beam; the assembled equation of motion for the unit cell is then derived in the following form:

2 3 0 q   p Mq + Kq + K q K q (6.7.5)
The matrix M and K are respectively the linear mass and stiffness matrix of the beam, and the nonlinear stiffness matrices ( p K , q K ) therein weakly contribute to the response of the system.

The quadratic terms p K do in (6.7.5) not contribute to the dispersion relation according to the first order perturbation theory, thus only the effect of the cubic terms q K on the system will be considered. The incremental equation of motion for the pre-stress medium (an incremental gradient of deformation is applied to the structure) is given by

    3 0 t t n n n n n         NL M M q K Κ q K q (6.7.6)
where , nn MKare the assembled mass and stiffness matrices at increment n; similar to the 1D situation treated previously, ,

t t t t n n n n    M M K K in (6.7.6
) are the assembled incremental mass and stiffness matrices respectively (this last matrix has been computed in section 2 based on discrete homogenization), depending on the form and magnitude of the gradient of deformation applied to the structure, which can be calculated in the same manner as in the initial configuration (no applied deformation), but for different lengths and directions of the lattice beams,   NL NL KK is the nonlinear stiffness matrix (cubic terms) arising from the nonlinear stress-strain relationship (we use in subsequent applications Cauchy-Green tensor in the extension mode only), and  is a small parameter used to quantify the degree of nonlinearity.

The first step in the analysis of the nonlinear dispersion relation is the introduction of the non dimensionless time t   , followed by the asymptotic expansion of the frequency and displacement

    2 0 1 0 2 1 0 3 2 , , 0 4 t t NL n n n n n n n n n AA                  Μ Μ K K q D (6.7.9)
Nontrivial solutions of (5.9) exists only for a vanishing determinant, viz

    10 , , 0 nn A    D ,
thus leading to a relation between

1 n
 and 0 n  .

In the case of a mere geometrical nonlinearity [24], the    is then given by:

2 1 0 3 4 2 i H ii i H ii A       NL 0 K Μ (6.7.11)
Finally, the frequency

i n
 is given by the relation

2 0 0 3 4 2
H i i i i nn iH n i i A       NL 0 K M (6.7.12)
In order to illustrate previous formulation, the effect of a pre-strain on a structure made of repetitive hexagonal cells (Fig. 6.5) undergoing an incremental gradient of deformation (compression solicitation) will be studied. The dispersion relation for the pre-strained structure has been analyzed in Fig. 6.24 without the effect of geometrical nonlinearities (thus we have neglected the term

  3 0 NL n n
Kq, as done in the homogenization method presented in section 2), considering the propagation of the longitudinal and shear waves within the structure along the edge of the first irreducible Brillouin zone. We compare for each propagation mode the results obtained by the homogenized theory (solid line) and the perturbation method (dashed line). The structure undergoes a compressive strain in the y-direction. In subsequent figures, the red line corresponds to the initial configuration, the green line to the configuration obtained by the application of 30% deformation in y-direction, while the blue line corresponds to the band diagram for an imposed deformation of 50% (in y-direction).

Fig. 6.24 Dispersion relation for the hexagonal lattice under compression for the longitudinal mode (left) and the shear mode (right) for the 2D hexagon. Comparison between perturbation method (dashed line) and homogenization theory (solid line) for three loading Fig. 6.24 highlights that the homogenized theory reproduces the perturbation method for the two modes of propagations (longitudinal and shear) along the edge of the first Brillouin zone. The two modes describe non-dispersive waves, since the frequency is directly proportional to the wavenumber. We can conclude that, in the low frequency range (for both longitudinal and shear modes), the dynamical properties of the pre-stressed material can be investigated directly from the homogenized effective medium. We next investigate the evolution of the band gap and of the geometrical nonlinearity on the dispersion relation; Fig. 6.26 highlights the weak effect of the geometrical nonlinearity on the dispersion relation, even for large values of the wave amplitude. We observe a small shift in the band gap (a partial band gap exists between the longitudinal and shear modes when moving from low to high amplitudes). In the initial non perturbed configuration, a drop of w= 1 rd/s in the band gap occurs for higher amplitudes. A new configuration is obtained after applying a compression deformation 30% yy E  and a drop of 3.8 rd/sec occurs in the band gap, while this drop decreases to 3.08 rd/sec after increasing the deformation up to 50% yy E  . In both cases of low and high amplitudes, Fig. 6.26 shows that the band gap between modes increases, when the regular hexagon is deformed to a re-entrant configuration for a strain level 50% yy E  . Longitudinal mode (solid lines) and shear mode (dashed line). Red line: initial configuration. Green line: 30% deformation gradient in y-direction. Blue line: gradient of deformation of 50% in y-direction.

For higher amplitudes, the frequency of the shear mode rises above the frequency of the longitudinal mode, due to the lower values of the initial frequency in the shear mode.

8 Conclusion

We analyze in this chapter the influence of large deformations on the propagation of acoustic waves in repetitive network materials. An incremental scheme for the update of the frequency and phase velocity of the homogenized medium has been developed, considering successively 1D, 2D and 3D repetitive network materials. The constructed scheme relies on the computation of the effective tangent stiffness accounting for the variation of the lattice geometry due to the kinematic loading imposed over the unit cell. Starting from the writing of the dynamical equilibrium equation for a Cauchy continuum, we write in the spirit of Lindstedt-Poincaré method the asymptotic expansions of the frequency, effective density and displacement components, leading to update formula for the frequency and phase velocity of waves In the present large strains regime, since the beam length is changing, one has to expand it versus the asymptotic parameter  as all other kinematic variables (these expansions are not repeated in

this subsection), b b b 2 b p b 0 1 2 p l l l l l         (A3)
The induced perturbation of the resulting efforts and moments is then obtained as

                      O(b) b b E(b) sz M . . 2 . . b b dl l E M                                                                        B B
     bb b E(b) sz . . . 2 . b b dl l                             B B B I C A B e e φ e e φ (A4)
Insertion of the expressions (A1) and (A4) together into the incremental equilibrium (6.2.9) entails the identification of the total tangent stiffness matrix

  R S b b b T,n o,n u,n ,n bB K K K K       (A5)
with ou K , K , K  therein respectively the linear stiffness, the initial displacement stiffness and initial stress stiffness expressed in closed form in the sequel.

The linear stiffness matrix writes

          b b b b b b b b S s z s z o,n 3 3 3 3 b b b E A 12E I 12E I K l l l                         e e e e e e (A6)
The initial displacement stiffness receives the expression . l

                                      b b S u,n o 2 b b b b b oz b sz bb 3 b o bb 2 b b b b 3o b sz 3 b . E A 1 K . . Ll l 1 .. l 12E I B. l . l 1 . l 12E I l                                                                   
                                 Ω C.P Be e B C.P (A7)
The stress stiffness receives the expression

          b b b b b b ,n z 2 b b b N 2T T K . . ll l                                         C.P e B C.P Ω C.P (A8)

Overview

Periodic cellular networks have been the object of extensive investigations. From the early sixties, a growing attention has been paid to the study of the mechanical properties of periodic networks due to the great advantage (static and dynamic) presented from their spatial periodicity [1][2][3][4] for instance, the presence of frequency bandgaps.

A lot of attention has been paid to wave bandgaps in periodic networks [5]; within some frequency interval, waves cannot propagate through the network, but are attenuated quickly [5].

The location and width of the existing bandgaps depend on the network topology and on the direction of wave propagation. The propagation of elastic waves in a linear framework has deserved a wide body of works [6][7][8][9][10][11], whereas only a few authors analyzed so far wave propagation in nonlinear media. The incorporation of nonlinear aspects of wave propagation in structures is necessary whenever large deformations occur [12][13][14][15], but it remains a considerable challenge. Two types of nonlinearities may be present in a broad sense, which can be classified as material nonlinearities [16][17] and geometrical nonlinearities [18]. The propagation of nonlinear waves in periodic structures is accompanied by a number of new phenomena that are different and can never be observed for linear media.

The presence of a nonlinearity in periodic structures results in the dependency of the wave propagation, phase and group velocities upon the amplitude of the wave; this phenomenon deserves the name amplitude dependent dispersion relation. It opens new possibilities for a passive tuning of the dispersion band structure through an amplitude-dependency of propagating waves, thereby going beyond a mere control of the dynamic and acoustic properties of repetitive structures by the design [17] or by application of an external stimulus.

This entails that solutions of the wave propagation equations are much more complex compared to harmonic plane solutions of the linear wave equations: nonlinear periodic structures support a variety of wave solutions depending on wave amplitude, waves interactions, and the type of nonlinearity, one can for instance mention solitary wave solutions for Boussinesq type equation, and shallow water waves for Burger's equation [19][20][21][22][23][24][25].

The study of nonlinear elastic waves has been limited so far in the literature to classical Cauchytype elasticity theory, which relies on the sole first order displacement gradient. Models based on Cauchy-type theory do not give realistic predictions on the properties of the medium such as the dispersion relation, since the Cauchy effective medium lacks internal length parameters. The Cauchy medium is non-dispersive, which means that waves propagate independently of the wavenumber [26]. It has however been proven by experiments that most waves are dispersive, that is, each wavenumber travels with a different phase velocity [27][28]. This explains the success of gradient-enriched theories in capturing dynamic behaviors overlooked by classical elasticity. In order to circumvent this drawback, an energetic method has been developed in this chapter to link the dispersive aspects of wave propagation to gradient elasticity theories in a nonlinear effective medium obtained by the homogenization of the underlying network microstructure. We shall focus in this chapter to the hexagonal lattice.

T GS L L                     U e U e (7.2)
. . 1

                   U e U e b b c z c M EI L L L   (7.3)
In those expressions, S is the beam cross section, L the beam length, I z the quadratic moment of the beam, . between both beam extremities.

Those expressions are written in a general case and exhibit a nonlinear elastic dependency of the forces and moment with respect to the kinematic variables. In order to apply the asymptotic homogenization method, some simplifications are made. We consider that all beams are initially rectilinear, so that no initial deformations and initial curvatures are present, implying that the corresponding kinematic variable vanishes in the expression of the normal force (the coupling between tension and bending is accordingly neglected).

Based on these general expressions of the forces and moments, the expression of the internal deformation energy of a single beam element is obtained, written in the local coordinates ( , )

 attached to the beam as Based on [30] and using Eq. ( 7.4), one can evaluate the continuous energy density of a hexagonal lattice made of extensible beams pictured in Fig. , ll KK and 3 l K . We further make the assumption that the three beams have equal rigidities for each deformation mode.

                      2 1 2 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 3 1 1 1 1 1 2 2 2 2 2 1 12 1 2 2 2                                                                           i i j j i i i i i i i i j j j j j j j j i i i i j j j j i i i i i i i i j j j j j j j j ES L EI L L L L W L EI L L                                     (
Expanding the finite differences in (7.4) up to the second order gradient of the continuous displacement field leads to a nonlinear second order effective continuum. After development, the energy density of the homogenized hexagonal lattice can be written in closed form versus of the slenderness ratio  and the angle  as (see the Appendix A):

                                1 3 2
E W L                                                    2 4 4 2
                                                        (7.5)
We next analyze one-dimensional and two-dimensional wave propagation within such repetitive networks, based on the constructed effective continuum based on (7.5), as will be done in the next section.

Dynamical equilibrium equation

For extensible beams, developments of Appendix A show that the energy density of a nonlinear second gradient medium takes the general form: The second order hyperstress follows similarly as:

  2 2 2         ijk ij k Ws S u x (7.8)
The equations of motion for a second gradient medium along the j x directions then write in index form as follows:

2 ij ijk j j j k * u , i,j, k S x x x 1, 2               (7.9)
The effective density therein is given in general by * From Eq. (7.9), we obtain two differential equations that describe the propagation of longitudinal waves polarized in the direction of incident wave, and of shear waves polarized in a direction perpendicular to the direction of the incident wave.

Throughout this work, one considers the material behavior at low frequencies, so equivalently when the wavelength of the wave is greater than the unit cell size, a range of frequencies in which it makes sense to have recourse to the homogenized effective medium. In order to show the validity of the homogenization method at low frequency, we compare the results obtained in the linear case to those given by Bloch theorem (for more detail see [10]). From expression (7.5) of the continuum strain energy, we obtain the specific expressions of the first order Piola-Kirchhoff stress and hyperstress tensors, based on the general definitions introduced in Eq. (7.7) and (7.8). Inserting these expressions into the equation of motion (7.9), we obtain the homogenized nonlinear wave equation:

2 2 4 2 * 1 2 3 2 2 4 2 u u u u u E E E x x x x t               (7.10)
where u is the longitudinal displacement, 1

E the linear effective modulus derived from the linear part of the energy density, 2 E the nonlinear effective modulus characteristic of the large deformation behavior and 3 E the second order effective modulus, characterizing the second order gradient behavior. These moduli are expressed versus the microstructural lattice parameters in Appendix B. The equation of motion (7.10) describes the quasi-static (dispersive) behavior of waves propagating within the periodic lattice. Let note that the term of the second order modulus in (7.10) is always negative; this means that the presence of a microstructure results in positive dispersion, so that increasing the wavenumber leads to increase the phase velocity.

Wave equation (7.10) represents a Boussinesq-type equation; an analogous mathematical model has been essentially formulated and studied for shallow water waves, to describe wave propagation in fluid [31]. The general properties of Boussinesq and Boussinesq-type equations have been intensively studied, including nonlinear effects in the dynamics of solids [32][33];

solutions of equations of this type have been studied in [34][35][36][37][38][39].

We shall note that Eq. (7.10) can be integrated using elliptic functions; single elliptical wave solutions for this equation are well known in the literature [40]. A set of elliptic functions depending on the degree of nonlinearity are presented explicitly in [41]; they are described by a universal coefficient s which does not depend on the mechanical properties of the material, but takes into account the shape, period and velocity of waves. Parameter s can be considered as a quantitative measure of how much the nonlinearity mode differs from the linear situation.

We shall next derive an explicit analytical stationary solution of the governing wave equation (7.10).

Using the change of variable z x t k

 

, where  is the frequency and k the wavenumber; we transform equation (7.10) into the ordinary differential equation for the new non-dimensional strain of the wave function [41]    u N z

: 2 2 2 0 N a N b N c z       (7.11)
where

2 2 2 3 3 1 , 2 l E E k a b E E     
and c is an integration constant.

The initial condition writes as in [41]:

, 0 for 0 2

     AN N z z (7.12) 1 0 0 N dz   (7.13)
with A the strain amplitude. The physical meaning of Eq. (7.13) is that the displacement cannot grow infinitively and this condition leads to the determination of constant c.

The solution of Eq. (7.11) can then be expressed as [41]:

  2 2 0 ( ) , () 22 21
()

k A A s N z sn z s Es Ks        (7. 14
)
where s is the so-called universal constant, the modulus of the elliptic functions, describing the degree of nonlinearity ( 01 s ). This parameter accounts for the period, shape, and velocity of the waves. Generally, it may be considered as a quantitative measure of the intensity of the nonlinear effects and of how much the nonlinear motion differs from the linear motion. For low values of s, there is a small effect of the nonlinearity on the constitutive law and thus the structure behaves approximately in a linear manner. As s increases, the influence of the nonlinearity on the constitutive law increases. In (7.14), sn (.) is the elliptic Jacobin sine, K(s), E(s) are the complete elliptic integrals of first and second kind respectively and 0 k is the propagation constant related to the wavenumber k as follows:

  0 2 kk Ks   (7.15)
The strain amplitude can be calculated from the following equation

  2 0 () 31 () bA Es Ks k    (7.16)
In the limiting case 0 s  (which corresponds to the linear case), the solution of (7.14) is reduced to the harmonic plane wave solution,

      cos sin 2 2 A A N z kz u kz k      with 0 kk 
In the opposite case, that is when 1 s  , Eq. (7.14) describes a localized solitary strain wave with the non-dimensional strain given by

        2 0 0 0 sech tanh 2      A A N z k z u z k z k .
Compression solitary waves can exist when A<0, and for A>0, dilatation solitary waves (tension wave) will be obtained.

It may be possible to develop an asymptotic solution for the nonlinear wave equation (7.10) in the case of a weak nonlinearity, corresponding to low values of s. The asymptotic development based on the Lidstedt-Poincaré method is explained in detail in [42], and the comparison of the asymptotic development with the exact solution is presented in [41,43]; the agreement is acceptable up to quite high values of the universal constant s.

We shall in the sequel consider the following non-dimensional system parameters: A relation between the phase and group velocities is obtained after a lengthy calculation starting from the dynamical equilibrium equation and the definition of the group velocity.

It can be seen from Fig. 7.6 that for a weak nonlinearity, the phase velocity v p exceeds its linearized value v p0 and describes a supersonic mode; for a higher nonlinearity, the wave changes from a supersonic to subsonic mode (approximately for s=0.7).

Sensitivity analysis for wave propagation in 2D nonlinear periodic hexagonal networks

The constitutive law will be formulated from the strain energy density evaluated for the hexagonal lattice in a 2D Cartesian plane (the same continualization method as in 1D is used); thereby, a nonlinear wave equation in 2D can be deduced based on the dynamical equilibrium equation; these equations will however not be written here. The methodology for determining the acoustic characteristics (the dispersion relation and the phase and group velocities) for periodical uniform 2D lattices is the same as for the previous 1D situation. One can extend the solution of the 1D wave equation, expression (7.14), to a 2D context adopting the same form of the nondimensional strain N in the y direction, but considering now two different strain amplitudes 00 , AB: We shall next analyze the effect of the angle  which parameterizes the hexagonal lattice (Fig. 7.2) and the slenderness ratio (namely the ratio of the beam thickness to beam length) on the frequency band structure. We firstly present the evolution of the frequency versus  for two different values of the wavenumber and for both the subsonic (for s=0.1) and supersonic modes (for s=0.99), in order to investigate the role of the chosen topology on the dynamical properties at the level of the computed effective continuum. ), resulting in a lower dispersion relation for the shear mode.

    2
We next analyze the impact of the slenderness ratio (the ratio of the beam thickness to its length) on the frequency band structure for subsonic and supersonic modes (for both cases of weak and high nonlinearities). Inspection of Fig. 7.9 shows that the limit (s=0.7) of the transition between the subsonic and supersonic modes does not change with parameters  and . Values smaller (resp. larger) than 0.7 characterize the domain of validity of the supersonic (resp. subsonic) modes; these results are in very good arguments with those obtained in Fig. 7.4. As expected, when passing from the reentrant configuration to the regular hexagonal lattice, a pronounced shift occurs in the dispersion relation; the same evolution is also observed by increasing the slenderness ratio  . results are in very good agreement with those obtained in [17], where the authors demonstrate the dependency of the dispersion curves on the wave amplitude in the nonlinear effective medium.

Anisotropic behavior of the hexagonal lattice based on the nonlinear effective media

In order to highlight the rather anisotropic behavior of the re-entrant version of the hexagonal lattice as well as the non-dispersive propagation in the linear case, we plot in Fig. 7.11 the phase velocity for both classical and reentrant configurations; we first consider linear wave propagation, obtained from the nonlinear wave Eq. (7.10) by neglecting the nonlinear elastic moduli E 2 and E 3 . We next analyze similar results for nonlinear waves; the phase velocity for the nonlinear effective medium is represented in Fig. 7.12 for both supersonic (s=0.1) and subsonic modes (s=0.8). density. Differences in the shape of the phase velocity plots when increasing the wavenumber from k=1 to k=4 indicates the dispersive behavior of the nonlinear effective medium.

Conclusion

We have analyzed in the present chapter the propagation of nonlinear waves in homogenized periodic nonlinear hexagonal networks, considering successively 1D and 2D situations. The main originality of this work is the derivation of the nonlinear wave equation from the homogenization of the initial repetitive network.

Wave analysis is performed on the basis of the expression of the identified effective strain energy density in the nonlinear regime. The obtained nonlinear second order gradient continuum has two propagation modes: an evanescent subsonic mode that disappears after a certain wavenumber and a supersonic mode characterized by an increase of the frequency with the wavenumber. For a weak nonlinearity, a supersonic mode occurs and the nonlinear dispersion curves lie above the linear dispersion curve. For a higher nonlinearity, the wave changes from a supersonic to an evanescent subsonic mode and the dispersion curves drops below the linear case and the frequency vanishes at certain values of the wavenumber. An important decrease in the frequency occurs for both subsonic and supersonic modes when the lattice becomes auxetic (negative Poisson's ratio), and the longitudinal and shear modes become then very close to each other.

The present methodology shall be extended in future works to analyze nonlinear wave propagation phenomena in other types of repetitive networks and microstructured media.

Appendix A: Evaluation of the strain energy density of the hexagonal lattice

The expression of the internal deformation energy is written as

                      2 1 2 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 3 1 1 1 1 1 2 2 2 2 2 1 12 1 2 2 2                                                                           i i j j i i i i i i i i j j j j j j j j i i i i j j j j i i i i i i i i j j j j j j j j ES L EI L L L L W L EI L L                                     (A.1)
Based on (A.1), the strain energy of the structure can be expanded as follows:

1 2 2 2 2 2 1 1 1 1 1 1 1 1 2 1 2 2 1 2 2 2 1 1 1 1 1 1 1 1 2 1 2 2 1 2 1 1 1 1 2 2 1 1 2 2 2 2 2 2 i i i i i i i i j j j j j j j j l l i i i i i i i i j j j j j j j j a a K K a a a a W 3 1 2 2 1 1 2 3 1 3 3 1 3 1 1 1 1 1 1 2 1 1 3 1 3 3 1 3 1 1 1 2 1 1 2 2 2 i i i i j j j j i i i i l f j j j j i i i i j j j j a K K a a a 2 3 1 2 2 1 1 1 1 2 1 2 2 1 2 3 1 3 3 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 2 i i i i i i i i f j j j j f j j j j i i i i i i i i c j j j j j j j j K K a a a K 2 3 2 2 1 1 1 1 2 1 2 2 1 2 2 1 2 2 1 2 2 1 1 1 1 3 1 3 3 1 3 3 1 3 3 1 3 1 2 2 1 2 2 i i i i i i i i c j j j j j j j j i i i i i i i i c j j j j j j j j a K a K (A.2)
Here and in the sequel, indices 1, 2 and 3 refer respectively to the first, second and third beam.

One denotes separately the contributions to the strain energy associated to each beam: the energy in extension writes

1 1 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 2 2 2 1 2 2 i i i i j j j j i i i i l f j j j j i i i i j j j j i i i c j j j a K K a a W a a K 2 1 1 1 1 1 1 1 1 1 2 i i i i i j j j j j (A.3)
The contribution to the energy in flexion writes

1 1 1 1 2 1 1 2 2 22 1 2 i i i i j j j j i i i i l f j j j j i i i i j j j j i cj a K K a a W a K 2 1 1 1 2 2 1 2 2 1 2 2 1 2 1 22 i i i i i i i j j j j j j j a (A.4)
The contribution to the energy in shear writes

3 3 3 2 2 1 1 2 3 1 3 3 1 3 1 1 3 1 3 3 1 3 2 1 1 3 1 3 3 1 3 3 1 31 1 1 1 1 2 1 1 2 2 2 2 1 2 i i i i j j j j i i i i l f j j j j i i i i j j j j i cj a K K a a W a K 2 1 1 1 3 3 1 3 3 1 3 3 1 3 1 22 i i i i i i i j j j j j j j a (A.5)
The three beams have here identical rigidities.

Using different finite difference schemes and a Bernoulli kinematic assumption (for more details see [30]), the previous expressions of W 1 , W 2 and W 3 simplify to As a summary of previous computations, the total energy of the hexagonal unit cell writes accordingly as

1 2 3 W W W W (A.9)
This entails the energy density by dividing W, in (A.9) by the area of the elementary unit cell, that is coordinates (for more details, see [30]). In the 1D situation, one obtains after development the expression of the energy density given in equation (7.9).

Appendix B: Elastic moduli

              3 1 2 4sin cos 4cos 3sin 3 4cos cos E E             (B.1)                   3
E E                           (B.2)                   4
E E L                          (B.3)
Appendix C: instability analysis based on Legendre-Hadamard ellipticity condition The condition of strong ellipticity will be used for the wave propagation problem; this condition implies that the modes of propagation of a dynamical equilibrium equation exist so that zero propagation speeds are ruled out. The Hadamard inequality for the strong ellipticity condition is equivalent to the following two independent inequalities, which result from similar inequalities written in a general context in [44] and specialized here to the 1D case:

    22 12
,, 0, 0

xx WW FF F F F F       (C.1)
with     12 , W F W F the strain energy for the first and second order continuum respectively, and 1 F   the transformation gradient in the present 1D context.

As explained in [44], the loss of the strong ellipticity for the first order gradient continuum (ignoring microstructural effects) leads to the existence of a critical value c F , such that

    2 1 0 c W F FF   
. This condition characterizes the existence of discontinuous strain gradients in the equilibrium solutions of elastic solids whose strain energy density is the macroscopic energy density   1 WF, so discarding microstructural effects described by the second order gradient energy density   2 WF. For certain forms of   2 WF, the strong ellipticity condition is restored fort the second order gradient contribution, so that it holds the additional inequality

  2 2 ,, 0 c xx W FF FF      (C.2)
For the specific expression of the strain energy in 1D previously obtained in equation (7.5), incorporating the linear, nonlinear and second gradient parts, and for a high nonlinearity (s close to 1), the nonlinear parts of the constitutive law become dominant, so that the two previous derivatives become:

      2 13 1 2 2 2 2 . ss WW F E E E F FF      (C.3)   2 2 3 ,, 0 s xx W E FF    (C.4)
The second derivative in (C.3) is negative for

      2 1 3 1 , wave 2 1 , 1 , / 0 2 s s c c x c c solitary WW E F k F u f k F E F F               

Extensible energy of plain wave fabric

Although the methodology exposed in the present work is quite general and can be applied to any 3D repetitive architectured material, we shall consider being specific plain weave textile pattern, pictured in Fig. 8.1.

The general methodology for determining the linear and nonlinear energy density of 3D structures is summarized in algorithmic format in the next subsection. . Where E, A, L and I are respectively the Young modulus, the cross section area, the length and the quadratic moment of the considered beam.

        2 2 2 2 2 1 . . . . 2 2 2                                        U e U e U
                                    U e U
For simplicity reasons, we consider that all beams are initially rectilinear, so that no initial deformations or initial curvatures are present, implying that the corresponding kinematics variables ,  cc  are nil in the expression of the normal force (i.e. neglect the coupling between tension and bending).

The expression of the internal deformation energy (see Appendix A) decomposes into extension, flexion and shear contributions expressing as follows;:

        1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 Linear case 1 1 2 2 i i i i i i i i j j j j j j j j k k k k k k k k extension flexion shear extension W W W W EA W u u u u v v w w L L                                                      1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 Linear case 2 Linear case 11 2 11 2 i i i i i i i i i i i i j j j j j j k k k k k k j j j j j j k k k k k k y z z z z flexion z y y y y EI uu LL W EI uu LL                                                                                          1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 3 Linear case 3 Linear case 12 1 1 2 2 2 12 1 1 2 2 2 i i i i i i i i i i i i j j j j j j j j k k k k k k k k j j j j j k k k k k z z z z z shear y yy EI L v v u u L W EI L w w u L                                                         1 1 1 1 2 i i i i j j j k k k yy u                                            
We presently keep a discrete form of the energy; the transition to the continuous energy density is achieved following the steps explained in the Appendix B.

b) Development of previous expressions using different finite difference schemes and a Bernoulli kinematic assumption.

1 1 1 1 1 1 2 2 2 2 2 2 ; 2 2 i i i i j j j j k k k k j j j j k k k k i i i i u a u u a u u u a u u a s s s s 1 1 1 1 1 1 2 2 2 2 2 2 ; 2 2 i i i i j j j j k k k k j j j j k k k k i i i i v a v v a v v v a v v a s s s s 1 1 1 1 1 1 2 2 2 2 2 2 ; 2 2 i i i i j j j j k k k k j j j j k k k k i i i i w a w w a w w w a w w a s s s s 1 1 1 1 1 1 1 1 1 1 1 1
;

i i i i i i j j j j j j k k k k k k j j k k i i z z z

W S W W W

The constitutive law can then easily be obtained; it includes both classical elasticity and second order elasticity, reflected by the existence of the stress and hyperstress tensors computed from the homogenized strain energy density: the first order Piola-Kirchhoff stress is given in component form as the sum of a linear term and a nonlinear contribution:

  linear nonlinear L NL ij ij ij i j WW u x             
The second order hyperstress tensor follows similarly as the partial derivative of the second gradient part of the energy density versus the strain gradient tensor:

  2 2 SG ijk ij k W S u x         where   , , 1, 2,3 i j k 
, and it only includes a linear contribution. This entails overall that the nonlinearity in the effective constitutive law is only included in the first order response, encapsulated in the first order Cauchy stress-strain relation.

Wave propagation in linear and nonlinear 3D homogenized media

The equations of motion of a second gradient medium along the j x directions are written in index form based on the previous identification of the Cauchy stress and hyperstress tensors as follows: where  is the density of the beam.

L NL 2 ij ij ijk j j j k * i linear Nonlinear S , i,j, u S k G 1, 2,3 x x x x                             ( 8
We obtain from Eq. (8.1) three differential equations that describe the propagation of longitudinal waves polarized in the direction of the incident wave, and of two shear waves polarized in a perpendicular direction, namely vertical shear waves in the same plane and horizontal shear waves propagating in a plane perpendicular to the direction of propagation (Fig. 8.2). The three propagations modes are coupled via the kinematics parameters of the rigidity matrices.

Wave propagation in the linear effective medium

We start our dynamical analysis by the linear case, so that one neglects the nonlinear contribution of Cauchy elasticity in Eq. (8.1), viz 0

NL ij   .
In order to compute the dispersion relation, phase and group velocities from the dynamical equation of motion, one adopts for the solution of the dynamical equilibrium Eq. (8.1) the form of a generalized plane wave. For a harmonic wave propagating along an axis in an infinite planar medium, the generalized displacement field with components ( , , ) q U V W  at a point r is assumed in the following form   , ( , , )

wt i q Q e Q U V W   kr (8.3)
where Q is the vector of the wave amplitude, k the wave vector and w the frequency.

Substituting the generalized plane wave Eq. ( 8.3) in the dynamical equilibrium equations delivers the wave equation motion on the following algebraic equation: Without the effect of second order gradient continuum, the obtained effective medium is nondispersive; due to the proportionality of the frequency to the wavenumber k in this situation, the phase and group velocities become equal.

, 1, [ ( ) 0 3 ] 2, U DV k k k W          (8.

Wave motion in the nonlinear 3D homogenized medium

The equations of motion derived from Eq. (8.1) are the equations of a nonlinear second order gradient medium(such that 0 NL ij   ), representing a 3D Boussinesq-type equation, which differs from the classical elasticity theory by the additional terms corresponding to the second order gradient continuum and the nonlinearity contribution of the beam energy at large strains. After insertion of the effective constitutive law for the stress and hyperstress tensors, the equation of motion (1) for the 1D situation representative of a microstructured beam deforming along its axis in tension and compression only takes the form

2 2 4 2 * 1 2 3 2 2 4 2 u u u u u E E E x x x x t               (8.5)
with x the axial coordinate, and the homogenized moduli 1 2 3 ,, E E E therein given versus the microstructural parameters of the unit cell by An explicit analytical stationary solution of the governing nonlinear wave Eq.(8.5) is obtained,

      22 2 1 2 3 2 8 3 2 1 ; ; 3 3 12 a a l l a aa K k a K k L k a E E E K k a L L L         with ,
using the change of variable 1 2 3     z k x k y k w t
 , where  is the angular frequency, k1, k2 and k3 are the wavenumber in the x, y, z directions respectively (explain in detail in the chapter 7).

The general solution of the dynamical equilibrium equation (8.5) The methodology exposed so far for solving the nonlinear wave Eq. (8.5) is quite general; it will be applied to be more specific to the plain wave textile pattern in the sequel.

Dispersion relations and group and phase velocities evolutions for plain weave

The dispersion relations, phase and group velocities are evaluated in the 3 directions of the Cartesian basis  

,, e e e in order to show the influence of the universal constant s (degree of nonlinearity) on the three propagation modes (L: longitudinal mode, SV: vertical shear mode, SH: horizontal shear mode).

We plot in Fig. 8.5 the dispersion relation versus the wavenumber k for different values of the universal constant and for the three propagation modes (L, SV, SH). Two sets of modes can be observed (Fig. 8.5): an evanescent subsonic modes for higher nonlinearity for which the wave velocity is lower than that in the linear case (corresponding to a nil value of the s parameter), and supersonic modes for weak nonlinearity, with waves propagating at higher velocity compared to the linear case. Inspection of the deformation energy in the present large strains regime (previous equation in step d) shows that the energy for wave propagation in the longitudinal direction is greater than for the two other propagation modes, thus resulting in a higher dispersion relation and phase velocity for the longitudinal mode (for both cases subsonic and supersonic). As for the linear case, one can observe from the dispersion relation in Fig. 8.5 a higher value of the longitudinal and vertical shear modes (and a lower value for horizontal shear).

The same results are obtained for wave propagation along 2 e direction due to the geometrical symmetry of the studied textile. Similar shapes as in Fig. 8.5 are obtained, in case of wave propagation in the vertical 3 e direction, but the longitudinal wave shows the lowest dispersion relation. It can be seen from Fig. 8.6 that for weak nonlinearity, the phase velocity increases, corresponding to a supersonic mode; for each propagation mode (L, SV, SH), the phase velocity exceeds the velocity obtained for the linear case, while for a higher nonlinearity, the wave passes from the supersonic to the subsonic mode and the phase velocity decreases, vanishing for specific values of the wavenumber k. For the same reason as discussed before, waves propagate with the highest phase velocity for the longitudinal mode.Fig. 8.6 also shows that SH and SV modes vanish early compared to the longitudinal mode, itself vanishing for a quite high value of the wavenumber k.

We shall note that the same results are obtained in the 1 e direction, while in 3 e direction, a permutation between the mode occurs, and the highest phase velocity occurs for the SH mode.

We evaluate in the sequel the group velocity g C versus the wavenumber k, in order to relate the dispersion relations obtained in Fig. 8.5 to the flow of energy within the effective medium. One can observe from Fig. 8.7 that the group velocity tends to zero when the wave angular frequency of all evanescent subsonic modes approaches its maximum value. The vanishing of mode occurs when g C   , which means that the total energy of the medium is completely dissipated and no flow of energy occurs any more. In the supersonic mode, as the wavenumber k increases, the accumulated energy increases in the medium (the group velocity increases), which entails an increase in the dispersion relation. Higher values of the group velocity in the longitudinal mode lead to higher energy flow (this constitutes a definition of the group velocity) in the medium and then to higher dispersion relations. It is further instructive to compare different effective continua in terms of the dispersion relations.

In comparison with the present second order gradient continuum, the classical continuum is not dispersive since the frequency is directly proportional to the wavenumber k. For the nonlinear micropolar or Cauchy medium, only the supersonic mode occurs, and the dispersion relation lies above the linear case. When the effect of nonlinearity is neglected (restricting to small strain energy), the solution of Eq.(8.1) consists of planar harmonic waves, and three non-dispersive modes propagate within the material with a frequency increasing with the wavenumber (Fig

We show in Fig. 8.9 the evolution of the phase velocity for a wave propagation in both  

12 , ee and  

8.3Anisotropic and dispersive behavior of the second gradient medium for 3D plain weave

The analysis of anisotropy and dispersive effects for the 3D textile is first done in the linear case regardless of the material and geometrical nonlinear aspects. We plot in Fig. 8.13 the phase velocity for two values of the wavenumber, k=1(solid line) and k=3 (dashed line). The irregular form of the shape of the phase velocity in the three modes shown in Fig. 8.13 highlights the anisotropic behavior of the 3D textile. The increase of wavenumber k entails the increase of the phase velocity shown in Fig. 8.13, leading to the dispersive behavior exhibited by the second order gradient continuum. Regardless the value of the angle  , the same responses can also be observed.

In order to highlight the non-isotropic behavior of the nonlinear 3D textile configurations, their dispersive characteristics as well as the influence of the degree of nonlinearity s on these features, we plot in Fig. 8.14, Fig. 8.15 and Fig. 8.16 the phase velocity (propagation in the x-y plane) for low, medium and high values of the wavenumber: k=1, k=1.5 and k=3. As expected, as the wavenumber increases (k=3), another propagation mode disappears, namely the SV modes (Fig. 8.16). The dispersive behavior of the medium leads to changes in the shape and values of the phase velocity plots. One can also note that all 3 propagation modes disappear as the wave number further increases for a high nonlinearity (s=0.999, s =0.95), whereas they do exist for a weak nonlinearity and for any wavenumber.

One shall note that the anisotropic behavior for the 3D textile doesn't change with the wavenumber k, but this anisotropic is less significant in the situation of weak nonlinearity.

288

The relative shapes of the iso-frequency contours for the textile (in plane propagation   12 , ee ) for low and high nonlinearity situations and for low and high wave amplitudes are compared in Fig. 8.17.

Fig. 8.17 Dispersion iso-frequency contour for the textile structure for supersonic and subsonic modes, a) low amplitude wave, b) high amplitude wave.

We can conclude that moving from a weak to high nonlinearity and from low to high wave amplitudes, the degree of anti-symmetry of the textile increases, as shown by the irregular shape of the iso-frequency plots in Fig. 8.17.

It has to be noted that the validity of previous analyses is limited to low frequency elastic wave propagation (corresponding to long wavelengths),for which one is entitled to use an equivalent homogeneous model to describe the actual heterogeneous medium and study wave propagation within the material. A nonlinear second gradient homogenized continuum has been constructed by homogenization of the initial repetitive network material, which provides suitable results when the size of the heterogeneities is lower enough in comparison to the size of the representative volume element (RVE).

Conclusion

The general objective of this chapter is the analysis of wave propagation phenomena within architectured media, relying on an effective substitution continuum obtained by homogenization.

The proposed methodology is quite general and applicable to any 3D repetitive network of beamlike structural elements, considering beams within the network undergoing large displacements and large rotations. Based on the writing of the equations of motion, we have analyzed nonlinear wave propagation in the obtained homogenized nonlinear second order gradient continuum. The resulting wave equations are of Boussinesq type, with solution identified as elliptic functions.

The influence of the degree of nonlinearity on the dispersion relations has been analyzed, highlighting subsonic and supersonic modes propagating respectively with a velocity lower (resp.

higher) than the velocity of linear non-dispersive waves. Subsonic and supersonic modes correspond respectively to regimes of high and low nonlinearity characterized by the so-called universal constant. The wave propagation analysis has been done in a systematic manner, leading to a rich sets of results providing a comparison between longitudinal and shear waves.

Important results of these analyses include a more marked anisotropy when the degree of nonlinearity increases and the disappearance of the horizontal and vertical shear modes when the wavenumber increases.

Although the considered second order substitution medium proves convenient to simplify the initial microstructure, it certainly brings simplifications in wave propagation analysis, so that the present study shall be complemented in future works by fully resolved finite element computations of wave propagation over the whole lattice within the identified unit cell.

Wave analysis in repetitive network materials prone to geometrical nonlinearities shall be given special interest in future developments relying on the identification of hyperelastic models based on adequate homogenization schemes.

extension flexion shear i W W W W      i i i (A.1) such as 0  flexion W           1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 1 1 1 1 1 3 6 11 2 2 2 1 1 2 2 1 1 2 2                                                   i i i i i i i i i i j j j j j j j j k k k k k k k k j j k k extension l extension extension extension extension la W k u u u u v v w w W W W L W k u u u a               1 1 2 2                                                                  i i i i i i i i i i i i i i j j j j j j k k k k k k j j j j j j j j k k k k k k k k extension extension la u v v w w W W k u u u u v v w w a           1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 10
                                                   i i i i i i i i i i i i j j j j j j j j k k k k k k k k j j j j k k k k extension extension la extension extension la W W k u u u u v v w w W a W k u u u u a       1 1 1 1 1 1 1 1 1 2 2 2 2 7 7 7 7 8                          i i i i j j j j k k k k extension v v w w W              1 11 2 2 2 11 2 2 2 i i i i i i i i i i i i i i j j j j j j j j k k k k k k k k j j j j j j k k k k k k cL z z z z shear cL y y y L k v v u u W L k w w u u                                                                          11 1 1 1 1 11 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 6 2 11 1 2 2 2 2 2 2 2 2 2 2 22 11 2 2 2 1 2 ii i i i i i i i i i jj kk j j j j j j j j k k k k k k k k j k shear shear shear y ca z z z z shear ca W W W a k v v u u W k w w                                                                             1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 12 2 2 2 2 2 2 2 1 22 i i i i i i i j j j j j j j k k k k k k k shear y y y y W a uu                                                              1
                                                         i i i i i i i i i i i i i i i j j j j j j j j k k k k k k k k j j j j j j j k k k k k k k k ca z z z z shear ca y y y y a k v v u u W a k w w u u           1 10 2                             i j shear W              1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
                                                   i i i i i i i i i i i i i i i j j j j j j j j k k k k k k k k j j j j j j j j k k k k k k k k ca z z z z shear ca y y y y a k v v u u W a k w w u u                      1 11
                                                                i i i i i i i i i i i i i j j j j j j j j k k k k k k k k j j j j k k k k shear ca z z z z shear ca y y W a k v v u u W a k w w         1 1 1 1 11 11 8 2 7 7 7 7 1 2                                  i i i i j j j j k k k k shear yy W uu Using a finite difference scheme, one can next write 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 1 1 2 2 2 2 2 2 1 2 2 k k k k k k k k extension l j j j j j j j j k k extension la j j u u u u v v w w L L L L W k L L L L s s L s s s s s s u u a W k a s s                                                                                   k k k k k k j j j j j j k k k k extension la j j j j u u v v w w a a a a a a a s s s s s s u u u u a a W k a a s s a s s                                                                                               2 2 2 2 2 2 2 2 4 4 4 4 2 2 2 2 k k k k j j j j v v w w a a a a s s s s                                                2
                                                                                 2 2 2 2 2 2 2 2 2 2 7 7 7 7 7 7 2 2 2 1 2 2 2 2 k k k k k k j j j j j j u u v v w w a a a a a a a s s s s s s                                                                      2 2 2 2 2 1 1 1 1 2 2 1 2 2 2 2 1 1 1 1 1 2 2 1 2 2 2 1 2 2 2 2                                                                                                i i i i j j j j k k k k i i i i i j j j j j k k k k k cL shear cL v u u v L L k L s s s s W w w u u w L L k L L s s s s s 2                    u u u v v cos v sin w w w sin v cos w             (B.5)
This entails the first and second derivatives

1 1 1 2 2 2 1 3 3 3 [1] [2] [3] [1] [2] [3] [1] [2] [3]                                                     i i i N i i i i i k k i i i u u u Y Y Y x y z u u u u u Y Y Y s s x y z u u u Y Y Y x y z (B.6) and [1] [2] [3] 2 [1] [2] 2 [1] [3] 2 [2] [3] [1] [2] [3] 2 [1] [2]                                                 i i i i i i i i i N ii k k u u u u Y Y Y Y Y x y z x y uu Y Y Y Y x z y z u u u u Y Y Y Y Y x y z uu ss 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 2 2 2 2 2 3 3 3 3 2 [1] [3] 2 [2] [3] [1] [2] [3] 2 [1] [2] 2 [1] [3] 2 [2] [3]                                                                                 i i i i i i i i xy u u Y Y Y Y x z y z u u u u Y Y Y Y Y x y z x y u u Y Y Y Y x z y z             (B.7) with   1..12  i the index of the beam, 1, 2,3 N indicates the space dimension; in the present 3D context so N3  , 12 1 1 2 2 12 [1] [1] [2] , [2] [3] [3]                      YY Y Y Y Y YY and 3 33 3 [1] [2] [3]       Y YY Y
represent the periodicities vectors. For the investigated example, 12 .

Overview

The analysis of wave propagation in hyperelastic media depends initially on the type of constitutive law. When considering microstructured solids, the effective constitutive law written in the large strains range reflects the impact of the microstructure, and can be obtained thanks to suitable homogenization schemes.

In recent years, different materials have been analyzed in the context of anisotropic finite-strain elasticity. These include composites, foam-like structures, 2D and 3D textile preforms and synthetic solids. Cellular solids, by contrast to compact materials, are two or three dimensional bodies divided into cells, the walls of which are made of a solid material capable of undertaking large elastic deformations without plastic failure or fracture. There are numerous examples of such network structures, including repetitive large scale deployable structures like antenna, 3D textiles, cellular materials and especially auxetic structures having excellent damping and impact absorption capabilities.

We shall employ the discrete asymptotic homogenization technique [1][2][3][4][5] which is perfectly suited to the discrete architecture of different types of networks, and which is extended to the nonlinear range in the present work. Due to the very small bending rigidity of the beams building such networks, it is reasonable to assume that the nonlinear response is essentially due to the change of the network configuration (the orientation and length of the beams change with ongoing deformation), thus we shall mostly account for geometrical nonlinearities at the microlevel. The geometrical nonlinear behavior of cellular structures and network materials was extensively studied by [6,7], considering the example of foams, using simplified pin jointed model for which the bending contribution of the skeleton struts was neglected. Wang and Cuitino [8] proposed another approach where axial, bending and twisting deformations at local level were considered. One study based on an homogenization technique was given in [9]. Linear effective models to analyze structures on the basis of a beam model were presented in [10,11], in which stretching and simultaneous bending occur. More recently, Janus-Michalska [12] extended this linear model to construct the stress-strain relation and strain energy function for the hyperelastic cellular material with arbitrary symmetry. An alternative approach was proposed by Vigliotti et al. [13] using a computational homogenization to derive a nonlinear constitutive model for lattice materials. The homogenization is the periodic network is done towards a Cauchy continuum at the mesoscopic level, whereby the nodal rotations are condensed and finally expressed versus the deformation applied over the unit cell, using the equilibrium equations. Details related to the main elements of the asymptotic homogenisation technique considering the micropolar framework can be found in [1,5]. T respectively, the first and second order terms according the power of  homogenized continuous functions.

                                                                         

Large strains effective response of network materials

We next expose the extension of the previous DH method to the nonlinear regime, in order to account for the large changes of configuration of the networks. The basic idea is to update the lattice geometry at each new increment of the external load applied to the unit cell boundary, from which a linearized elastic computation will be done over the same load increment. Since details of the method have been presented in [14,15], we only recall the main thrust of the method in the present contribution. The main steps of the DH method leading to the nonlinear response of the homogenized continuum are written in algorithmic format in Box 1. Note that although the main source of nonlinearities at microscopic level is the modification of the network geometry, the obtained constitutive law is a nonlinear relation between stress and strain.

A dedicated code has been constructed from the proposed algorithm to solve for the nodal kinematical unknowns (displacements) of each beam within the repetitive unit cell. The code is written in symbolic language uses an input file including the initial reference unit cell topology and mechanical properties, and delivers as an output the homogenized mechanical response in both the linear and nonlinear regimes (the classical moduli).

For each kinematic increments

  k Gn E ;
For each iteration k inside the increment loop;

1. Compute the effective mechanical properties in the linear regime based on the linear discrete homogenization framework [1,5].

2. Define the incrementally imposed strain applied over the RUC. 8. Repeat steps 1-7 up to the maximum applied strain and curvature over the unit cell.

Box.1. Algorithm for the nonlinear discrete homogenization of repetitive lattices.

Identification of a hyperelastic strain energy density for the hexagonal lattice, the reentrant lattice and plain weave textile

Since the DH method is predictive, it can be conceived as a virtual testing method (instead of doing real measurements, which can be costly) to provide a database of uniaxial loading response to identify a strain energy density for an assumed hyperelastic effective homogeneous material.

We shall calibrate a strain energy function of two preselected hyperelastic models for the three investigated lattices (Fig. 9.2). In the present 1D context of microstructured beams operating under pure tensile loadings, the sole degree of freedom is the scalar displacement along the beam, variable () ux, with a spatial gradient denoted by the scalar quantity ,x u (the comma denotes the partial derivative).

The following two forms of the hyperelastic function Note that there is no unique choice of the mesoscopic (homogenized) constitutive law; observe that the adopted choice of the strain energy density there above means that we restrict the mesoscopic stress (the second Piola-Kirchhoff stress) to be at most a quadratic function of the Lagrangian strain (its conjugated strain).

The constitutive law can be derived from the form taken by the strain energy density of the hyperelastic model; the first Piola-Kirchhoff stress tensor is computed as the partial derivative The function   T A, B,C there above is built as the quadratic measure of the error between the DH stress component DH T and its analytical counterpart, obtained from the hyperelastic potential.

The material constants of the model are identified from a least square method (9 sampling points are used), relying on uniaxial tension as the kinematic loading imposed over the lattice unit cell.

the function Lsqcurvefit in the Optimization Toolbox of MATLAB has been used [14,15] to identify the coefficients (A, B and C) of the strain energy density for these three lattices (Fig. 9.2) based on the incremental scheme developed in section 9.2; they are listed in Table 9.1. We shall in the next section rely on the two different forms of the strain energy density to analyze nonlinear wave propagation.

Analysis of nonlinear wave propagation in the homogenized hyperelastic continua

Different types of nonlinear wave propagation equations are considered in this work: harmonic plane waves based on the perturbation method [29], solitary waves for the Boussinesq type equation [32], or shock waves for Burger's equation.

We shall consider the following non-dimensional system parameters:  the Young modulus, density and length of the beam structures respectively.

We first evaluate the constitutive law for the two forms of the previously identified strain energy functions.

Wave propagation analysis for the form 1 of the hyperelastic effective medium energy

Considering first form1 of the nonlinear strain energy without the second gradient terms, we can write the dynamical equilibrium equation as: The first order equation at order   0 O  describes linear wave propagation in the effective (linear) medium.

We take the plane harmonic wave as a solution of the   The dispersion relation for the three investigated lattices involves an amplitude-dependent frequency in the context of a nonlinear effective medium, based on Eq.9.15.  ; this behavior can be observed for the three investigated lattices.

Wave propagation analysis for form2 of the hyperelastic energy

Recall that the strain energy density is selected as a cubic function of the linearized strain     It appears from Fig. 9.6 that for all values of parameter s between 0 and 1, the supersonic mode always occurs; these results are in very good arguments with those obtained in Fig. 9.5.

Effect of the second gradient terms on wave propagation for the hyperelastic medium

We shall in the sequel consider the linear second gradient effective rigidity into the strain energy in order to study their impact on the propagation of nonlinear waves.

The form1 and form2 of the nonlinear constitutive laws accounting for the second gradient terms can be written respectively as follows: where A, B,C are the same coefficients written in Table 9.1, and  accounts for the second gradient terms calculated in details in [33] (see Table 9.2).

Wave propagation analysis for form1 with additional second gradient terms

The dynamical equilibrium equation for the second gradient effective medium can be written as follows: Inserting the constitutive law (Eq.9.23) into the dynamical equilibrium, Eq.25, further leads to  

u u u OB x u u u u u u O B C x x x x x                                         (9.27)
The first order equation at order   0 O  describes linear wave propagation in the effective medium with second gradient terms which the plane harmonic wave is solution.

The second order equation at order  

1 O  delivers a relation between the corrected frequency 1 

and the linear frequency 0  (the same as expression (9.15)). For a weak nonlinearity, a supersonic mode occurs and the dispersion curves lie above the linear dispersion curve (when v = v 0 ). For a higher nonlinearity, the waves changes from a supersonic to an evanescent subsonic mode and the dispersion curves drops below the linear case, and further vanishes for certain values of k. The supersonic and subsonic modes are observed for the three networks, with variation in the dispersion curve, especially the vanishing of the subsonic mode for the hexagonal lattice appears earlier in comparison to the two other lattices.

In the nonlinear case, the introduction of the second order gradient parameters leads to two propagation modes (subsonic and supersonic), whereas only the supersonic modes occur for a pure Cauchy continuum. Fig. 9.9 shows the evolution of the frequency with the degree of nonlinearity s for the hexagon lattice, the re-entrant lattice and the textile plane weave pattern. Fig. 9.9 Evolution of the frequency of the longitudinal mode versus parameter s for the hexagonal network (red), the re-entrant network (green) and plain weave textile (blue).

Fig. 9.9 makes clearly apparent that the limit of the supersonic mode depends on the type of network. The s-limits of the supersonic mode are 0.78, 0.9 and 0.99 for the classical hexagon, the re-entrant hexagon and the textile structure respectively. For values of s smaller than these limits, the supersonic mode occurs; when s exceeds these limits, there is a transition from the supersonic mode to the evanescent subsonic mode. It is further apparent from Fig. 9. 9 that the evanescent subsonic mode occurs for a higher value of the degree of nonlinearity for the textile plane weave.

Dispersion analysis for the first and second order gradient hyperelastic models

In order to highlight the dispersive behavior of the nonlinear effective medium, we plot in

Conclusion

We analyzed in this chapter nonlinear wave propagation occurring within microstructured beams including a repetitive network undergoing configuration changes under pure tensile loadings, leading to geometrical nonlinearities. Three types of repetitive microstructures have been considered in order to exemplify the analysis of (nonlinear) wave propagation: the hexagonal network, its re-entrant version, and plain wave textile.

The effective nonlinear constitutive law have been identified from a micromechanical scheme in terms of the strain energy density expressed as a nonlinear function of the small strain tensor; first and second order grade 1D homogenized continuum have been thereby identified in the nonlinear range, based on a strain driven incremental scheme written over a reference unit cell taking into account the variation of the lattice geometry. The coefficients of the selected constitutive models have been identified for three specific network materials based on the proposed homogenization scheme.

The dynamical equations of motion have been next written for two different forms of the constitutive law -successively discarding and including second order terms -originating from the selected strain energy density, leading to specific dispersion relations. The inviscid Burgers equation is obtained as a specific wave propagation equation for the first order effective continuum, whereas Boussinesq equation emerges when adopting a nonlinear strain gradient effective continuum. The presence of second gradient order terms in the nonlinear equation of motion leads to the presence of two different modes: an evanescent subsonic mode for a high nonlinearity, vanishing beyond a given value of the wavenumber k, and a supersonic mode for a weak nonlinearity.

The generalization of such nonlinear wave propagation analyses to 2D situations will be performed in future work, based on suitable hyperelastic constitutive laws obtained from the homogenization of the existing microstructure.

Chapter 10: Conclusion and perspectives 10.1 Summary

The analysis of wave propagation within 2D and 3D periodic network materials has been done, taking into account microstructural scale effects and a viscoelastic behavior of the network, relying on a homogenization method specific to discrete structures to construct substitution media of the initial discrete network. Two different strategies for the enrichment of the effective continuum have been considered, namely micropolar and second order gradient continua. The dynamic analysis of random fibrous networks is also presented for 2D case. The propagation of linear and nonlinear waves within the network represented by the substitution continuum has been analyzed in two successive parts of the document.

Our contribution focuses on two important issues: the impact of the enrichment of the effective continuum on the dispersion relations has been analyzed for 2D and 3D network materials, considering micropolar and second gradient effective continua. This methodology has been especially applied to metamaterials showing an auxetic behavior, which show an enhanced partial band gap at low frequency. The second important issue is the analysis of nonlinear wave propagation within repetitive networks undergoing configuration changes under large applied kinematic loadings, and reflecting geometrical nonlinearities. In the second part of the work, an incremental scheme for the update of the frequency and phase velocity of the homogenized medium has been developed, considering successively 1D, 2D and 3D repetitive network materials. We have furthermore constructed hyperelastic media by virtual tests relying on the homogenization method at the basis of this thesis, and obtained different types of wave propagation equations (Burgers and Boussinesq equation) and solutions (shock waves, solitary waves, and harmonic planar waves), the acoustic properties of which has been studied.

Research contributions

The research presented in this dissertation provides the following original contributions:

 A general methodology to analyze dispersion effects in network materials showing microstructural effects, based on the construction of second order gradient effective medium;  Wave propagation analysis in dissipative networks also accounting for microstructural effects;

to composite materials made of such repetitive networks as the reinforcement embedded into a softer elastic or viscoelastic matrix.

Although the research presented here offers in our opinion significant advancements in linear and nonlinear dispersion analyses, there are a number of exciting opportunities available for further exploration. The tools and investigations presented have brought dispersion analysis in microstructured solids showing internal lengths and / or nonlinear effects to the point of experiment design and implementation to tests the obtained numerical predictions. Several fundamental theoretical investigations and possible experimental sets up are discussed next as potential options for future work.

Experimental investigation of wave dispersion and attenuation

It remains to validate the predictions of the simulations by suitable experiments; this shall be done with regards to two main aspects:

-Measurements of band diagrams;

-Measurements of the transmission loss factor.

There are several existing experimental set up to measure the transmission loss factor due to the tortuosity of the tested microstructures and / or the internal damping due to the viscoelastic material behavior. We envisage designing an experimental set up allowing to measure the acoustic properties of structures submitted to stationary waves (Kundt tube) tested in situ under uniaxial or even biaxial loading conditions. In parallel to this, we shall attempt to extend Biot's model to higher gradient continuum models for porous microstructures to have at hand a suitable modeling framework to interpret the obtained measurements.

Regarding nonlinear aspects, the design of sensitive nonlinear network materials should be at the forefront of future experimental analyses.

Wave propagation in generalized continua constructed by homogenization

In many situations, second gradient media are not sufficient to account for microstructural effects, especially in situations for which a microstructure is present at different (more than one) scale levels; this is true for both cases of network materials and composites. We shall accordingly extend the class of constructed generalized continua to micro stretch, micro dilatation and micromorphic continua, constructed by suitable homogenization method, either from analytical or numerical inspiration. Micromorphic effective continua are further prone to lead to full band gaps at low frequencies; it will accordingly be of interest to develop micro stretch and micromorphic continua of microstructured solids in order to compute their band diagram.

Elaboration of a rationale for the design of metamaterials

Beyond this aspect, elaboration of design rules of metamaterials giving rise to significant band gaps requires the elaboration of a suitable methodology combining homogenization schemes and topology optimization. Especially, the construction of proper higher order schemes (second order gradient effective models) that does not depend on the microstructure remains a challenging task.

In many cases, second gradient effects are quite limited and conceiving microstructures for which boundary layer effects (associated to second order gradients) are significant in comparison to the unit cell size remains challenging. We have merely focused on a specific class of metamaterials, namely auxetic, but it is clear that the family of conceivable metamaterials is much wider. Large deformations give rise to instabilities inducing band gaps, as preliminary results in this work have

shown. The conception of microstructures with a state of anisotropy varying with the imposed deformation and prone to such instabilities leading to negative tensile or bulk moduli remains a challenging perspective.

Impact of disorder and defects on acoustic properties

In this work, we have considered perfectly periodic structures, in the absence of defects. It is nevertheless clear that the presence of disorder like topological fluctuations (encountered e.g. in biological membranes), or damage and microcraks will alter both the static and dynamic properties in such imperfect networks; this remains a challenging task to be addressed in the future. Some wok has been recently initiated in this direction with the group of Dr. S. Phani (D pt of Mechanical Engineering, UBC, Vancouver) in order to investigate from theoretical and experimental point of view the impact of existing cracks on wave propagation features.

Dynamic homogenization

The prediction of the acoustic properties based on the homogenized continuum is only descriptive of the low frequency dynamic behavior of heterogeneous periodic media, in conditions when the scale separation between the displacements and the unit cell size holds true. In the regime of higher frequencies, this scale separation is lost and the inertia associated with the local motion become quite significant (the effect of inertia can be neglected for low frequency) for that the proposed homogenized method cannot be applied.

The dynamic homogenization taking into accounts micro-inertia effects and space-time interactions for heterogeneous periodic media under dynamic loading and covering the higher frequency regime remains one of the most challenging issues to be addressed in the future.

Title: Computation of the effective static and dynamic properties of network materials accounting for microstructural effects and large deformations'

Abstract: Micropolar and second gradient effective continua are constructed as two different strategies to account for microstructural effects. The influence of additional degrees of freedom or higher order displacement gradients on the dispersion relations is analyzed in both situations of elastic and viscoelastic behaviors of the material. Generalized effective continua lead to dispersive waves, as observed in experiments. In the second part of the thesis, we analyze the influence of large deformations on the propagation of acoustic waves in repetitive network materials. Both theoretical and numerical methods are developed in order to assess the influence of finite strains developing within the networks on the evolution of their band diagrams. An incremental scheme for the update of frequency and phase velocity of the computed homogenized medium is developed based on a perturbation method for 1D, 2D and 3D structures, considering with a special emphasis auxetic networks. This scheme shows an important effect of the applied finite deformation on the frequency and phase velocity of the propagating waves. A perturbation method for nonlinear periodic structures is developed to extend Bloch's theorem to cover both geometrical and material nonlinearities. Hyperelastic first and second order gradient constitutive models of different network materials are identified based on dedicated homogenization methods, from which specific wave equations are formulated -Burgers and Boussinesq equations -the dispersion properties of which are analyzed.

9. Chapter 9 : 10 .

 910 Nonlinear wave propagation analysis in hyperelastic 1D continuum materials constructed by homogenization ................................................................................ 9.1Overview ............................................................................................................................. 9.2 Incremental scheme for the computation of the effective hyperelastic effective models .. 9.2.1 Large strains effective response of network materials ................................................. 9.3 Identification of a hyperelastic strain energy density for the hexagonal lattice, the re-entrant lattice and plain weave textile .................................................................................................. 9.4 Analysis of nonlinear wave propagation in the homogenized hyperelastic continua......... 9.4.1 Wave propagation analysis for the form 1 of the hyperelastic effective medium energy ............................................................................................................................................... 9.4.2Wave propagation analysis for form2 of the hyperelastic energy ................................ 9.5 Effect of the second gradient terms on wave propagation for the hyperelastic medium .... 9.5.1 Wave propagation analysis for form1 with additional second gradient terms ............. 9.5.2 Wave propagation analysis for form2 with second gradient terms .............................. 9.6 Dispersion analysis for the first and second order gradient hyperelastic models ...............9.7 Conclusion .......................................................................................................................... References ................................................................................................................................ Chapter 10: Conclusion and perspectives ..................................................................... 10.1 Summary ........................................................................................................................... 10.2 Research contributions ..................................................................................................... 10.3 Recommendations for future work ................................................................................... 10.3.1 Experimental investigation of wave dispersion and attenuation ................................ 10.3.2Wave propagation in generalized continua constructed by homogenization .............. 10.3.3Elaboration of a rationale for the design of metamaterials ......................................... 10.3.4Impact of disorder and defects on acoustic properties ................................................ 10.3.5Dynamic homogenization ...........................................................................................

Fig. 1 . 1 18 Fig. 1 . 2 Full 20 Fig. 2 . 1 36 Fig. 2 . 2 40 Fig. 2 . 3 42 Fig. 2 . 4 48 Fig. 2 . 5 51 Fig. 2 . 7 52 Fig. 2 . 8 55 Fig. 2 . 59 Fig. 2 .Fig. 3 . 7 Fig. 3 . 8 Fig. 3 . 9 Fig. 3 .Fig. 3 .Fig. 3 .Fig. 4 . 3 Fig. 4 . 5 Fig. 4 .)Fig. 5 .Fig. 6 . 1 Fig. 6 .Fig. 6 .xxEFig. 6 .Fig. 6 .Fig. 6 .Fig. 7 . 5 .Fig. 7 . 6 .Fig. 7 ..Fig. 7 .Fig. 8 . 1 Fig. 8 . 9 Fig. 8 .Fig. 8 .Fig. 8 .Fig. 8 .Fig. 8 .

 11181220213622402342244825512752285525923738393334345456166666757677818988888 Fig. 1.1 Periodic structure composed of a repetitive unit cell ..................................................................... 18 Fig. 1.2 Full Band gap in periodic structures. ............................................................................................. 20 Fig. 2.1 Kinematic and static beam variables. .............................................................................................. 36 Fig. 2.2 Shift factor with respect to the reference lattice. ........................................................................... 40 Fig. 2.3 Representative unit cell of the investigated lattices: (a) hexagonal lattice, (b) re-entrant lattice   0   c) Diamond chiral lattice d) Hexachiral lattice. .............................................................................. 42 Fig. 2.4 Damping ratio and frequency band structure for low and high damping situations respectively for the classicalhexagonal lattice and the reentrant lattice along the edge of the first Brillouin zone for the a) Rotary mode, b) Longitudinal mode, c) Shear mode. ................................................................................. 48 Fig. 2.5 Frequency band structure and damping ratio in damping situations e = 500Mpa.sec for the a) classical hexagonal, b) the reentrant hexagon lattice along the edge of the first Brillouin zone. Red line: Rotary mode -Green line: longitudinal mode -Blue line: shear mode. ...................................................... 49 Fig. 2.6 Frequency band diagram and damping ratio for the four considered lattices, a) Rotary mode, b) Longitudinal mode, c) Shear mode. ............................................................................................................ 51 Fig. 2.7 Damping ratio for the four considered lattices in a damped medium e = 500Mpa.sec for the three modes of propagation. a) Rotary mode, b) Longitudinal mode, c) Shear mode. .............................. 52 Fig. 2.8 Phase velocity for the three modes for 2 values of µe, a) Hexagonal lattice, b) Hexachiral lattice, c) Re-entrant lattice, d) Diamond chiral lattice. .............................................................................................. 54 Fig. . Eige a es of a t pi al ell a e sho i ta ula fo at i sta ts t= τ i ed a d t = i g ee . The three rows correspond to the three points O, A and B in k space, while the three columns correspond to the three dispersion branches for the hexagonal lattice. ................................................... 55 Fig. 2.10 Unit cell for a) the regular hexagonal lattice, b) Hexachiral lattice. ............................................. 59 Fig. 2.11 Dispersion relation in a damped medium, for the a) Hexagonal lattice, b) Hexachiral lattice. Red li e: lo gitudi al ode. G ee li e: shea ode. Co pa iso et ee Blo h's theo e dashed li e a d homogenization theory (solid line) for two modes of propagation ............................................................ 60 Fig. 2.12 a) Diamond chiral lattice with the selected primitive unit cell, b) First Brillouin zone in the reciprocal lattice. The points O, A, and B and the vectors Y1, Y2 (direct lattice) and b1,b2 (reciprocal lattice) are as defined in the Table below. ................................................................................................... 64 Fig. 3.1 Schematic diagram showing (from left to right) the wave propagation modes: the longitudinal mode, the shear vertical mode, the horizontal shear mode. ..................................................................... 78 Fig. 3.2 Elementary unit cell of, (a) 3D hexagonal lattice (left column), (b) plain weave fabric (right column). ....................................................................................................................................................... 79 Fig. 3.3 Retardation times in traction (a) in x-direction for the 3D hexagonal lattice, (b) in z-direction for the 3D hexagonal lattice, (c) in x-direction for the plain weave, (d) in z-direction for the plain weave. .... 80 Fig. 3.4 Dispersion relation and damping ratio evolutions for the 3 lattices with a direction of propagation θ=π/6, =π/6 for the 5 following modes. (a) A1μ arbitrary wave, (b) Aβμ arbitrary wave, (c) longitudinal wave, (d) vertical shear wave and (e) horizontal shear wave. ..................................................................... 83 Fig. γ.5 Dispersion relations evolutions with a direction of propagation θ=π/6, =π/β for the plain weave textile structure. ........................................................................................................................................... 84 Fig. 3.6 Frequency band structure, in a damped medium with e  =100. (a) regular hexagonal lattice at di e tio θ=π/ , γ=π/ , Te tile latti e at di e tio θ=π/ , γ=π/ , e-entrant lattice at direction θ=π/ , γ=π/ a d d Te tile latti e at di e tio θ=π/ , γ=π/ . .................................................................. 85

Fig. 1 . 1

 11 Fig. 1.1 Periodic structure composed of a repetitive unit cell One of the most attractive features in acoustic metamaterials is the generation of wave attenuation phenomena: within some frequency intervals, waves cannot propagate through the periodic but are attenuated quickly and completely within the periodic cells. The corresponding frequency range forms the so-called frequency band gap. Therefore, a lot of research activities have fostered on creating band gap artificially by constructing acoustics metamaterials with negative values of the constitutive parameters, for instance by adding local resonators within the (periodic) structures: these include split ring resonators for negativity, Helmholtz resonators for negative modulus and a membrane resonator for a negative dynamic mass. These local resonators have significant effect on the frequency band gap and allow exploring the desired band gap in the low frequency regime. The review article by Hussein et al. [2] presents a thorough history of acoustic metamaterials. In

Fig. 1 . 2 Full

 12 Fig. 1.2 Full Band gap in periodic structures.

15 ) 4 .

 154 Order the previous equilibrium equation according to the successive powers of  5. Convert the discrete equilibrium equations after homogenization to continuum self-equilibrium equations of forces and moments in the form of Riemann integrals (Appendix C)

Fig. 2 . 2

 22 Fig.2.2 Shift factor with respect to the reference lattice. Accordingly, the constitutive equations for the equivalent 2D anisotropic micropolar viscoelastic continuum write in matrix format as

Fig. 2 . 3

 23 Fig.2.3 Representative unit cell of the investigated lattices: (a) hexagonal lattice, (b) re-entrant lattice   0   c) Diamond chiral lattice d) Hexachiral lattice.

  1 the mass of the set of lattice beams, A cell being the area of the periodic cell; the micro-inertia terms receives the general expression

12 (

 12 vector and is a complex frequency function (of the wave vector) that permits wave attenuation in time. In the limiting case of no damping, it holds i   , and the usual form of plane waves is recovered. Substituting equations (2.34) and (2.33) into the equation of motion (2.25), one obtains an equation of wave motion of the form wave without attenuation in the x-y plane, the propagation constants along the x and y directions are 11 k   and 22 k   . For a plane wave without attenuation in the Cartesian coordinate system in 2-D, the propagation constants along the x and y kk obtained by solving the eigenvalue problem in (2.35) represents plane wave propagation at the frequency .

Fig. 2 . 4

 24 Fig.2.4 Damping ratio and frequency band structure for low and high damping situations respectively for the classical hexagonal lattice and the reentrant lattice along the edge of the first Brillouin zone for the a) Rotary mode, b) Longitudinal mode, c) Shear mode.

Fig. 2 . 5

 25 Fig. 2.5 Frequency band structure and damping ratio in damping situations e = 500Mpa.sec for the a) classical hexagonal, b) the reentrant hexagon lattice along the edge of the first Brillouin zone. Red line: Rotary mode -Green line: longitudinal mode -Blue line: shear mode.

Fig. 2 .

 2 Fig.2.5 shows a partial band gap occurring between the branches; this partial band gap is more significant between the R and L modes than between L and S modes for the regular hexagonal lattice, while the partial band gap for the re-entrant lattice occurring between R and L modes is

Fig. 2 . 6

 26 Fig.2.6 Frequency band diagram and damping ratio for the four considered lattices, a) Rotary mode, b) Longitudinal mode, c) Shear mode.

Fig. 2 . 7

 27 Fig. 2.7 Damping ratio for the four considered lattices in a damped medium e = 500Mpa.sec for the three modes of propagation. a) Rotary mode, b) Longitudinal mode, c) Shear mode. Results in Fig. 2.7 show that the variation of the slenderness ratio  leads to an important change of the damping ratio for the 4 lattices and for the 3 modes of propagation; as the slenderness ratio increases, the damping ratio increases and shifts in the frequency band structure occur. This can be explained based on the homogenized viscosity matrix: a decrease of the slenderness ratio results indeed in an increase of the viscosity matrix.

Fig. 2 . 8

 28 Fig.2.8 Phase velocity for the three modes for 2 values of µe, a) Hexagonal lattice, b) Hexachiral lattice, c) Re-entrant lattice, d) Diamond chiral lattice.

Fig. 2 . 9

 29 Fig.2.9 Eigenwaves of a typical cell are shown in tabular form at instants t=βτ (in red) and t = 0 (in green). The three rows correspond to the three points O, A and B in k space, while the three columns correspond to the three dispersion branches for the hexagonal lattice.

Fig. 2 .

 2 Fig.2.10 Unit cell for a) the regular hexagonal lattice, b) Hexachiral lattice. Hexagonal lattice Hexachiral lattice 1 2 31 41 , k k q e q q e q  

1 N gives a relation between parameters 12 (

 112 , ,) kk delivering the dispersion relation for each damped lattice. We shall note that we can restrict the values of the wavevector to the edges of the irreducible part of the first Brillouin zone.

  2.11 the dispersion relation for longitudinal and shear modes (these 2 modes only exists at low frequency) b)

Fig. 2 .Fig. 2 .

 22 Fig.2.11 Dispersion relation in a damped medium, for the a) Hexagonal lattice, b) Hexachiral lattice. Red lineμ longitudinal mode. Green lineμ shear mode. Comparison between Bloch's theorem (dashed line) and homogenization theory (solid line) for two modes of propagation Fig.2.11 shows that in the region corresponding to the low frequency range and high wave number, the homogenized method is applicable, since the dispersion curves obtained by Bloch theorem and by the homogenization method are very close to each other.

Fig. 2 .

 2 Fig.2.12 a) Diamond chiral lattice with the selected primitive unit cell, b) First Brillouin zone in the reciprocal lattice. The points O, A, and B and the vectors Y1, Y2 (direct lattice) and b1,b2 (reciprocal lattice) are as defined in the Table below. The vector basis of the reciprocal lattice (also periodical) in the wavevector space whose basis vectors * j e are given by:

FM

  are respectively the normal and transverse forces expressed at the extremity node E (efforts at the origin node O are opposite), are the bending moments about the y and z axis.The vectors ,, uu and  in (3.2) through (3.10) are the displacement, the velocity, the rotation and the rotation rate respectively, and e  ,  are the viscosity coefficients in extension and shear.

  notations, vectors and tensors are denoted using boldface symbols. Those expressions of the forces and moments in previous expressions (3.2) through (3.10) are next involved in the homogenization of the initially discrete lattice towards an equivalent micropolar continuum.

and 1 J

 1 are the components of the rotary acceleration vector. The effective density and micro-inertia therein are given in are, respectively, the mass of the set of lattice beams and the inertia of the cell, A cell being the area of the periodic cell.

  The eigenvalue problem of Eq. (3.21) yields a characteristic equation of order 12 which has 3 of their conjugate solutions representing the 3 modes (illustrated in Fig.3.1): the longitudinal mode (designated by L), the vertical shear mode (labeled SV), the horizontal shear mode (SH in short), and two arbitrary modes due to the presence of the rotation degree of freedom. In the Cauchy medium, we have only 3 modes of propagating (L, SH and SV).

Fig. 3 . 1

 31 Fig.3.1 Schematic diagram showing (from left to right) the wave propagation modes: the longitudinal mode, the shear vertical mode, the horizontal shear mode.

  22) in which s represents the branch number. In equation(3.22), one identifies the natural frequency factor s  , viz, the following quantities

Fig. 3 . 2

 32 Fig.3.2 Elementary unit cell of, (a) 3D hexagonal lattice (left column), (b) plain weave fabric (right column).

  3.3 the effective computed retardation times for the homogenized medium versus the corresponding microscopic retardation time of the beams, in case of uniaxial tension in x and z directions for both structures. They are evaluated for each deformation model as the ratio of the effective viscosity to the effective elastic coefficient from (3.16).

Fig. 3 . 3

 33 Fig.3.3 Retardation times in traction (a) in x-direction for the 3D hexagonal lattice, (b) in zdirection for the 3D hexagonal lattice, (c) in x-direction for the plain weave, (d) in z-direction for the plain weave. The inspection of Fig.3.3 shows that the retardation time for the homogenized medium depends on the retardation time for each beam for the different solicitations. For the two examples and for any type of solicitations, the homogenized retardation time and the beam retardation times are proportional. The same results are found for the case of traction in y-direction and for in-plane xy shear (results not shown here).

Fig. 3 . 4

 34 Fig.3.4 Dispersion relation and damping ratio evolutions for the 3 lattices with a direction of propagation θ=π/6, =π/6 for the 5 following modes. (a) A1μ arbitrary wave, (b) Aβμ arbitrary wave, (c) longitudinal wave, (d) vertical shear wave and (e) horizontal shear wave.The acronyms H, R and T used here and in the sequel are abbreviations of the terms regular hexagonal, re-entrant lattices and plain wave fabric respectively. Whatever the direction of propagation, we see that the location of the branches in the frequency band structure drops with the viscosity e  , while the damping ratio increases with viscosity. It can be seen in Fig.3.4(c)

Fig. 3 . 5

 35 Fig.3.5 Dispersion relations evolutions with a direction of propagation θ=π/6, =π/β for the plain weave textile structure. The disappearance of the horizontal shear wave in Fig. 3.5 clearly demonstrates that the plain weave textile only shows an in-plane wave propagation behavior. The band structure evaluated for the three structures is shown in Fig.3.6 with a direction of propagation θ=π/6, =π/6, and for an in-plane propagation for the plain weave textile.

Fig. 3 . 6

 36 Fig. 3.6 Frequency band structure, in a damped medium with e  =100. (a) regular hexagonal lattice at direction θ=π/6, =π/6, (b) Textile lattice at direction θ=π/6, =π/6, (c) re-entrant lattice at direction θ=π/6, =π/6 and (d) Textile lattice at direction θ=π/6, =π/β.As shown in Fig.3.6, a non-complete band gap occurs for the hexagonal lattice between the different branches; the excess in partial band gap exhibited by the re-entrant hexagonal lattice shows special acoustics properties for this auxetic metamaterial. The complete band gap occurring in the textile structure for a direction of propagation θ=π/6, =π/6 highlights the possibility to generate new fabric materials with significantly better acoustics properties. This

Fig. 3 . 7

 37 Fig.3.7 Frequency band structure for a variation of retardation times. (a) A2 arbitrary mode for hexagonal lattice, (b) A1 arbitrary mode for hexagonal lattice, (c) A2 arbitrary mode for Textile, (b) A1 arbitrary mode for the plain weave textile.

  k  is called anomalous, and in the opposite case ( ) ( ) p g s s c k c k  , the media is called normal dispersive. from the dynamical equilibrium equation and the definition of the group velocity. We present in Fig.3.8 the modulus of the phase velocity for the 5 branches shown in Fig. 3.4; 2 values of viscosity e 200  and e0  are chosen and different directions of the wavevector function of the wavenumber k are considered.

Fig. 3 . 8

 38 Fig. 3.8 Phase velocities for the five wave modes, for e0  and e 200  , for the three investigated examples. (a)regular hexagonal lattice for θ=π/6, =π/6, (b)re-entrant lattice for θ=π/6, =π/6, (c)Textile for θ=π/6, =π/6, d)Textile for θ=π/6, =π/2.

  the k-space, with damping   e 200  and without damping   e0  , for the three modes of wave propagation (L, SV, SH)

Fig. 3 . 9

 39 Fig. 3.9 Group velocity for the 3D hexagone (in x-direction)

Fig. 3 .Fig. 3 .

 33 Fig. 3.10 Group velocity for the 3D hexagone (in y-direction)

Fig. 3 .

 3 Fig. 3.12 Transmission loss for the periodic lattices with the 2 arbitrary modes of propagation. (a) elastic medium, (b) viscoelastic medium e  =200.

Fig. 3 .

 3 Fig.3.13 Unit cell for the regular hexagonal 3D lattice.

Fig. 3 .

 3 Fig.3.15 the frequency band structure by applying Bloch's theorem and the homogenized method, at the edge of the first Brillouin zone. It can be shown how the three branches starting from the origin O in the homogenized theory coincide with those obtained from Bloch theory in the long wavelength limit (O-A-L-H-A region) and for a frequency equal to 1 rad/sec.

Fig. 3 .

 3 Fig.3.15 Comparison of the frequency band structure between the generalized Bloch theorem (blue line) and homogenized theory (red line) along the edge of the first Brillouin zone.

  with the two indices1 and 2 denoting the origin and extremity node of each beam respectively, and   ii u , v the two components of the displacement field in a planar situation. Here and in the sequel vectors and tensors are denoted using boldface symbols. Each beam within the lattice works in traction-compression under the action of normal forces N O and N E (respectively exerted at the two extremity nodes), and in flexion under the action of the transverse efforts T O , T E and moments M O , M E , as pictured on Fig. 4.1.

Fig. 4 . 1

 41 Fig. 4.1 Kinematic and static parameters of a beam element. Based on the equations of beam theory, one shall first express the efforts and moments accounting for the second order displacement gradient, acting at the extremity node of a viscoelastic beam in the Cartesian reference basis, as    

1 and 2

 2 respectively refer to the origin and extremity node of each beam, L εb the beam length,  u and  u are the displacement and the displacement velocity vectors, respectively, e  the extensional viscosity and b I  the quadratic moment of the beam. The resultant and moment expressed in (4.1), (4.2), and (4.3) depend upon the relative displacement and velocity of the extremity nodes of each beam within the lattice.

  of the bending and stretching stiffness can be expressed versus the slenderness parameter be the same for all lattice beams.

  the considered axis e 1 or e 2 , and i  the shift factor (equal to 1  ) for nodes belonging to a neighboring cell, and nil for nodes located inside the considered cell.

  b the beam length and w b the virtual rotation rate, selected as the average of the rotation rates of both beam extremity nodes, and assumed to be regular enough:

  ijk pq pqr pq pqr S successively the stress and hyperstress tensors, and their conjugated kinematic quantities, namely the first and second displacement gradients and their time derivatives, the first and second displacement velocity gradients. The constitutive tensors ,

4 . 2

 42 can further be simplified basing on symmetry properties of the studied lattices: it has indeed been shown[54] that for centro-symmetrical lattices the pseudo-tensors e B   and v B   vanishes. This leads to the following important simplified writing Evaluation of the characteristic lengths of the effective second gradient continuumIn order to identify the internal lengths associated to the different deformation modes (extension and in-plane shear in the present 2D case), we rely on the format of the elastic strain energy density, including stress and hyperstress contributions factor out the successive components of the small strain first order tensor: considering for instance the first stress components 11 e , its expression versus the similar then identify from previous term the internal length in extension in the first direction, the square of which being defined as the ratio of the second order modulus to the classical first order modulus, viz manner, one elaborates the shear energy contribution ( 12 is a first order shear modulus)

For a harmonic wave 12 [

 12 propagating along an axis in an infinite planar second gradient medium, the generalized displacement field with components , UV at a point r is assumed in the following form[57] (each component is a complex number), and a complex frequency function (function of the wave vector) that permits wave attenuation in time. In the limiting case of no damping, it holds =±iω, and the usual form of the plane wave is recovered. Substituting equations (4.47) and(4.48) in the equation of motion(4.45) and (4.46) delivers the following algebraic equation

  For a plane wave without attenuation in the x-y plane, the propagation constants along the x and y kk obtained by solving the eigenvalues problem in (4.49) represents plane waves propagating at the frequency .

Fig. 4 . 2

 42 Fig. 4.2 Representative unit cell of the investigated lattices: (a) Hexagonal lattice, (b) pantograph (inextensible beams),(c) Diamond chiral lattice, (d) re-entrant lattice (θ<0).

  One obtains a vanishing Cauchy stress: 0 e v   , so that elasticity only exists for the second order contribution. The elastic hyperstress tensor is obtained in terms of its independent components as

Fig. 4 . 3

 43 Fig. 4.3 Damping ratio for two values of the damping coefficient µe=10 and =20 e  for a) the diamond chiral lattice for θ=π/4ν b) the diamond chiral lattice for θ=π/6νc) the re-entrant lattice for θ=π/4ν d) the re-entrant lattice for θ=π/6. We plot in Fig. 4.3 the damping ratio  versus the wave number for different viscosity coefficients e  . Fig. 4.4 shows the frequency band structure of the diamond chiral and re-entrant

Fig. 4 . 4

 44 Fig. 4.4 Dispersion relation for low and high damping situations for the a) diamond chiral lattice for θ=π/4ν b) the re-entrant lattice for θ=π/4ν c) the diamond chiral lattice for θ=π/6 and d) the reentrant lattice for θ=π/6 .Results show shifts in the frequency band diagrams (the damping frequency decreases) due to the presence of damping; these shifts are more important with an increase of the viscosity coefficient for the longitudinal and shear modes. This behavior has also an impact on the damping ratio diagrams: when the viscosity coefficient increases, the damping ratio increases proportionally to achieve a maximum value of unity for the 2 modes. These results are in very good agreement

Fig. 4 . 7

 47 Fig. 4.7 Comparison of the dispersion relation for the 4 lattices in a damped medium (µe=20) for : a) Longitudinal waves and b) Shear waves. Results in Fig. 4.5 and Fig. 4.6 indicate that the diamond chiral and pantograph lattices present low natural frequencies in both longitudinal and shear modes due to the lower intensity of the (first order) stress and hyperstress components. The same observation can be made for the re-

Fig. 4 . 6 .

 46 Fig. 4.6. It can be seen in Fig. 4.6 that the longitudinal branches for the four studied lattices are more sensitive to the viscosity coefficient and presents a higher value of damping, in comparison to the shear mode. This can be explained based on the first and second order homogenized viscosity matrix: we find that the stress and hyperstress components coefficients are greater for the longitudinal mode in comparison to the transverse mode ( 11

Fig. 4 . 8 Fig. 4 . 9 Fig. 4 .

 48494 Fig. 4.8 Dispersion relationsversus θ for a) longitudinal waves for the diamond chiral lattice, b) longitudinal waves for the re-entrant lattice, c) shear waves for the diamond chiral lattice and d) shear waves for the re-entrant lattice.

Fig. 4 .

 4 Fig. 4.11 Dispersion relation for longitudinal and shear wave for the Cauchy and second grade continua: a) diamond chiral lattice; b) re-entrant lattice.

  and k the modulus of the wave vector. The dissipative phase velocity is given as follows: the damping frequency. Let recall that the secant slope of the line connecting the origin point with the point of interest on the curve in the dispersion diagram gives the phase velocity, while the tangent to the dispersion curve at any point gives the group velocity. The results in Figs. 4.12 and 4.13 show shifts in the phase velocity due to presence of damping; as the viscosity coefficient increases, these shifts are more pronounced.

Fig. 4 .

 4 Fig. 4.12 Phase velocity for the longitudinal and shear wave, for the diamond chiral lattice for a) θ=π/4 and b) θ=π/6.

Fig. 4 .Fig. 4 .

 44 Fig. 4.13 Phase velocities for the two wave modes of the pantograph lattice

  lattice beams have a viscoelastic behavior described by Kelvin-Voigt model and the homogenized second gradient viscoelasticity continuum model which has first and second order elasticity coefficients reflecting both the initial lattice topology, anisotropy and microstructural features in terms of geometrical and micromechanical parameters. The continuum models enriched with the higher-order gradients of the displacement and velocity introduce characteristic lengths parameters which account for microstructural effects at the mesoscopic homogenized level. The dynamical equations of motion for the equivalent second order continuum have been written to analyze the wave propagation characteristics of four different lattices (the chiral diamond, the classical and reentrant lattices, and the pantograph). A comparative study of the dispersion relations and damping ratio evolutions for the longitudinal and shear waves has been done. The developed model allows analyzing both the effects of damping and internal length scale through the second order gradients on the wave propagation characteristics. The developed homogenization method is valid in the range of low frequencies.

Fig. 5 . 1

 51 Fig. 5.1 Typical window of analysis for a random fibrous networkThe mutual interactions of fibers provide the network connectivity; its non-uniformity can be captured by the fiber number density N (defined as the number of fiber centers per unit area) or rather the contour density of the network, scalar quantity

22 AFig. 5 . 2 Fig. 5 .

 22525 Fig. 5.2 Variation of classical elastic moduli with normalized internal length 0 b lL for a constant density

Fig. 5 .

 5 Fig. 5.4 a) Variation of the anisotropy measure versus window size for the classical moduli and for the b) couple stresses moduli

  1 the mass of the fibrous microstructure, and WOA A the area of the WOA. Note that we have the reduced the set of initially three dynamical equations written in (5.3) to a set of two truly independent dynamical equations.

  in the left-hand side of Eq. (5.7) vanishes; the obtained positive roots obtained characterize the dispersion relations for planar wave propagation. Two modes of wave propagation exist, namely the longitudinal mode (designated by L), and the shear mode (labeled S). The structure of the coefficients of the wave motion matrix in(5.7) shows different powers of the wave vector components; by comparison, for the Cauchy medium, the frequency and the wavenumber only have the same quadratic powers in the wave motion equation, thus the medium is non dispersive.



  the Young modulus, density of fibers and window size respectively.

Fig. 5 .

 5 5 displays the band structure of the random fibrous medium with normalized density 100, for different values of l b . Obviously, an increase of the partial band gap with frequency occurs between the two modes, until it becomes constant when moving from the non-affine deformation regime (at small l b ) to the affine deformation regime (at large l b , when the macroscopic deformation becomes very close to the microscopic deformation).

Fig. 5 . 5

 55 Fig. 5.5 Frequency band structure for the random fibrous medium versus wavenumber and internal length for an effective density , for the longitudinal and shear modes. Continuous line: propagation for the longitudinal mode. Dashed line: shear mode We also observe in Fig. 5.5 the stabilization of the shear mode beyond the value 0.8 kL  in the affine regime (when

Fig. 5 .

 5 Fig.5.6 illustrates the evolution of the modulus of phase velocities in the form of polar plots, for three different values of internal bending length b l in order to highlight the effect of this parameter on the anisotropic dynamic behavior and the dispersive characteristics of the network. Two different values of the wave number are selected in Fig.5.6 in the analysis of the dispersive behavior of the couple stress medium.

Fig. 5 . 6

 56 Fig. 5.6 Modulus of the phase velocities for three values of internal bending length for an effective density 0 50 DL  , for the two modes of propagation, a) 0 / 20 b lL  , b)

Fig. 5 . 7

 57 Fig. 5.7 Frequency band structure versus wave number for different network densities a) Non affine regime, b) Affine regime.

Fig. 5 . 8 Fig. 5 . 9

 5859 Fig. 5.8 Modulus of phase velocities for three values of density in the affine regime ( 0 0 / .02 b lL  ) for the two modes of propagation, a)

Fig. 5 .

 5 Fig. 5.10 Frequency band structure versus window size for two directions of propagation, a) Non affine regime b) Affine regime The obtained results (Fig.5.10) show an increase of frequency with window size for both affine and non-affine regimes; the width of the partial band gap between modes increases with frequency. An increase in frequency occurs due to an increase in the effective couple stress modulus   11 22 , DD density when increasing window size (no change of the classical properties

Fig. 5 . 4 LLFig. 5 .

 545 Fig. 5.11 Modulus of the phase and group velocities for three values of the window size in the affine regime 0

  therein the classical fourth-order elastic tensor, ijklmn D the sixth-order tensor of elastic moduli, and ijlmn M the fifth-order coupling tensor between the first and second order elastic responses which does not vanish for non-centrosymmetric microstructures.For microstructures exhibiting central symmetry -which we assume here and in the sequel -, the fifth-order coupling elastic stiffness tensor M vanishes, so that the previous constitutive law takes

Fig. 5 .

 5 Fig. 5.13 Variation of internal length versus WOA size

Fig. 5 .

 5 Fig. 5.15 Modulus of the phase and group velocities for random fibrous medium for two different values of the wavenumber in the SG medium, a) Non-affine regime, b) Affine regime. Red line:

Fig. 5 .

 5 Fig. 5.16 Effect of the window size and density on the frequency band structure, a) and c) non affine regime, b) and d) affine regime for wavenumber k=0.5. Inspection of Fig.5.16 shows the same influence of density and window size on the dispersion relation occurs for the second gradient and couple stress models in both regimes; a significant increase in frequency for L and S modes occurs when the density or the window size increases when moving from couple stress to second gradient effective continuum. The dispersive behavior of the second gradient medium and the anisotropic characteristic of the random fibrous mediumdo not change when varying the density or the size of the window of analysis. It is important to note that the second gradient effective medium is dispersive for the longitudinal mode which is not the case for the effective couple stress model (as mentioned in section 3).

Fig. 5 .

 5 Fig. 5.17 Polar plot of the modulus of the phase velocity for random fibrous media with a wavenumber k=1. a) Affine regime. b) Non affine regime, Continuous and dashed lines: phase velocity for D=100 and 0 150 DL  respectively. c) Affine regime. d) Non affine regime. Continuous and dashed lines: phase velocity for

6. 1

 1 OverviewStructures having a periodical distribution of their geometry and material properties of their constituents present interesting wave propagation properties, like the existence of frequency band gaps, local resonances, responses directionality due to their anisotropy, left-handedness, cloaking, or negative acoustic refraction. These unusual acoustic properties are due to material and structural heterogeneities associated to periodic modulations of the stiffness and inertial properties, resulting e.g. from modifications of the microstructural configuration. Moreover, the field of acoustic metamaterials has raised a considerable interest due to the possibility to tailor their microstructure to obtain various interesting effects like local resonances, partial or full band gaps, or cloaking, see[1] and references therein. Soft metamaterials have the capability to sustain large deformations, and as a consequence offer promising opportunities of adjusting the acoustic characteristics through the deformation.

  b l  , to a macroscopic lattice length constitutes a small parameter versus which all geometrical and kinematic variables will be expanded, hence b lL    . More details relative to the lattice parameterization are given in the chapter 1 (case of neglecting the rotation and the viscous part).

Fig. 6 . 1

 61 Fig.6.1 Set of repetitive lattices parameterized by a small parameterMaintaining the reference area or volume fixed one considers the limit situation of a continuous density of unit cells obtained when the small parameter tends to zero. In this limit, a continuum, equivalent in a certain sense to the initial lattice, is obtained. To obtain this limit behavior, one does mathematically study the equilibrium of the lattice and the dependence of the governing equations versus the introduced small parameter. Asymptotic expansions of the nodal position, tensions and external forces are written and inserted in the equilibrium equations, preferably expressed in weak form. Taylor series expansion of the displacements and possibly rotational degrees of freedom are next inserted into these equilibrium equations. The discrete sums are finally converted in the limit of a continuous density of beams into Riemann integrals, thereby highlighting continuous stress and strain measures.

  b l and unit beam director b e , both resulting from the beam vector b b b l  B e , successively computed as bb l  B and b b b l  B e

  reference unit cell. The mesoscopic incremental equilibrium equations are next

Fig. 6 . 3 Fig. 6 . 4

 6364 Fig.6.3 Computation of the incremental Cauchy stress A dedicated code has been constructed from the proposed algorithm, in order to solve for the nodal kinematical unknowns (displacements of each beam) within the repetitive unit cell. The code uses an input file the reference unit cell topology and micromechanical properties, and delivers as an output the homogenized mechanical properties (classical moduli and Poisson's ratio) and the nonlinear stress-strain response for a given deformation path imposed over the RUC.In view of the analysis of nonlinear wave propagation, we first express the increment of Kirchhoff stress

1 O

 1  is similar to the linear kernel of order   0 O  . Cancelling the secular terms (the terms multiplied by the effective incremental Young's modulus n E :

Fig. 6 . 5

 65 Fig.6.5 Three studied repetitive lattices. a) 2D hexagonal, b) Milton lattice, c) Hexachiral lattice

Fig. 6 . 6

 66 Fig.6.6 Dependency of the phase velocity on the deformation. Comparison between the incremental scheme and Parnell method [23] for a beam under tension incorporating a hexagonal repetitive microstructure   sec

Fig. 6 . 8

 68 Fig.6.8 Phase velocity under traction versus the applied tensile deformation for the hexagonal, reentrant and Milton lattices. Higher values of the phase velocity are obtained for the re-entrant configuration, while the lowest values occur for the Milton lattice; such effects of the unit cell topology are essentially attributed to the (tangent) rigidity matrix and to the effective density of the considered microstructure.

  vectors condensing the three independent components of the Cauchy stress tensor and small strain tensor built as the symmetrical part of the gradient of the displacement field, vector u .Introducing as for the 1D case the perturbation parameter  in the relation ˆtt nn cc   , which defines a perturbed modulus ˆt n c , leads to the incremental equation of motion in tensor format as

  the frequencies for the longitudinal and shear modes respectively, n  is the frequency of the increment n and the coefficients ij a are the components of the incremental stiffness matrix. The expression of the phase velocity is given by p c k   , where : k  k is the modulus of the wavevector written as the complex number 12 k ik  k ; the real part 1 k represents the attenuation in the x-y plane, and 2 k is the imaginary part of the phase constant.

Fig. 6 .Fig. 6 .

 66 Fig. 6.10 Phase velocity in the longitudinal (red) and shear (blue) modes for the hexagonal lattice for different compression loads (here measured positively).a) 0 yy E  , b) 0.25 yy E  , c)

Fig. 6 .

 6 Fig. 6.13 Phase velocity for the hexachiral lattice for the longitudinal and shear modes for three different biaxial load levels,a) 0 xx yy E E  , b) 0.1 xx yy EE , c) 0.3 xx yy EE .

Fig. 6 .Fig. 6 . 6 . 5

 6665 Fig.6.14 illustrates that the hexachiral lattice shows a nearly isotropic behavior for both the longitudinal and shear modes according to the corresponding circular phase velocity plot in the initial configuration (no deformation is applied). An anisotropic behavior (small degree of anisotropy) appears when applying an axial deformation. The lower degree of anisotropy in comparison to the hexagonal lattice (Fig.6.11) is due to the smaller modification of the hexachiral lattice geometry under an axial deformation.



  5.2) where l  , SV  and SH  are the longitudinal, vertical shear and horizontal shear modes respectively for the new configuration, *  and * 0 are the density of the new and initial contribution, and ij a are the coefficients of the tangent stiffness matrix resulting from the imposed incremental gradient of deformation.

Fig. 6 .

 6 Fig.6.15 Incremental dispersion relation under an incremental axial load for the 3D hexagonal lattice. The solid line corresponds to the longitudinal mode, the dashed line to the vertical shear mode and the dotted line to the horizontal shear mode.

Fig. 6 .

 6 Fig. 6.16 phase velocity for different axial load for the 3D hexagonal lattice.a) 0 xx E  , b) 0.05 xx E  , c) 0.15 xx E  , d) 0.175 xx E  . The anisotropy of the 3D hexagonal lattice in the initial state ( 0 xx E  ) is due to the nonsymmetrical geometry (different dimensions of the unit cell along directions x, y and z), this

Fig. 6 .

 6 Fig.6.17 proposed geometric description of the investigated 3D lattices: (a) 3D re-entrant lattice proposed and (b) pyramid shaped unit cell proposed by Zheng et al., 2011.

Fig. 6 .

 6 Fig. 6.18 Incremental dispersion relation under an incremental axial load (x direction) for the pyramid lattice for different values of the wave number. a) Propagation in direction χ=π/6νӨ=π/6, b) in plane propagation χ=π/βνӨ=π/4.

Fig. 6 .

 6 Fig.6.19 Incremental dispersion relation under an incremental axial load (x direction) for the 3D re-entrant lattice, a) direction of propagation χ=π/6νӨ=π/6, b) in plane propagation χ=π/βνӨ=π/4.

Fig. 6 .

 6 Fig.6.20 and Fig.6.21 illustrate the anisotropic behavior and its evolution for the two auxetic structures (pyramids and 3D re-entrant): the degree of anisotropy increases when increasing the level of the applied deformation in extension.

Fig. 6 .

 6 Fig. 6.20 phase velocity for different axial load for the pyramids structure,) 0 xx E  , b) 0.05 xx E  , c) 0.1 xx E  , d) 0.175 xx E  .

Fig. 6 .

 6 Fig. 6.21 Phase velocity in the longitudinal mode for different axial loads for the 3D re-entrant structure.a) 0 xx E  , b) 0.05 xx E  , c) 0.1 xx E  , d) 0.175 xx E  .

Fig. 6 .

 6 Fig.6.22 Influence of the Poisson ratio and density on the frequency band structure for the 2D hexagon under compression. Continuous (resp. dashed) lines corresponds to edge A of Brillouin zone (resp. edge B). Blue lineμ evolution of density versus Poisson's ratio.When the hexagonal passing from its regular form to re-entrant configuration (passing from positive values of poison ratio to negative values), the density of the material increases, while a shift in the shear mode and a rise in the longitudinal mode occur, providing an increase in the width of partial band gap.

Fig. 6 .

 6 Fig.6.23 Influence of Poisson's ratio and density on the frequency band structure for the γD pyramid under compression. Continuous line corresponds to edge A of Brillouin zone while the dashed line corresponds to the edge B. Results on Fig.6.23 indicate significant effects of a change of Poisson's ratio on the three modes of propagation: a small rise in the 3 modes is observed, while a significant shift in the density occurs.

1 O  equation ( 6 

 16 . The index i therein and in subsequent relations represents the i th branch of the dispersion relation. The corrected frequency 1 n

Fig. 6 .Fig. 6 .

 66 Fig.6.25 depicts the incremental dispersion diagram relying on the perturbation method, taking into consideration the geometrical nonlinearity ( () q NL K), for low (dotted line) and high (dashed line) amplitude wave excitations. In comparison, the solid line corresponds to the absence of geometrical nonlinearity. In the longitudinal mode, a small shift occurs in the dispersion relation as the wave amplitude increases, whereas the shift is more pronounced for the shear mode. For the two modes, the effect of the geometrical nonlinearity becomes weaker when an incremental compressive loading is applied to the structure; this behavior is caused by the inverse

Fig. 6 .

 6 Fig.6.26 Incremental dispersion relation. Low wave amplitude (left) and high wave amplitude (right). Longitudinal mode (solid lines) and shear mode (dashed line). Red line: initial configuration. Green line: 30% deformation gradient in y-direction. Blue line: gradient of deformation of 50% in y-direction.

  propagating through the network, as the lattice configuration changes under the imposed kinematic loading over the unit cell. Especially, the frequency of propagating waves depends on the geometrically nonlinear stiffness matrix of the nonlinear network. The influence of the state of nonlinear deformation within the effective continuum is reflected by the dependency of the angular frequency and phase velocity versus the tangent modulus associated to the new lattice ), we have introduced the projection operators P and C expressing as

7. 2 Fig. 7 . 1

 271 Fig.7.1 Kinematics of a beam element within the repetitive network.     22 2 22 .. .1 2 2 2 bb b cc z F ES EI L L L 



  the rotation of the central node of the beam, and , bb  ee are respectively the beam director and the transverse unit vector. The 2D displacement vector U in expressions (7.1) through(7.3) is parameterized by a curvilinear coordinate denoted s= (s 1 , s 2 ), with components ( any quantity  

L

  the beam length. Expression (7.4) include extensional, flexural and shears contributions to the energy.

Fig. 7 . 2

 72 Fig.7.2 Unit cell of the 2D hexagonal lattice The unit cell of the 2D hexagonal lattice has three beams shown on Fig. 7.2, linking together the pair of nodes (i, j-1)-(i, j), (i, j)-(i-1, j+1) and (i, j)-(i+1, j+1), have shear rigidities denoted

2 Ws

 2 are the first and second order contributions of the strain energy density in small strains, and 3 Ws is the contribution of the strain energy density accounting for the occurrence of large strains.The strain energy density entails the first order Piola-Kirchhoff stress in component form as follows:

1 M

 1 the mass of the set of lattice beams within the unit cell and cell A the area of the periodic unit cell in 2D.

Fig. 7 . 3 7 . 3

 7373 Fig.7.3 Dispersion relation for the hexagonal structure in the linear case (low frequency). Blue lineμ longitudinal mode. Red lineμ shear mode. Comparison between Bloch's theorem (dashed line) and homogenization theory (solid line) for both longitudinal and shear modes

  group velocities respectively, with parameters , E  and L therein the Young modulus, the density and length of the beam respectively.

Fig. 7 . 5

 75 Fig. 7.5 Group velocities for different values of s for re-entrant hexagon with

Fig. 7 . 6

 76 Fig.7.6 Variation of the phase velocity for different values of s.

.

  Nontrivial solutions of the equation of motion exist (inserting the displacements from (7.17) into the equation of motion) if the determinant of the obtained equation vanishes; the positive roots characterize the dispersion relations for planar wave propagation for different values of s. Two modes of propagation exist, namely the longitudinal mode (designated by L) and the shear mode (labeled T). The direction of wave propagation denotes the orientation of the wave vector with respect to the x axis in the Cartesian basis.

Fig. 7 . 7

 77 Fig.7.7 Frequency band structure versus angle, for the supersonic mode s=0.1 and subsonic mode s=0.99, for the longitudinal and shear modes. a) k=1 -direction of propagation equal to π/4, b) k=1 -direction of propagation equal to π/β.We observe in Fig.7.7 that  has a significant effect on the frequency band structure: the transition from negative values of  (re-entrant case) to positive values leads to an important increase in the frequency (for the two directions of wave propagation) for both subsonic and supersonic modes, and the longitudinal and shear modes move apart from each other. Inspections of Fig.7.7 shows that negative values of  are recommended for the hexagonal lattice, since both longitudinal and shear modes coincide for such auxetic lattices. The difference between the longitudinal and shear modes results from the difference between their deformation energy.By varying parameter s from weak to high nonlinearity, one can show that the stress and hyperstress components coefficients are always greater for the longitudinal mode in comparison to the transverse mode due to the inequalities ( 11

Fig. 7 . 8

 78 Fig.7.8 Frequency band structure versus the angle and the slenderness ratio for a direction of propagation equal to π/β and a wavenumber k =1, a) Supersonic mode s=0.1, b) Subsonic mode s=0.8. Continuous line: longitudinal mode; dashed line: shear mode. We observe from Fig.7.8 that the frequency increases when passing from the re-entrant configuration to the regular one for the specific angle 16   

Fig. 7 . 9 .

 79 Fig. 7.9 Frequency band structure versus the degree of nonlinearity s, for a direction of propagation π/4 and a wavenumber k =1μa) dashed line: 6   

Fig. 7 .

 7 Fig. 7.10 Frequency band structure versus wave amplitude, for a direction of propagation π/β and a wavenumber k =1, a) subsonic mode, s=0,999, b) supersonic mode, s=0.1. Continuous line: longitudinal mode -dashed line: shear mode.Green line: 0 5 A  -Red line: 0 1. A  Results in Fig.7.10.a and Fig.7.10.b for different values of s demonstrate that the dispersion curves rise when increasing the wave amplitude, for both longitudinal and shear modes. These

Fig. 7 ..

 7 Fig. 7.11 Phase velocity of the hexagonal lattice in the linear effective medium for different configurations, a) 6    , b) 0   and c)

Fig. 7 .

 7 Fig. 7.12 Phase velocity for the hexagonal lattice in the nonlinear effective medium with,a) 6    and s=0.1, b) 6   

  next made from the local basis (attached to the beam) to the Cartesian basis: the local components (A.10), one obtains the energy density expressed in Cartesian

6 )

 6 in which case the strong ellipticity condition is not satisfied, despite that the second order derivative in (C.4) is always positive; this implies the instability of the propagation of wave and the appearance of a vanishing propagation speed . Previous loss of ellipticity condition entails the existence of a specific wave number c k associated to the vanishing frequency, satisfying the condition function of frequency and wave number   , fk  such that ellipticity is lost for the critical value of the transformation gradient c F . This accordingly explains the vanishing of subsonic modes beyond certain values of the wave number k, observed on Fig.7.4. Note the stabilizing effect of the second order gradient energy density when combined with the first order gradient contribution. For small values of s (close to 0), the influence of the nonlinearity is weak, and the constitutive law of the structures is dictated by the linear parts only; the previous second order derivatives in (C.1) become (the nonlinear modulus 2 E is neglected in (C.3): Both quantities in (C.5) and (C.6) are accordingly always positive, which demonstrates the propagation of supersonic modes obtained in Fig 7.4 for low values of s (the propagation velocity does not vanish).

Fig. 8 . 1

 81 Fig.8.1 Schematic representation of a 2D plain weave fabric (left) and the chosen unit cell (right). The evaluation of the energy density for any repetitive beam network is based on the summation of the energy contributions associated to each beam within the identified unit cell of the repetitive lattice. It relies on the following steps: a) For each beam b, one writes the expressions of forces and moments exerted on the beam extremities in the nonlinear framework of 3D Bernoulli beams (we neglect the twisting moment), considering a large displacement large rotations framework



  the rotation of the central node of the beam, U a three-dimensional displacement field parameterized by a curvilinear coordinate denoted s = (s 1 ,s 2 ,s 3 ), with components ( , , )  U u v w .

  the curvilinear abscissa along the beam. c) The transition is next made from the local basis attached to the beam (curvilinear coordinates) to the Cartesian global basis; one then obtains the expression of the total deformation energy for the unit cell of the network (for more details, see Appendix A and chapter 7) d) The energy density is finally evaluated by dividing the total energy, elaborated as the sum of the energies of all beams within the unit cell in the Cartesian basis, by the area of the elementary unit cell that is   22  S = 4 a L cos , considering the more specific plain weave unit cell. It accordingly holds in full generality the decomposition of the energy density of the unit cell into three contributions reflecting successively the linear part of the first order elasticity, the linear part of the second order elasticity and the nonlinear part associated to first order elasticity, viz / linear SG nonlinear

. 1 )

 1 In previous equations, u i are the components of the acceleration vector, and the effective density is given in general by 1 the mass of the set of lattice beams, A cell being the area of the periodic unit cell. For the considered plain weave textile structure, the effective density is computed as:

Fig. 8 .

 8 Fig.8.2 a) Direction of propagation in 3D space defined by the two angular variables , b) Schematic diagram showing (from left to right) the wave propagation modes: the longitudinal mode, the shear vertical mode and the horizontal shear mode. In the Eq. (8.1), the first term reflecting the sole effect of Cauchy type elasticity describes a nondispersive (so called quasi-static) behavior of the homogenized (Cauchy) medium. The second term describes the nonlinear physical and geometrical effect; it affects the dispersive behavior of the medium. The introduction of second gradient terms in Eq. (8.1) brings an essential change in the behavior of the medium from non-dispersive to dispersive.

4 )

 4 Fig.8.3 the dispersion relation, phase and group velocities in the longitudinal direction (e 1 direction).

Fig. 8 .

 8 Fig.8.3 a) Dispersion relation, b) phase velocity, c) group velocity, in the longitudinal direction (e1). The red line corresponds to the longitudinal mode; the green line corresponds to vertical shear while the blue one corresponds to horizontal shear mode. One can observe in Fig.8.3a that frequency increases nonlinearly with the wavenumber k, thus demonstrating dispersive effects for the second order effective medium. This feature is further attested by the difference between the phase and group velocities (Fig.8.3b and Fig.8.3c).

Fig. 8 .

 8 Fig.8.4 shows the evolution of the frequency and the phase velocity versus the wavenumber in the two planes (e 1 , e2) and (e 1 , e3). As expected, the frequency and the phase velocity increase with the wavenumber. It has been shown recently that the plain weave textile shows very good inplane mechanical properties, but suffers from some weakness in the out-of-plane direction; this behavior can be clearly observed from the weak values of the horizontal shear mode for waves propagating in the (e 1 ,e 2 ) plane (Fig.8.4a and Fig.8.4b).The dispersion relation for the textile with a direction of propagation in the plane (e 1 , e3) is pictured in Fig.8.4c and Fig.8.4d. One can observe that the horizontal shear mode is higher compared to the longitudinal and vertical shear modes.

Fig. 8 . 4

 84 Fig.8.4 Dispersion relation and b) phase velocity for the textile plain weave structure for wave propagation in (e1,e2) plane, c) Dispersion relation and d) Phase velocity for wave propagation in the plane (e1,e3). The red surface; longitudinal mode; green surface: vertical shear mode; blue surface: horizontal shear mode.

  al kk the flexural rigidities of the beams with length a and L respectively (see Fig.8.1), and , al KK the tensile rigidities of the beams with length a and L respectively. Conversely, neglecting these terms delivers the conventional 3D Cauchy medium equations. Due to nonlinearity terms in (8.5), the harmonic plane Ansatz is not a solution of this nonlinear wave equation, but one shall instead resort to elliptic functions as shown in the sequel.The extension of the 1D equation (8.5) to 3D networks exhibiting deformation modes in different directions relies on the computation of the energy density by the continualization techniques exposed in previous section 2. The three differential equations of motion are similar in form to Eq. (8.5), and they show couplings between longitudinal waves, horizontal shear waves and vertical shear waves, through the kinematics parameters of the rigidity matrices.

Fig. 8 . 5

 85 Fig.8.5 Dispersion relation for different values of s along the direction. Longitudinal (red), Vertical Shear (green) and Horizontal Shear (blue) modes. Solid lines: s=0.95, dotted lines: s=0.999, dashed lines: s=0.1, bold lines: s=0.4.Black line: dispersion relation in the linear case.

Fig. 8 . 6

 86 Fig.8.6 Phase velocity versus wavenumber k for different values of s along the direction. Longitudinal (red), Vertical Shear (green) and Horizontal Shear (blue) modes. Solid lines: s=0.95, dotted lines: s=0.999, dashed lines: s=0.1, bold lines: s=0.4.

Fig. 8 . 7

 87 Fig.8.7 Group velocity versus wavenumber k for different values of s along the direction. Longitudinal (red), Vertical Shear (green) and Horizontal Shear (blue) modes. Solid lines: s=0.95, dotted lines: s=0.999, dashed lines: s=0.1, bold lines: s=0.4. In the sequel, the dispersion relation, the phase and group velocities of the investigated continuum are pictured in Fig.8.8 for in-plane propagation in planes   12 , ee and   13 , ee to

Fig. 8 . 8

 88 Fig.8.8 Dispersion relation, a) supersonic modes for two values of s=0.1 and s=0.4, b) evanescent modes for s=0.95 and s=0.999, for wave propagation in the plane, c) supersonic modes for two values of s=0.1 and s=0.4, d)evanescent modes for s=0.95 and s=0.999 for wave propagation in the   13 , ee plane.

Fig. 8 .

 8 Fig.8.8 reveals that the studied continuum presents two sets of dispersion relations: dispersive modes described by the subsonic modes vanishing beyond a specific wavenumber for the three modes of propagation (L, SV, SH) (Fig.8.8b and Fig.8.8d), and supersonic modes in which the frequency increase by the increasing of the wavenumber (Fig.8.8a and Fig.8.8c).

Fig. 8 . 9 13 ,

 8913 Fig.8.9 Phase velocity versus wavenumber, with different values of s, a) supersonic mode for s=0.1 and s=0.4 with wave propagate in plane, b) subsonic mode for s=0.95 and s=0.999 with wave propagation in plane, c) supersonic mode for s=0.1 and s=0.4 with wave propagation in the   13 , ee plane, d) subsonic mode for s=0.95 and s=0.999 with wave propagation in   The shape of the phase velocity in Fig.8.9 reflects the dispersion relation obtained in Fig.8.8. Two sets of modes are obtained: supersonic modes (Fig.8.9a and Fig.8.9c) and evanescent subsonic modes (Fig.8.9b and Fig.8.9d respectively), corresponding respectively to low and high values of the nonlinearity parameter. It clearly appears that the subsonic modes vanish beyond certain values of the wave numbers 12 k , k in the  

Fig. 8 .Fig. 8 . 44 

 8844 Fig. 8.10 Group velocity versus wavenumber, for different values of s, in plane propagation, a) supersonic modes for 1 g c k     , b) subsonic modes for 1 g c k     c) supersonic modes for

Fig. 8 .

 8 Fig.8.12 a) Dispersion relation, b) phase velocity, c) group velocity, for different values of s in 3D space propagation.

Fig. 8 .

 8 Fig. 8.13 Phase velocity in a 3D linear textile versus with 2    , for a) L mode, b) SV mode, c) SH mode. The continuous line corresponds to k=1, the dashed line to k=3.

Fig. 8 .

 8 Fig. 8.14 Phase velocity versus  with

Fig. 8 .

 8 Fig. 8.16 Phase velocity plots versus  with 2    for k=3 in the supersonic mode for 4 values of the degree of nonlinearity, a) s=0.999, b) s=0.95, c) s=0.4 and d) s=0.1.

  moment of the beam, b b ,  ee respectively the unit director and transverse unit vectors for each beam. The truss under consideration is made of beams and is completely defined by the initial positions of the nodes and their connectivity. Each beam links two nodes and is oriented so that it has an origin node O(b) , an end node E(b) , and a central node C(b) , Fig. 9.1. The shift parameter  traduces the fact that the end node belongs to the next neighbouring cell. Each extremity node has two displacements in the two principal directions and one rotation in plane   i, j , as pictured in Fig.9.1 together with the efforts and moments.

Fig. 9 . 1 2 )

 912 Fig.9.1 Kinematic and static parameters of a lattice beam

   in the present 1D situation, representative of classical elastic Cauchy materials, which are coined Form1 and Form2 here and in the sequel:  Form1: The strain energy density takes the form[24]

6 )

 6 context, all stress measures coincide, so especially it holds the identity   We shall in the sequel consider different suitable forms of the constitutive law written in terms of nonlinear relations between the first Piola-Kirchhoff stress and the transformation gradient. We shall also as a matter of simplification of notations omit the x dependency. Note that the strain energy density can easily be expressed versus the stretch, due to the relationF   .The material parameters are identified based on a combination of virtual tensile test performed over the unit cell of the three considered lattices; their identification proceeds from the minimization of the following function with respect to the set of material parameters A, B,C

Fig. 9 . 2

 92 Fig.9.2 Three investigated lattices: a) classical hexagonal lattice, b) re-entrant hexagonal lattice and c) textile plane weave.

8 )

 8 Inserting the constitutive law into the dynamical equilibrium, Eq.8, leads to in the analysis of the nonlinear dispersion relation in the continuum medium is the introduction of the dimensionless time t   and the parameter to enforce the non-linearity in the dynamical equation, CC   , thus leading to the asymptotic expansion of frequency and axial displacement (here truncated to the first order), successively developments (Eq.9.10) into the nonlinear wave equation and ordering versus the successive powers of the small parameter  produces a set of equation as follows:

1 O 1 O 1 

 111 is the wavenumber and A the wave amplitude. Subsequent substitution of the expression of wave solution into the    term results in the partial differential equation is the corrected frequency based on the nonlinear terms. The frequency is then updated versus the wave amplitude as follows:

Fig. 9 .

 9 Fig.9.3 illustrates the dispersion relation based on Eq.9.15 for the classical and reentrant hexagonal networks and the textile plane weave structure.

Fig. 9 . 3

 93 Fig.9.3 Linear (red) and nonlinear (green) dispersion curves for a) the hexagonal network, b) the re-entrant network and c) the textile plane weave. Dispersion shifts occur for the longitudinal wave through the introduction of the nonlinear parts represented by the corrected frequency 1

Fig. 9 .Fig. 9 . 5

 995 Fig.9.4 (a) Shape of the solitary strain wave and (b) Shock wave The strain amplitude can be calculated from the following equation   2 0 3 () 31 () CA Es Ks k    (9.22)

Fig. 9 .

 9 Fig.9.6 shows the frequency for the supersonic longitudinal mode for the three lattices versus the degree of nonlinearity s.

Fig. 9 . 6

 96 Fig.9.6 Frequency band structure versus the degree of nonlinearity s for the hexagonal lattice (in red), the re-entrant lattice (in green) and textile plane weave (in blue).

( 9 . 26 )

 926 Using the same perturbation method described in 9.4.1 leads to the following equations:

Fig. 9 .

 9 Fig.9.7 illustrates the dispersion relation for the classical and reentrant hexagon networks and the textile plane weave structure with and without second order terms, in order to investigate its impact on the dispersion relation for both the linear and nonlinear media.

Fig. 9 . 7 1  29 ) 30 )

 9712930 Fig.9.7 Linear (red) and nonlinear (green) dispersion curves with (continuous) and without (dashed line) second gradient terms. Dispersion shifts occur for longitudinal waves due to the nonlinear terms introduced in the wave equation, accounted for by the corrected frequency 1 ; this behavior can be observed for the three investigated lattices. The influence of the second gradient term is more pronounced for the hexagonal lattice in both linear and nonlinear cases due to the higher values of  . Furthermore, inspection of Fig.9.7 reveals that the classical continuum is not dispersive since the frequency is proportional to the wavenumber k (dashed red lines).

Fig. 9 . 8

 98 Fig.9.8 Dispersion relation with different values of s based on Boussinesq type equation for a) the hexagonal lattice, b) Re-entrant lattice and c) Textile Plane weave lattice.The dispersion relations are pictured in Fig.9.8 for different values of the nonlinear parameter s for the hexagonal network. In the sequel, parameters v and v 0 refer to the phase velocities in the nonlinear and linear medium respectively.

Fig. 9 .

 9 Fig.9.10 the phase and group velocities in the longitudinal mode for the hexagonal (red line) network; the re-entrant network (blue line) and textile plane weave (green line). In the following, results for the phase velocity (corresponding to the continuous line) and the group velocity (corresponding to the dashed line) are shown.

Fig. 9 .

 9 Fig.9.10 Evolution of the phase and group velocities of the longitudinal mode for the three investigated lattices.

Fig. 9 .

 9 Fig.9.10 shows that the low frequency behavior (for longitudinal mode) for the three investigated lattices is dispersive, due to the noticeable difference between the phase and group velocities in the effective nonlinear medium. Form1 and Form2 of the constitutive law lead to a variation between the phase and group velocities, which become more pronounced as the wave number is increased.

  317

Wave propagation in linear network materials accounting for scale effects 2. Chapter 2: Micropolar dissipative models for the analysis of 2D dispersive waves in periodic lattices Summary The

  computation of the dispersion relations for dissipative periodic lattices having the attributes of metamaterials is an actual research topic raising the interest of researchers in the field of acoustics and wave propagation phenomena. We analyze in this chapter the impact of wave damping on the dispersion features of periodic lattices, which are modeled as beam-lattices. The band diagram structure and damping ratio are computed for different repetitive lattices, based on the homogenized continuum response of the initially discrete lattice architecture, modeled as Kelvin-Voigt viscoelastic beams. Three of these lattices (reentrant hexagonal, chiral diamond, hexachiral lattice) are auxetic metamaterials, showing negative Poisson's ratio. The effective viscoelastic anisotropic continuum behavior of the lattices is first computed in terms of the homogenized stiffness and viscosity matrices, based on the discrete homogenization technique.

7, we analyze nonlinear wave propagation in homogenized 2D nonlinear second order gradient continuum; a Boussinesq type wave equation is formulated, the solutions of which are characterized. Chapter 8isan extension of the method presented in chapter 7 to 3D repetitive networks, selecting 3D textile plane weave as an application of the proposed methodology. We analyze in Chapter 9 the acoustic properties of microstructured beams including a repetitive network undergoing configuration changes leading to geometrical nonlinearities. Finally, we conclude the thesis in Chapter 10 by summarizing the achievements and limitations of the work, and describe several directions with promising opportunities for continued work.

I.

Table 2 -

 2 1 Geometrical and mechanical parameters of the four studied lattices

Type Geometrical parameters of the unit cells Mechanical properties Re-entrant

		l		 50 mm 15 mm t  , 1mm, 		30   	s E	1400MPa, 1000kg/m   3   	0.3,
	Hexachiral				l		 50 mm 15 mm t  ,		1mm	s E	1400MPa, 1000kg/m   3   	0.3,
	Diamond chiral				l		 50 mm 15 mm t  ,		1mm	s E	1400MPa, 1000kg/m   3   	0.3,
	Hexagonal	l		 50 mm 15 mm t  ,		1mm, 	30  	s E	1400MPa, 1000kg/m   3   	0.3,

expressions of forces and moments at beam level

	32 2 00 2 LL -12EI 6EI 32 2 3 2 3 22 EA -EA 00 12EI 6EI 00 LL LL 6EI 4EI -6EI 2EI 00 LL LL -EA EA 00 00 LL 12EI -6EI 12EI -6EI 0 -0 L L L L 6EI 2EI 6EI 4EI 0 0 -L L L L                 	C=	32 2 2 00 LL -12 I 6 I 32 2 3 2 3 22 A A 0 0 -12 I 6 eI 00 LL LL 6 I 4 I 6 I β eI 00 -LL LL eA A 00 -0 0 LL 1 ee ee e ee e e ee 2 I -6 eI 1β eI -6 I 0 -0 L L L L 6 I β I 6 I 4 I 0 0 -L L L e e e e L                 
	Appendix B:

Table 2 -

 2 2 basis vectors of the direct and reciprocal diamond chiral lattices and coordinates of the edge of IBC. L denotes the length of each beam.

	Topology	Direct lattice	Reciprocal lattice	coordinate of the edge of
				IBC
	Diamond chiral	1
	lattice	2

  

	xx x	xy y	xz z	u
	yx x	yy y	yz z	v
	zx x	zy y	zz z	w
	xy mm xx m x y z xz	yz	zy	I		x x
	yx x m	yy y m	yz z m	zx	xz	I	y y
	zx x m	zy y m	zz z m	xy	yx	I	z z

Table 4 -

 4 1 Geometrical and mechanical parameters of the four lattices Type

Geometrical parameters of the unit cells Mechanical properties

	Re-entrant	l	50mm,	t	1mm,	30

Table 4 -

 4 2 Elastic and viscous internal lengths for the studied lattices

Internal length Diamond chiral Hexagonal lattice Reentrant lattice Pantograph

	l	e xx	L	2	2	2	2 2

12 , EE and 3 E in the nonlinear wave equation

	The elastic moduli 12 , EE and 3 E involved in the nonlinear wave propagation equation (7.10) are
	given respectively versus the lattice microstructural parameters after straightforward
	computations as follows:

  3. Compute the incremental Second Piola-Kirchhoff stress tensor[14, 15] Check convergence at iteration k; if it is attained, go to next step.5. Compute the incremental deformation gradient and its Jacobean. 6. Update Cauchy stress at increment ( n1  ) by a push forward of its Lagrangian

							  k n SE   k S T, n Gn K:   
	4. counterpart from n  to n1  
	  k n 1 		J	 F S 1 n n . 	(k) n	  S	  k n	  T n .  F	J	1 n n  F S F (k) n n T . .				J	1 n n  F	.	    k n  S	T n F .	
										( k ) n								  k n
	7. Update the network configuration from n  to n1   .

Table 9 -

 9 1 Coefficients of the three forms of the hyperelastic strain energy potential

			Form 1			Form2
		A	B	C	A	B	C
	Hexagonal	-0.1	40.1	1469.3	0.0604	22.218	316.5386
	Re-entrant	-0.4939 244.775 10220.08	0.0999	156.3594	1896.04
	Plain weave	-0.0929 124.129 330.3384	0.0103	116.0903	101.9276

Table 9 -

 9 2 Coefficients of the second order terms for the three investigated structures

		Hexagonal	Re-entrant	Plain weave
		701.48	25.98	4.4

Calcul des propriétés effectives statiques et dynamiques de matériaux architectures prenant en compte les effets microstructuraux et les grandes déformations'

  Résumé : Nous analysons les propriétés dynamiques de milieux architecturés périodiques et de réseaux fibreux aléatoires en petites et grandes déformations, à partir de méthodes d'homogénéisation dédiées, afin de calculer leurs propriétés statiques et dynamiques. Des modèles effectifs de type micropolaire et du second gradient sont élaborés afin de prendre en compte l'impact de la microstructure sur le comportement effectif. L'influence des degrés de liberté additionnels sen rotation additionnels et des gradients d'ordre supérieur du déplacement sur les relations de dispersion sont analysés pour des comportements élastique et viscoélastique du matériau constitutif. Les milieux continus généralisés ainsi construits conduisent à des effets dispersifs, en accord avec les observations. Dans la seconde partie du travail, nous analysons l'influence des grandes déformations sur la propagation des ondes élastiques dans des milieux architecturés périodiques. Des méthodes théoriques assortis de schémas numériques sont développés afin de prédire l'influence des déformations finies générées au sein des structures sur l'évolution de leur diagramme de bande. Un schéma incrémental d'évolution de la fréquence et de la vitesse de phase du milieu continu homogénéisé est établi, à partir d'une méthode de perturbation établie pour des structures 1D, 2D et 3D, en considérant plus particulièrement des structures auxétiques. Les calculs montrent un effet important de l'état de déformation appliquée et de la densité effective sur l'évolution de la fréquence et de la vitesse de phase des ondes. Une méthode de perturbation spécifique aux structures périodiques nonlinéaires est développée afin de généraliser le théorème de Bloch pour couvrir les non linéarités tant géométriques que matérielles. Des modèles hyperélastiques du premier et du second gradient de différentes structures sont identifiés par des tests virtuels reposant sur une méthode d'homogénéisation dédiée, qui permettent de formuler des équations d'onde spécifiqueséquations de Burgers etde Boussinesqdont les propriétés dispersives sont analysées. Mots clefs: matériaux architecturés ; propagation d'ondes ; méthodes d'homogénéisation ; effets d'échelles ; grandes déformations ; viscoélasticité

Key words: network materials; wave propagation; homogenization methods; scale effects; large deformations; viscoelasticity Titre : '

1000kg/m  

, ee planes and for the three propagations modes (L, SV, SH).

, ee plane.

, ee plane, and for

k , k in   13 , ee plane for higher

2 2

ACKNOWLEDGEMENTS

A thesis is the result of the association of many people. I have written these few lines in order to express my gratitude to them. Nevertheless, the exercise remains perilous because it is difficult to thank absolutely everyone, that the "forgotten" excuse me in advance. This thesis was carried out in Co-direction between Lebanese university and the Laboratory of Energy and Theoretical and Applied Mechanics (LEMTA) at the University of Lorraine. I wish to dedicate my first words to my supervisors. First of all, to Mr. Jean-François

List of Tables

Appendix D: First Brillouin zone for the diamond chiral lattice

The first Brillouin zone is defined as a primitive unit cell of the reciprocal lattice, and it can be constructed as follows:

with those found in [57,58]. The same influence of the viscosity coefficient on the damping ratio and frequency band structure can also be observed for the pantograph and the hexagonal lattices (see Fig. 4.2). The vanishing of certain modes of propagation (for the re-entrant lattice) for a certain wavenumber is due to the total dissipation of energy in the presence of damping. There are several examples in the literature showing the cut-off phenomenon, by which energy completely vanishes: this can be seen e.g. in reference [58] a cut off branch appears when the authors increase the viscosity coefficient. In [57,18] one can see the same cut off phenomenon when the viscosity coefficient is increased; the cut off phenomenon also appears in [59]. The energy decreases exponentially, this decrease depends on the real part of the complex frequency function, which is ( Re( )

; this entails that energy will be dissipated faster when the viscosity coefficient is increased (which increases the damping ratio  as shown in Fig. 4.3). One can also conclude that the damping ratio and the natural frequency (case of non-damping) for the re-entrant lattice has higher values compared to the diamond chiral lattice, thus the diamond chiral lattice can filter frequencies more efficiently than the re-entrant lattice. We also observe for the pantograph lattice that the phase velocity for the longitudinal and shear modes is identical and is not influenced by the direction of wave propagation θ; this follows logically from the obtained dispersion relations (Fig. 4.10). implies that the pantograph will have the best absorption properties amongst the considered lattices since the waves will propagate with lower velocity and will thus be better absorbed. The presented results entail that the pantograph lattice exhibits the best acoustic characteristics.

Effect of internal length on the dispersion relations

Changing the lattice topology and its geometrical parameters leads to a variation of the characteristic lengths associated to the different deformation modes, since the expression of these lengths follow from the computed effective moduli.

We present in Fig. 4.15 the influence of the magnitude of the internal length in extension and shear on the natural frequency, for the diamond and reentrant lattices.

The two main classes of generalized continuum theories are the higher-grade theories, in which the gradients of strains or the higher order gradients of the displacement are incorporated; the higher-order theories incorporate additional degrees of freedom and constitute a second class. The reader is referred to [14] who proposed a historical overview of generalized continuum theories.

The development of the non-linear theory of elasticity returns to the seminal work of Cosserat brothers that did not get the attention it deserved for a long time. At the beginning of the 1960s, prominent authors became interested in Cosserat theories [15][16], and a special case of the Cosserat continuum theory was investigated by [17], in which the rotation of the rigid Cosserat triad is defined in terms of the displacement gradients, deserving the name of couple stress theory. In a recent work, [18][19] identified the couple-stress moduli of vertebral trabecular bone, based on a prototype model for the 3D internal architecture, whereby the network of trabecular is modeled as a porous material with an idealized periodic structure made of 3D open cubic cells;

this was followed by a more recent contribution of the same authors [20] devoted to the identification of strain gradient models for composite materials, including three-dimensional random porous polymer scaffolds, composite reinforced by inclusions, and woven composites.

In this chapter, we will use the couple stress theory and strain gradient models as a modeling framework in order to analyze the size effects at the level of windows of analysis of the random fibrous network.

The chapter is organized as follows: the generation of a random set of fibers in a 2D context within windows of analysis is explained in section 2, together with the method used for the identification of the couple stress moduli based on the equivalence of strain energy. The influence of important parameters characterizing the network on the computed effective mechanical moduli of the couple stress medium is studied in section 3, such as the bending length, the density of fibers, and the size of the window of analysis. In section 4, we present the second gradient homogenized medium of random fibrous network. The acoustic properties of random fibrous networks are analyzed based on strain gradient models as an alternative framework in section 5, which provides a comparative analysis of the same features obtained by couple stress theory in the first part of the chapter. Finally, a summary of the work is exposed in section 6 together with perspectives for future work.

the phase velocity does not depend on the wavenumber, but only on the direction of wave propagation.

We shall plot the dispersion relation by restricting the values of the wave vector to the edges of the irreducible part of the first Brillouin zone (O, A, B).

In order to exemplify the methodology of the incremental dispersion relation, the longitudinal and shear frequencies are identified in a first step for the hexagonal lattice; the response of this lattice is computed under a compression in y direction, in order to highlight the variation of the band gap when passing from a regular configuration to a re-entrant onefor which the structure is auxetic, and under simple shear.

In order to investigate the evolution versus strain of the incremental dispersion relation and of the band gap for the hexagonal structure (changing from a regular to an auxetic configuration under a compressive load), we plot in Fig. 6.9 the dispersion relation for both longitudinal and shear modes under varying compressive and shear loads. A variation of the load (compression or shear) significantly alters the dynamical response of the structure (Fig. 6.9); especially, the frequency of the dynamical response under compression increases significantly, even for small load changes.

Those expansions are then substituted into the equation of motion (5.6) to yield ordered equations

8)

The perturbation method is applied to equation at order  

Summary

The general objective of this chapter is the analysis of wave propagation phenomena within architecture media, relying on an effective substitution continuum obtained by homogenization. The proposed methodology is quite general and applicable to any 2D repetitive network of beamlike structural elements, considering beams undergoing large transformations. Based on the writing of the equations of motion of a nonlinear second order gradient continuum, we analyze the nonlinear wave propagation in the obtained homogenized nonlinear second order gradient continuum. The resulting wave equations are of Boussinesq type, the solution of which being elliptic functions. The influence of the degree of nonlinearity on the dispersion relations is analyzed, highlighting subsonic and supersonic modes propagating respectively with a velocity lower (resp. higher) than the velocity of linear non-dispersive waves. Subsonic and supersonic modes correspond respectively to regimes of high and low nonlinearity characterized by the socalled universal constant s. The existing anisotropy of wave propagation becomes more marked when the degree of linearity increases.

We represent in Fig. For a weak nonlinearity, a supersonic mode occurs and the dispersion curves lie above the linear dispersion curve (v = v 0 ). For a higher nonlinearity, the waves change from a supersonic to an evanescent subsonic mode at approximately s=0.7 and the dispersion curve drops below the linear curve and the frequency vanish for certain values of k (each depending on parameter s) (see Appendix C for more details).

In a nonlinear situation, the second order gradient anisotropic continuum has two propagation modes (subsonic and supersonic), whereas Cauchy or micropolar continua only have supersonic modes, for which an increase of the frequency with the wavenumber occurs.

One can observe from Fig. 7.5 that the group velocity tends to zero when the wave frequency of all evanescent subsonic modes approaches their higher value. The vanishing of mode occurs when g v   , which means that the strong ellipticity condition is not satisfied [44] (instability occurs for high values of s). In the supersonic mode, as the wavenumber k increases, the accumulated energy in the medium (increasing in the group velocity) leads to an increase in the dispersion relation.

Chapter 8: Nonlinear dynamical analysis in 3D textiles based on second order gradient homogenized media

Summary

We expose in this chapter an extension of the method presented in the chapter 7 to 3D repetitive network, we have chosen the 3D textile plane weave as an application, noted that the proposed methodology is quite general and applicable to any 3D repetitive network of beam-like structural elements, considering beams undergoing large transformations. Based on the writing of the equations of motion of a nonlinear second order gradient continuum, we analyze the nonlinear wave propagation in the obtained homogenized nonlinear second order gradient continuum. The resulting wave equations are of Boussinesq type, the solution of which being elliptic functions. The influence of the degree of nonlinearity on the dispersion relations is analyzed, highlighting subsonic and supersonic modes propagating respectively with a velocity lower (resp. higher) than the velocity of linear non-dispersive waves. The three modes of propagation (longitudinal, vertical and horizontal shear) are compared in terms of dispersion relations, phase and group velocity diagrams. The existing anisotropy of wave propagation becomes more marked when the degree of linearity increases. The horizontal and vertical shear modes disappear successively when increasing the wavenumber.

nonlinearities, while the frequency of the supersonic modes increase due to the accumulation of energy with the wavenumber. The difference between the phase velocity considering the different propagation cases and the transition of the phase velocity in the longitudinal mode from higher values to lower values are aspects that have been discussed before.

The evolutions of the group velocity versus the wavenumber in the two planes of propagation are investigated in Figs (8.10 and 8.11). We shall note that the group velocity depends on the phase velocity, so for each propagation modes we have s sets of group velocity. The shape change of the phase velocity when passing from k=1 to k=1.5 involve the dispersive behavior for the second gradient medium. For k=1.5 (Fig. 8.15), we observe the disappearance of the horizontal shear mode for the higher nonlinearity case (s=0,999 and s=0.995), whereas the three modes do coexist for the weak nonlinearity case; this behavior also appears from inspection of Fig. 8.5. We analyze in this chapter the acoustic properties of microstructured beams including a repetitive network material undergoing configuration changes leading to geometrical nonlinearities. The effective constitutive law is evaluated successively as an effective first and second order nonlinear grade 1D continuum, based on a strain driven incremental scheme written over the reference unit cell, taking into account the changes of the lattice geometry. The dynamical equations of motion are next written, leading to specific dispersion relations. The inviscid Burgers equation is obtained as a specific wave propagation equation for the first order effective continuum, while the Boussinesq equations obtained when considering a strain gradient effective continuum. The presence of second gradient order term in the nonlinear equation of motion leads to the presence of two different modes: an evanescent subsonic mode for high nonlinearity that vanishes after certain values of wave number k, and a supersonic mode in case of a weak nonlinearity. This methodology is applied to different microstructures, including the regular and reentrant hexagon, and plain weave.

Appendix A: Expressions of the deformation energy of each beam in large strains

A novel procedure for predicting the effective nonlinear elastic responses of repetitive lattices in the framework of the mechanics of micropolar continua through a combined linear and nonlinear discrete homogenization scheme shall be presented briefly (the reader is refereed to [14,15] for more details). The nonlinear stress-strain response will be computed incrementally for 2D structures subjected to different loading cases (uniaxial, biaxial, simple shear), taking into consideration changes of the structure geometry. The predictive nature of the employed homogenization technique allows the identification of a strain energy density for hyperelastic models of the structure at the mesoscopic level. The identified hyperelastic constitutive models will then be involved in the analysis of wave propagation in repetitive network materials represented by the constructed effective substitution medium.

A lot of attention has been paid to the propagation of elastic of elastic waves in linear framework [16][17][18][19], whereas a few authors analyzed so far wave propagation in nonlinear media. The propagation of elastic waves in nonlinear composite materials is accompanied by a number of new phenomena such as amplitude dependent dispersion relations, or the occurrence of subsonic and supersonic modes that can never be observed in homogeneous linear media [20][21][22][23]. These features entail that the solutions of the wave propagation equations are more complex compared to the linear case, and they depend on the form of the dynamical equilibrium equations derived from the constitutive law. For example, the solitary surface waves discovered by John Scott

Russell [24] in 1834 have been developed as solution of the Boussinesq equation [25], the Korteweg & de Vries (KdV) equation [27], the Benjamin-Bona-Mahony (BBM) equation [26],

the Camassa-Holm (CH) equation [28]. The shock displacement wave can be used for nonlinear dynamical problems in the form of Burger's equationν as an alternative, the perturbation method can be used for certain kinds of constitutive laws [29].

The study of nonlinear elastic waves has been limited so far in the literature to classical Cauchytype elasticity theory, which relies on the sole first order displacement gradient. Models based on Cauchy-type theory do however not give realistic predictions of the medium properties, such as the dispersion relation, since the Cauchy effective medium lacks internal length parameters. The Cauchy medium is non-dispersive, which means that waves propagate independently of the wavenumber. It has however been proven by experiments that most waves are dispersive, that is, each wavenumber travels with a different phase velocity [30][31].

The outline of this chapter is as follows: section 2 is devoted to the description of the discrete homogenization method in a large strains context, at both microscopic and mesoscopic scales.

The homogenization technique and the expression of forces and moments in the framework of 2D

Euler-Bernoulli beams are exposed. The algorithm used for the incremental procedure for the update of variables at the lattice level accounting for the evolution of the network geometry will be described. In section 3, virtual simulations based on the developed discrete homogenization technique will be used for the calibration of a strain energy density of a hyperelastic model for three different lattices, leading to two different forms of the strain energy density function. Wave propagation analysis will be done in section 4, based on the identified strain energy hyperelastic functions. In section we study the effect of the linear second gradient terms on the dynamical behavior of the hyperelastic medium. Dispersive analysis of the three hyperelastic models is presented in section 6. We conclude by a summary of the work and perspectives of developments in section 7.

Incremental scheme for the computation of the effective hyperelastic effective models

The adopted computational method of the effective nonlinear response of lattice materials relies on a two steps methodology: the ground state effective moduli are first evaluated in the initial small strains regime, followed by an the evaluation of the nonlinear subsequent response, based on the update of the lattice configuration (geometry) as it is subjected to an increased kinematic loading imposed over the identified unit cell.

Note that for this type of equation, harmonic plane waves cannot be considered as solution, due to the non-vanishing of the secular term.

Using the change of variable The solution of the dynamical equilibrium equation can then be expressed in terms of the displacement as

()

A A s u y f y dy sn h s dy Es Ks

in which function f(y) describes solitary waves propagation and u(y) shock waves, where s is the universal constant describing the degree of nonlinearity ( 01 s ), sn(.) the elliptic Jacobin sine, and K(s), E(s) are the complete elliptic integrals of the first and second kind respectively, 0 2 k hy  and 0 k is the propagation constant related to the wavenumber k as follows:

Wave propagation analysis for form2 with second gradient terms

The equation of motion based on the constitutive law for the second order medium (9.24) writes as follows: The problem under consideration in the PDE (9.28) represents a Boussinesq-type equation: an analogous mathematical model has been essentially studied for shallow water waves to describe the internal wave in the layered of fluid [25]. The general character of Boussinesq and Boussinesq-type equations has been intensively studied including nonlinear effects in the dynamics of solids in [34,35]; solutions of equations of this type have been studied in [36][37][38][39][40].

Using the change of variable z x t k   , where  is the frequency and k the wavenumber, and after a simple transformation, the PDE (9.28) is transformed into an ordinary differential equation (ODE in short) for the new non-dimensional strain of the wave function, function   u Nz z    :

 Computation of wave propagation characteristics in random fibrous media based on second order gradient and micropolar effective continua;  Wave propagation analysis in nonlinear elastic continua constructed by the homogenization of the response of nonlinear repetitive networks;  Efficient and general computational analysis tools for assessing nonlinear dispersion effects in repetitive network materials;  A general perturbation method accounting for geometrical nonlinearities for analyzing wave propagation at higher frequencies in network media with a nonlinear constitutive law.

The analyses performed in this work and the principal obtained results are summarized in the synthetic diagram below, relying on a classification based on two main aspects: the type of effective continuum associated to microstructural effects and the impact of the strain amplitude (small vs. large deformations).