Pierre-Yves Guerder

Samy Missoum

Krishna Muralidharan

Prof Bou Olivier

Matar

Prof Philippe Vladimir Pernod

Preobrazhensky

Keywords: phononic crystals, nonlinear elastodynamics, numerical simulations cristaux phononiques, élastodynamique non-linéaire, simulations numériques

This work is dedicated to the theoretical and numerical study of nonlinear phononic crystals. The studied nonlinearities are those due to the second (quadratic) and third (cubic) order elastic constants of the materials that constitute the crystals. Nonlinear effects are studied by the means of finite element methods, used to simulate the propagation of an elastic wave through the crystals.

A first research project concerns the study of a bone structure, namely the dispersion of elastic waves in a structure composed of collagen and hydroxy apatite alternate constituent layers. Simulations showed that it exists a strong link between bones hydration and their ability to dissipate the energy. The second study relates to an elastic resonator. A structure composed of steel inclusions in a silica matrix shows a switch behavior when the cubic nonlinearities of steel are taken into account. This strong nonlinear effect appears when the amplitude of the incident wave reaches a threshold. A full analytical model is provided. The last study demonstrates the design of composite materials with both strong cubic nonlinearities and weak quadratic nonlinearities. The derivation of the mixing laws of the elastic parameters of a nonlinear material inside a linear one is performed up to order three. Equations show a strong amplification of the nonlinear parameters of the material for some concentrations. Numerical simulations allow to conclude that the above mentioned resonator can be produced.

For this thesis, an innovative tool based on the Discontinuous Galerkin (DG) finite element method is developed for the simulation of elastic wave propagation, in linear and nonlinear systems and in finite and semi-infinite media. The implementation of this DG code for 2D and 3D simulations benefits from the efficient exploitation of modern computer infrastructure (GPU units, clusters) using the property of massive parallelization of DG algorithms.

I dedicate this thesis

to my wife Pauline.

Introduction

On appelle cristaux phononiques les structures constituées par l'ajout d'inclusions périodiques -suivant une maille cristalline -dans un matériau hôte -appelé matrice. Ces structures ont la propriété de posséder des bandes interdites, c'est-à-dire des plages de fréquences pour lesquelles l'onde acoustique est totalement réfléchie et donc interdite dans le cristal, indépendamment de la direction de propagation. Pour ces fréquences, le cristal se comporte comme un miroir acoustique parfait [1]. La fréquence centrale de la bande interdite est déterminée par la taille, la périodicité et le taux de remplissage des inclusions. Physiquement, la bande interdite est due à la diffraction des phonons aux interfaces entre la matrice et les inclusions [2]. La plupart des cristaux suivent les conditions de résonance de Bragg et Mie. La figure 1 donne un exemple de cristal phononique 2D à maille carrée et de structure de bande pouvant être associée à un tel cristal. Sur la figure, r est le rayon des inclusions, a la distance entre les centres de deux inclusions voisines. La zone de Brillouin irréductible possède trois direction ΓX, XM et MΓ. Les fréquences de résonance fondamentales de Bragg sont V /(2a) pour la direction ΓX et V /(2 √ 2a) pour la direction ΓM, où V est la vitesse acoustique moyenne dans le cristal. Cette vitesse dépend du taux de remplissage, r/a [2]. La création de défauts, par l'ajout ou le remplacement de certaines inclusions du cristal, rendant celui-ci irrégulier, autorise au contraire certaines fréquences à exister dans la bande interdite. Ces défauts agissent comme les impuretés dans les semi-conducteurs dopés. Ils permettent en partie d'accorder les cristaux [1].

L'histoire de l'étude des cristaux phononiques a commencé en 1979 avec les premières observations expérimentales par Narayanamurti et al. dans un super-réseau en GaAs/AlGaAs [3]. Les premiers travaux théoriques sont les calculs de structures de bandes pour des structures élastiques composites par Sigalas et Economou en 1992 [5]. La première bande interdite complète fut observée en 1998 par Montero de Espinosa et al. dans une plaque d'aluminium avec des arrangements périodiques de trous cylindriques remplis de mercure. Les fréquences concernées étaient de l'ordre de 1 MHz [START_REF] Montero De Espinosa | Ultrasonic band gap in a periodic twodimensional composite[END_REF]. Plus récemment, grâce à la diversité des matériaux utilisés, les bandes interdites ont atteint les « très hautes fréquences » (VHF : 30-300 MHz) et même les « ultra hautes fréquences » (UHF : 300-3000 MHz) [2]. Gorishnyy a ainsi proposé un cristal phononique hypersonique (au-dessus d'un GHz) en utilisant des trous d'air dans de l'epoxy en 2007 [START_REF] Gorishnyy | Hypersonic phononic crystals[END_REF].

Les domaines dans lesquels les cristaux phononiques ont des applications potentielles sont les communications radio-fréquences ainsi que l'imagerie ultrasons pour la médecine et le contrôle non destructif. Ces cristaux pourraient également être utilisés pour miniaturiser les lentilles acoustiques, réaliser de l'adaptation d'impédance et découpler la taille des transducteurs de leur ouverture [2]. De plus, l'imagerie avec une résolution plus fine que la limite de diffraction est envisagée, en utilisant les composantes évanescentes de l'onde. L'invisibilité acoustique pourrait être atteinte en employant les cristaux xiii comme guides d'ondes [START_REF] Lu | Phononic crystals and acoustic metamaterials[END_REF]. L'interaction phonon-photon permettrait la modulation et le refroidissement optique [2]. Les cristaux phononiques laissent également espérer une amélioration du taux de conversion directe d'énergie pour les effets thermoélectrique et thermophotovoltaïque [START_REF] Chen | Engineering nanoscale phonon and photon transport for direct energy conversion[END_REF]. Enfin, aux échelles micrométriques, leur utilisation comme filtre pour isoler des structures résonantes des vibrations extérieures et du bruit ambiant a été démontrée [2]. En effet, ceux-ci permettent de fixer ces structures au substrat de manière rigide tout en leur fournissant un environnement sans vibration, ce qui ouvre la voie à de nouveaux système mécaniques de haute précision [1].

Un des enjeux majeurs dans l'étude de tels cristaux est l'accordabilité, c'est-à-dire la possibilité de créer un cristal ayant la structure de bande voulue, statiquement -une fois pour toutes, par l'ajout de défauts -ou dynamiquement -pouvoir changer la structure de bande grâce à un paramètre physique (pression/contrainte [START_REF] Robillard | Tunable magnetoelastic phononic crystals[END_REF], température [START_REF] Bou Matar | Band gap tunability of magneto-elastic phononic crystal[END_REF], champ électrique [START_REF] Robillard | Tunable magnetoelastic phononic crystals[END_REF]...). Ces techniques ont l'inconvénient de nécessiter un contact physique avec le cristal et requièrent également un stimulus très important pour un petit résultat [START_REF] Robillard | Tunable magnetoelastic phononic crystals[END_REF]. En 2009, Robillard et al. a montré la possibilité de contrôler la structure de bandes en utilisant un champ magnétique extérieur au cristal. Cette technique utilise l'effet de magnétostriction géante et modifie les constantes élastiques non-linéaires de certains matériaux comme le Terfenol-D, de plus de 50 %, ce qui permet une accordabilité importante et sans contact [START_REF] Robillard | Tunable magnetoelastic phononic crystals[END_REF]. Bou Matar et al. a étendu ce travail en expliquant comment adapter les structures de bandes en utilisant en plus les effets de transition de réorientation de spin. L'étude modélise ainsi un cristal phononique piézomagnétique [START_REF] Bou Matar | Band gap tunability of magneto-elastic phononic crystal[END_REF].

Dans cette thèse, nous allons nous concentrer sur l'étude de différents effets non-linéaires pouvant apparaître dans ces cristaux. Nous introduirons aussi plusieurs méthodes numériques permettant de simuler efficacement la propagation d'ondes dans ces structures en 1D. Nous présenterons également une méthode innovante et hautement parallélisable, nommée méthode de Galerkin Discontinue et l'exploiterons pour l'étude de cristaux en 2D.

Étude analytique et numérique des super-réseaux 1D

Cette section donne tout d'abord une description mathématique des super-réseaux 1D puis introduit plusieurs résultats analytiques préliminaires avec les exemples associés.

Les cristaux phononiques que nous étudions ici sont des super-réseaux constitués d'une succession de couches, comme sur la figure 2. Les deux types de couches, 1 et 2, ont des épaisseurs respectives d 1 et d 2 , et des impédances respectives Z 1 et Z 2 .

...

Z 1 Z 2 Z 1 Z 2 ... ← d 1 → ← d 2 → ← d 1 → ← d 2 →
       P = Z i A(ω)e -jkx -Z i B(ω)e jkx V =
A(ω)e -jkx + B(ω)e jkx .

(

La méthode des matrices de transfert consiste à écrire la solution de l'équation (2) sous forme d'un produit de matrices. Ce formalisme nous permettra de résoudre les équations de propagation des ondes dans des structures plus complexes, comme des bicouches avec ou sans défauts. Nous écrivons ainsi :

P V x = F i H i A B 0 , avec F i = Z i -Z i 1 1
et H i = e -jk i x 0 0 e jk i x .

(

) 3
Nous en déduisons la loi de dispersion pour un bicouche infini :

cos(kd) = 1 2 Tr((F 2 H 2 F -1 2 )(F 1 H 1 F -1 1 )) = cos(k 1 d 1 ) cos(k 2 d 2 ) - 1 2
Z 1 Z 2 + Z 2 Z 1 sin(k 1 d 1 ) sin(k 2 d 2 ). (4)
Grâce à cette équation (4), nous pouvons calculer les courbes de dispersion des structures utilisées dans les chapitres 2 et 3. Ces courbes sont représentées sur la figure 3.

Nous obtenons aussi l'équation de propagation d'une onde à travers N bicouches :

t 0 = M T 1 1 r , avec M T 1 = (F -1 L F 1 )C N 1 (F -1 1 F 0 ) (5)
où t est le coefficient de transmission, r le coefficient de réflexion et : Pour N bicouches avec un défaut central :

C 1 = α β β * α * , with                  α = (cos(k 2 d 2 ) -jS sin(k 2 d 2 )) e -jk 2 d 2 β = jD 12 sin(k 2 d 2 )e jk 1 d 1 S = 1 2 Z 1 Z 2 + Z 2 Z 1 D 12 = 1 2 Z 1 Z 2 -Z 2 Z 1 . ( 6
)
t 0 = M T 2 1 r , avec M T 2 = M T L H C M T R (7)
avec :

M T L = (F -1 C F 1 )C N L 1 (F -1 1 F 0 ), M T R = (F -1 L F 2 )C N R 2 (F -1 2 F C ) et C 2 = (F -1 2 F 1 H 1 )(F -1 1 F 2 H 2 ) ( 8
)
Ces équations nous permettent de tracer et comparer les spectres en transmission du cristal phononique qui sera étudié au chapitre 3. Ces spectres apparaissent sur la figure 4. Enfin, dans le cas d'une structure composée de N bicouches quart-d'onde, nous retrouvons l'expression du coefficient de transmission en énergie (ou en puissance) donnée par Beaky et al. dans [START_REF] Beaky | Two-Dimensional Photonic Crystal Fabry-Perot Resonators with Lossy Dielectrics[END_REF] :

T = t 2 = 4n 2N (1 + n 2N ) 2 où n = - Z 1 Z 2 . ( 9
)
Ce chapitre se termine par le rappel des équations implémentées dans les méthodes numériques utilisées dans cette thèse, à savoir la méthode des éléments spectraux (SEM), la méthode des différences finies (FDTD) et la méthode de densité d'énergie spectrale (DES).

Dispersion des ondes élastiques dans une structure osseuse

Cette section présente l'étude de la dispersion des ondes dans une structure osseuse composée de couches alternées de collagène et d'hydroxy apatite. Nous relions la dissipation de l'énergie à l'hydratation des os, grâce à un modèle chimique et à des simulations numériques.

Nous commençons par établir un modèle chimique du collagène en considérant sa capacité à absorber et désorber l'eau en fonction de la contrainte (donc de la déformation) qui lui est appliquée. Les équations de la thermodynamique nous permettent d'écrire la loi de Hooke comme suit :

ε = σ Y + η(X 1 -X 0 1 ). ( 10
)
où ε est la déformation, σ la contrainte, Y le module d'Young, η le coefficient d'expansion chimique et X 1 -X 0 1 le changement de composition (X 1 = X 1 (ε) et X 0 1 = 0.5). Ainsi, nous pouvons écrire σ = Y (ε)ε où Y est une fonction non-linéaire de ε.

La structure étudiée est modélisée par un système de masses et de ressorts représenté sur la figure 5. Dans ce modèle, les 32 masses (A) représentent les couches d'hydroxy apatite et les 32 masses (B) les couches de collagène. Les ressorts ont un comportement linéaire ou non-linéaire selon que le matériau à leur gauche est (A) ou (B). La structure complète est constituée d'une chaîne de 40 cellules élémentaires (CE) à laquelle est imposée une condition de périodicité. Les constantes des matériaux sont données par le tableau 1. Comme indiqué dans ce tableau, le module d'Young Y du collagène est une fonction non-linéaire de la déformation ε. Le développement en série de Taylor de cette fonction est donnée ci-dessous. De cette expression, nous obtenons la contrainte σ = Y (ε)ε et l'énergie élastique E(ε) = ε 0 σ(ε)dε. Ainsi :

Y (ε) = a 0 + a 2 ε 2 + a 4 ε 4 + a 6 ε 6 + a 8 ε 8 . E(ε) = a 0 2 ε 2 + a 2 4 ε 4 + a 4 6 ε 6 + a 6 8 ε 6 . ( 11
)
Les coefficients de ces développements figurent dans le tableau 2. Ces formules sont valides pour |ε| 0.09 et les coefficients sont obtenus par mise en concordance (fit) et extrapolation de données expérimentales figurant dans [START_REF] Bowman | The tensile behavior of demineralized bovine cortical bone[END_REF]. Tableau 2 -Coefficients du développement en série de Taylor de Y (ε) et E(ε) pour le collagène.

Le diagramme de dispersion de la structure considérée est donné sur la figure 3a. La structure présentée est analysée par des simulations numériques utilisant la méthode FDTD-DES. Celle-ci permet d'obtenir un diagramme de dispersion constitué des spectres en transmission des différents vecteurs d'onde.

Pour de faibles amplitudes de déplacement, ce diagramme de dispersion est quasi-superposable avec celui obtenu par la méthode des matrices de transfert. Le comportement du système est donc quasi-linéaire. En revanche, pour de plus grandes amplitudes de déplacement, les pics correspondant aux modes transmis s'élargissent considérablement, leur fréquence augmente et un remplissage des bandes interdites est observé. Le comportement du système devient fortement non-linéaire.

La figure 6 montre que si nous limitons le développement en série de Taylor de l'énergie élastique à l'ordre En conclusion, nous avons montré que la non-linéarité des couches de collagène dans les os, dues à l'absorption et à la désorption d'eau, donne naissance à des processus de diffusion quatre-phonons qui conduit au remplissage des bandes interdites. Ces processus constituant des manières pour le système de dissiper l'énergie, cette dissipation est la conséquence directe de l'hydratation du collagène.

Étude d'un résonateur élastique non-linéaire

Nous démontrons dans cette section le comportement d'interrupteur d'un résonateur élastique non-linéaire constitué d'inclusions d'acier dans une matrice de silice. Ce cristal phononique est modélisé et étudiée par une méthode spectrale et son comportement expliqué par un modèle théorique complet.

Nous considérons la structure décrite par la figure 7. Il s'agit d'un cristal phononique 1D constitué d'inclusions d'acier dans une matrice de silice. Les dimensions des constituants sont données par la figure. La figure 7a représente une structure « parfaite » c'est-à-dire sans défaut : les cinq inclusions ont la même largeur, égale au quart de la longueur d'onde du signal sinusoïdal avec lequel la structure sera excitée. La figure 7b représente la même structure avec un défaut : l'inclusion centrale a une largeur double. Le diagramme de dispersion de cette structure est donné sur la figure 3b.

La méthode des matrices de transfert fournit les spectres en transmission pour la structure parfaite et la structure avec un défaut. Ces spectres sont présentés sur la figure 4. Le pic central visible sur la figure 4b correspond à la fréquence de résonance du défaut, notée f res et possède un facteur de qualité Q.

La propagation d'une onde élastique dans la structure est ensuite simulée au moyen de la méthode des éléments spectraux (SEM). Pour ce faire, une onde sinusoïdale d'amplitude variable est appliquée du côté gauche du cristal (position « source » sur la figure 7) et l'onde transmise est enregistrée sur la droite (position « récepteur »). La fréquence f c de ce signal source est choisie légèrement inférieure à f res et se situe donc dans une bande interdite. Le défaut central a un comportement non-linéaire. La relation entre la contrainte et la vitesse est donc :

∂τ ∂t = k(t) ∂v ∂z où k(t) = ∂τ ∂ε = C 11 1 -2Γε -3δε 2 . ( 12
)
Dans l'équation [START_REF] Beaky | Two-Dimensional Photonic Crystal Fabry-Perot Resonators with Lossy Dielectrics[END_REF], k(t) est défini comme dans [START_REF] Naugolnykh | Nonlinear Wave Processes in Acoustics[END_REF], v est la vitesse particulaire, τ la contrainte, C 11 est le module d'Young linéaire, ε = ∂u/∂z est la déformation et Γ et δ sont les paramètres non-linéaires quadratique et cubique, respectivement. Ces quantités sans dimension valent Γ = 0 et δ = 27.37 pour notre système.

Le défaut central constitue aussi une cavité élastique résonante qui piège et accumule l'onde incidente. Lorsque la déformation dans la cavité atteint une valeur seuil, les non-linéarités présentes provoquent un adoucissement du matériau et la fréquence de résonance du défaut, initialement f res diminue et atteint la fréquence d'excitation f c . La cavité entre alors en résonance et l'onde est parfaitement transmise à travers le cristal. Lorsque l'amplitude de l'onde incidente diminue, l'inverse se produit : la fréquence de résonance de la cavité retrouve sa valeur initiale et le signal source n'est plus transmis à travers le cristal.

Nous obtenons ainsi un système possédant un comportement d'interrupteur contrôlé par l'amplitude de la source : il est non passant pour les amplitudes de source faibles et passant pour les amplitudes de source fortes. Ce comportement hystérétique est montré sur la figure 8. Sur cette figure, les points noirs (croix et cercles) décrivent les deux états stables du système (passant en haut à droite et non passant en bas à gauche). La courbe bleue, dite en « S », correspond au modèle théorique de notre système.

En effet, si nous définissons la puissance acoustique par P i = Z i v 2 i /2, nous pouvons démontrer que les puissances d'entrée et de sortie de notre système sont liées par l'équation :

P out P in = 1 1 + P out P 0 -∆ 2 où P 0 = Z 1 ω res ω 2 4|χ|Q 2 . ( 13
)
Dans cette équation, χ est le paramètre non-linéaire, fonction de Γ et δ, ∆ = 2Q (f res -f c ) / (πf res ) est l'écart entre la fréquence d'excitation et la fréquence de résonance. P 0 est la puissance de seuil au-dessus de laquelle les effets non-linéaires rendent le cristal passant. Notre modèle théorique nous permet de prédire le comportement du système en fonction des valeurs de chaque paramètre, à savoir Q, χ et ∆. Cependant, nous avons observé que la présence de non-linéarités quadratiques (lorsque Γ > 0.01δ) avait pour conséquence de briser le processus hystérétique, sans que nous puissions, à l'heure actuelle, expliquer pourquoi.

C'est pour cette raison que les chapitres 4 et 5 portent sur l'étude des lois de mélange des paramètres élastiques quadratique et cubique. Nous allons montrer que, dans le cas d'un système constitué d'inclusions non-linéaires dans une matrice linéaires, pour certaines concentrations, les paramètres cubiques peuvent être augmentés d'un facteur important et les non-linéarités quadratiques diminuées ou peu augmentées.

Lois de mélange pour les paramètres élastiques quadratique et cubique

Nous dérivons ici jusqu'à l'ordre trois les équations des lois de mélange des constantes élastiques pour tout système composé d'inclusions non-linéaires dans une matrice isotrope, aux faibles concentrations, grâce à une procédure d'homogénéisation. L'amplification des paramètres non-linéaires est étudiée.

Nous commençons par étudier un système masse-ressort, représenté sur la figure 9 comportant N ressorts parmi lesquels N 1 ont une constante de raideur

K 1 et N 2 une constante K 2 . Nous définissons le rapport Q = K 1 /K 2 et la concentration c = N 2 /N 1 .
Nous écrivons la relation contrainte-déformation comme suit :

σ = MAε 1 + Γ 2 εB + δ 2 ε 2 C , ( 14
)
avec A le paramètre linéaire et B et C les paramètre non-linéaires quadratique et cubique, respectivement :

A = 1 1 + c(Q -1)
and

B = 1 -c + c Γ 2 Γ 1 Q 2 (1 + c(Q -1)) 2 and C = 1 -c + c δ 2 δ 1 Q 3 (1 + c(Q -1)) 3 . ( 15
)
Par dérivation de ces expressions, nous obtenons les concentrations optimales c opt qui permettent de

c opt (n = 2) = 1 Q -1 - 2 Γ 1 Γ 2 Q 2 -1 et c opt (n = 3) = 1 2(Q -1) - 3 2 Γ 1 Γ 2 Q 3 -1 soit c opt (n) ≈ 1 (n -1)Q . ( 16
)
Nous remarquons que ces concentrations optimales sont différentes pour les cas quadratique et cubique. Nous pourrons donc amplifier les paramètres cubiques sans trop changer les quadratiques.

À la concentration optimale, le paramètre linéaire A vaut 1/2 (cas quadratique) ou 2/3 (cas cubique). À la concentration optimale pour les non-linéarités quadratiques, B et C valent respectivement :

B(c opt (n = 2)) ≈ Γ 2 Γ 1 × Q 4 et C(c opt (n = 2)) ≈ δ 2 δ 1 × Q 2 8 . ( 17
)
À la concentration optimale pour les non-linéarités cubiques, B et C valent respectivement :

B(c opt (n = 3)) ≈ Γ 2 Γ 1 × 2Q 9 et C(c opt (n = 3)) ≈ δ 2 δ 1 × 4Q 2 27 . ( 18
)
Nous remarquons que le facteur d'amplification des non-linéarités cubiques est bien plus important que celui des non-linéarités quadratiques. Nous pouvons donc accroitre considérablement le rapport C/B et ainsi favoriser les non-linéarités cubiques.

Nous avons également calculé les lois de mélange des paramètres non-linéaires quadratique et cubique pour des fluides et nous obtenons les mêmes facteurs d'amplification et les mêmes concentrations optimales.

Nous allons à présent calculer les expressions des lois de mélanges des paramètres non-linéaires quadratiques et cubiques pour des inclusions d'un matériau non-linéaire 2 dans une matrice linéaire 1, comme présenté sur la figure 10 qui définit également les notations utilisées pour les constantes.

En 3D, la relation entre le tenseur des contraintes T et celui des déformations ε s'écrit ainsi : et en 2D :

T = 2µ ε + λ Tr( ε) Î + A ε2 + B Tr( ε2 ) Î + 2B ε Tr( ε) + C Tr 2 ( ε) Î + E Tr( ε3 ) Î + 3E ε2 Tr( ε) + 2F Tr( ε) Tr( ε2 ) Î + 2F ε Tr 2 ( ε) + 4G ε Tr( ε2 ) + 4H Tr 3 ( ε) Î . ( 19
)
T = 2µ ε + λ Tr( ε) Î + e Tr( ε2 ) Î + 2e ε Tr( ε) + 3f Tr 2 ( ε) Î + 4g ε Tr( ε2 ) + h Tr 3 ( ε) Î + l ε Tr 2 ( ε) + l Tr( ε) Tr( ε2 ) Î . ( 20
)
Grâce à une méthode d'homogénéisation, nous définissons un matériau équivalent à l'ensemble du système. Celui-ci possède des constantes linéaires et non-linéaires effectives dépendant des paramètres des inclusions et de la matrice ainsi que de la concentration. Nous nous intéressons seulement au cas 2D et aux constantes linéaires K eff et µ eff , à la constante quadratique f eff et à la constante cubique h eff . Ces deux constantes non-linéaires sont en effet celles qui agissent principalement sur les ondes longitudinales de compression. Nous écrivons donc la relation contrainte-déformation ainsi :

T = 2µ eff ε + (K eff -µ eff ) Tr εI + 3f eff Tr 2 εI + h eff Tr 3 εI . ( 21
)
Dans ce cas, nous pouvons définir deux formules donnant respectivement la concentration optimale et l'amplification des paramètres (f eff /K eff )/(f /K 2 ) et (h eff /K eff )/(h/K 2 ). Dans ces formules, d désigne la dimension (2 ou 3) et g l'ordre (2 ou 3) tandis que x et y désignent les rapports :

x = µ 1 /K 1 et y = K 1 /K 2 .
Ainsi, si x → 0, x > 0 et y → ∞, le facteur d'amplification non-linéaire est :

Amp(x * ) ≈ 4g 9-3g y g-1 (18 -5g) g-1 (6 -g) g+1 , ( 22
)
et la concentration optimale :

c opt (x * ) ≈ g 2y × 18 -5g 6 -g . ( 23
)
Nous exploitons maintenant ces résultats dans un exemple réel. Considérons un système composé d'inclusions d'un polymère poreux non-linéaire dans une matrice de PDMS linéaire. Le polymère poreux est constitué de 99,13 % de polymère et de 0,87 % d'air.

Matériau

ρ 0 (kg.m -3 ) v l (m/s) v t (m/s) K 2D (

Étude numérique des cristaux phononiques non-linéaires 2D

Cette section exploite les résultats de la section précédente pour montrer que l'amplification des nonlinéarités cubiques et la réduction des non-linéarités quadratiques peuvent être observées dans un système propagatif.

Nous présentons d'abord la méthode de Galerkin Discontinue et son implémentation dans le framework Hedge. La méthode de Galerkin Discontinue est une méthode d'éléments finis permettant la résolution d'équations différentielles en utilisant des polynômes d'ordre arbitraire (comme dans la méthode des éléments finis) et en utilisant des flux aux frontières entre les éléments (comme dans la méthode des volumes finis).

La particularité de cette méthode est d'autoriser des solutions discontinues aux interfaces entre les éléments, ce qui la rend hautement parallélisable, puisque la solution sur un des éléments ne dépend pas des solutions sur ses voisins. Elle est donc particulièrement adaptée à l'exécution sur les grilles de calculs (clusters) et les processeurs graphiques (GPU).

Nous avons implémenté cette méthode dans un opérateur d'élastodynamique pour le framework Hedge, en utilisant l'équation [START_REF] Boechler | Bifurcation-based acoustic switching and rectification[END_REF] qui tient compte des non-linéarités quadratiques et cubiques.

Nous avons validé le fonctionnement de cet opérateur, tant pour sa partie linéaire que pour la gestion des non-linéarités. Pour la partie linéaire, nous avons comparé le résultat fourni par cet opérateur au résultat fourni par le programme de référence EX2DDIR qui calcule analytiquement la propagation d'une onde élastique dans un demi-espace possédant une surface libre, depuis une source ponctuelle directionnelle [START_REF] Berg | Modeling the earth for oil exploration[END_REF]. Pour les non-linéarités quadratiques et cubiques, nous avons vérifié que les valeurs des harmoniques deux et trois augmentait avec la distance parcourue conformément aux formules présentes dans la littérature [START_REF] Adler | Determination of the Nonlinearity Parameter B/A for Water and m-Xylene[END_REF][START_REF] Liu | Theoretical and experimental study of the third-order nonlinearity parameter C/A for biological media[END_REF]. Comme le montrent les figures 12 (a) et (b), l'accord entre les résultats obtenus et théoriques est très bon. Dans cette structure, représentée sur la figure 13, la concentration de polymère poreux dans le PDMS est de 0,67 %, concentration censée encourager l'amplification des non-linéarités cubiques au détriment des non-linéarités quadratiques.

En enregistrant l'évolution des deuxième et troisième harmoniques en fonction de la distance parcourue par l'onde depuis le point source, nous avons pu estimer la valeur des constantes non-linéaires d'ordre deux et trois grâce aux formules de la littérature précédemment utilisées pour la validation du code.

Nous avons ainsi obtenu une amplification du paramètre quadratique par un facteur 4,7 alors que

Conclusion et perspectives

Dans cette thèse, nous avons fourni des résultats analytiques pour l'étude des cristaux phononiques linéaires (chapitre 1) puis exploré le domaine des cristaux phononiques non-linéaires et observé de nouveaux phénomènes comme la dispersion des ondes élastiques par l'interaction quatre-phonons (chapitre 2) ainsi que la commutation (switching) d'une onde élastique par un résonateur non-linéaire (chapitre 3). Ces phénomènes sont dus aux non-linéarités des matériaux constituant les cristaux.

Nous avons étudié comment ces non-linéarités évoluaient lorsque des inclusions non-linéaires sont présentes dans une matrice linéaire et avons fourni des équations permettant d'optimiser la non-linéarité quadratique ou cubique d'un tel matériau composite (chapitre 4). Nous avons enfin montré la possibilité d'utiliser ces résultats dans un système propagatif en les appliquant à une structure 2D (chapitre 5).

Pour la suite de l'étude, nous proposons d'exploiter les résultats des lois de mélange dans le résonateur 1D pour disposer d'un défaut central possédant de fortes non-linéarités cubiques et de faibles non-linéarités quadratiques. Un modèle 2D puis 3D (plaque contenant des inclusions cylindriques) de ce résonateur pourrait ensuite être proposé puis cette structure pourrait être réalisée physiquement par micro-usinage, testée et caractérisée.

Introduction

Phononic Crystals and Their Applications

Abstract Phononic crystals (PnCs) are the analogue of photonic crystals for the acoustic waves. By adding scattering periodic inclusions, following a crystalline lattice, inside an homogeneous host material, some ranges of acoustic frequencies are totally reflected and thus forbidden in the resulting material. This creates acoustic band-gaps.

Initially, the constructed crystals were limited to frequencies below 1 MHz, but the micro-crystals made by micro-fabrication can reach very high and even ultra high frequencies, with an example above 1 GHz [2].

The key issue in the field of PnCs is tuning the acoustic band-gaps in order to tailor the behavior of acoustic waves and create new functions, such as selective frequency filtering and wave guiding, and new devices, such as wavelength demultiplexers. This can be done by using the nonlinear properties of the materials that constitute the crystals as well as using the geometrical nonlinearities of the structures [1].

The initial purpose of this thesis was to demonstrate, numerically, that one can dynamically tune the band gaps and localized modes in PnCs by using magnetoelastic materials. It then evolved in studying the nonlinear effects in various phononic structures. This thesis does not only concern the analytical study of nonlinearities in PnCs but also the methods used to simulate numerically and efficiently the acoustic wave behavior in such structures. So, some facts and recent breakthroughs about the resolution of the acoustic wave equation will be explained: a state-of-the-art highly parallelizable finite element method called the Discontinuous Galerkin method will be presented.

General Overview

Principle PnCs are periodic arrangements of inclusions inside an elastic or viscoelastic material or a fluid, for example metal in air, polymer in water, air in epoxy, etc. With such a structure, band-gaps may appear and are independent of the direction of propagation of the incident elastic wave. In this case, the PnC behaves like a perfect non absorbing acoustic mirror for the rejected frequencies [1].

The central frequency of the band-gap is determined by the size, periodicity, filling and arrangement of the inclusions. Physically, the phononic band-gap is due to the diffraction of the elastic wave at the interface between the matrix and the inclusions [2].

Creating defects by replacing or removing some of the inclusions in a PnC, by making the arrangement irregular, allows certain frequencies to exist within the band-gap. Those defects behave like impurities in doped semiconductors. Thus, it is possible to tailor the acoustic properties of the crystal [1].

Typical PnCs obey the Bragg and Mie resonance conditions. Figure 14 represents a square lattice crystal. On this figure, r is the inclusion radius and a is the distance between the centers of the inclusions.

The irreducible Brillouin zone has three directions, ΓX, XM and MΓ. The fundamental Bragg resonance frequencies are V /(2a) for the ΓX direction and V /(2 √ 2a) for the ΓM direction, where V is the average acoustic velocity in the crystal. This velocity depends on the filling ratio, r/a [2].

PnCs are elastic or acoustic depending on whether the host material can (solid) or cannot (gas or liquid) support transversely polarized waves [START_REF] Lu | Phononic crystals and acoustic metamaterials[END_REF].

History

The very first known experimental observation of PnCs was in 1979, when Narayanamurti et al. investigated the propagation of high-frequency phonons through a GaAs/AlGaAs superlattice. Their superlattice can be considered as a 1D PnC [3].

As for theoretical work, Sigalas and Economou demonstrated for the first time in 1992 that frequency gaps appear for elastic waves in periodic arrangements of spheres with high density compared to that of the host material [4]. Then, Kushwaha et al. calculated full band-structure for periodic, elastic composites in 1993 [5].

However, the first complete phononic band-gap was only observed in 1998 by Montero de Espinosa et al., using an aluminum alloy plate with a square periodic arrangement of cylindrical holes filled with mercury.

A band-gap appeared in the frequency range between 1000 and 1120 kHz [START_REF] Montero De Espinosa | Ultrasonic band gap in a periodic twodimensional composite[END_REF].

Recently, the variety of materials used in the fabrication of phononic devices allowed great improvements on the reachable frequencies for the band-gap, making it increase from about 1 MHz to very high frequencies (VHF: 30-300 MHz) and even ultra high frequencies (UHF: 300-3000 MHz) [2].

Gorishnyy demonstrated an hypersonic (above 1 GHz) PnC using air scattering inclusions in epoxy in 2007. The band-gap was measured by Brillouin light scattering [START_REF] Gorishnyy | Hypersonic phononic crystals[END_REF].

Applications

The domains in which PnCs have potential applications are radio-frequency communications and ultrasound imaging for medicine and nondestructive testing. Focusing devices made with PnCs could miniaturize acoustic lenses, adapt impedance and decouple the transducer size from the aperture [2].

Moreover, sub-diffraction-limited resolution by transmitting the evanescent components of the wave and acoustic shielding could be reached by using them as wave-guides [START_REF] Lu | Phononic crystals and acoustic metamaterials[END_REF]. Using the photon-phonon interaction would allow modulation and optical cooling [2]. They could also improve direct energy conversion by thermoelectric and thermophotovoltaic effects [START_REF] Chen | Engineering nanoscale phonon and photon transport for direct energy conversion[END_REF].

At the micro-scale, PnCs are used to isolate resonating structures, such as Coriolis force gyroscopes, mechanical resonators, filters and oscillators, from external vibrations and noise [2]. Therefore, they allow the rigid attachment of these devices to the substrate in a vibration-free environment, which allows to make high-precision mechanical systems [1].

Tunable Phononic Crystals

As discussed above, adding defects in the crystalline structure is an easy way to tailor the acoustic band-gaps of the PnC. Tunability of PnCs can also be achieved by changing the geometry of the inclusions through stress or thermal effects [START_REF] Bou Matar | Band gap tunability of magneto-elastic phononic crystal[END_REF].

Another way is to vary the elastic characteristics of the constitutive materials by applying external stimuli such as electric field, temperature or stress. However those techniques either require physical contact or the application of very large stimuli for a very small result [START_REF] Robillard | Tunable magnetoelastic phononic crystals[END_REF]. In optics, Soljačić et al. published in 2002 a description of a nonlinear photonic crystal capable of performing optimal bistable switching. Their analytical model accompanied with numerical simulations describe a resonator consisting of a cavity with a nonlinear optical index. When increasing the wave amplitude, the nonlinear frequency shift brings the wave frequency to the resonance frequency of the cavity. This leads to a higher transmission of the wave through the crystal. When decreasing the wave amplitude, the transmitted power delineates an hysteresic curve and the system reaches full transmission [START_REF] Soljačić | Optimal bistable switching in nonlinear photonic crystals[END_REF]. This kind of behavior has not yet been observed in acoustics even if the approach proposed by Robillard and Bou Matar seems to prove the feasibility of an acoustic bistable switch using piezomagnetic PnCs.

Granular Phononic Crystals

The discussion above presented general facts and breakthroughs about PnCs. Despite progresses have been made to discover their properties and try to tailor their band structures, very few attempts have been done to understand their nonlinear behavior, except in the specific case of granular PnCs which will be presented now.

Nonlinear granular PnCs are composed of statically compressed chains of particles, confined in a guide, that interact nonlinearly through Hertzian contacts [START_REF] Boechler | Bifurcation-based acoustic switching and rectification[END_REF]. Uncompressed granular crystals are called "sonic vacuum type crystals" and are incapable of transmitting linear elastic waves [START_REF] Yang | Frequency-and Amplitude-Dependent Transmission of Stress Waves in Curved One-Dimensional Granular Crystals Composed of Diatomic Particles[END_REF].

History

How waves propagate in statically and dynamically loaded granular media have lead to an extensive research work since the 1980s, for example with Nesterenko, who studied the existence and interaction of solitons in packed spherical granules, using numerical methods and Shukla et al. who modeled the wave propagation in these media and computed it using experimental data [START_REF] Nesterenko | Propagation of nonlinear compression pulses in granular media[END_REF][START_REF] Shukla | Experimental and computational modeling of wave propagation in granular materials[END_REF].

In 2004, Daraio et al. designed 1D and 3D crystals made of 0.03 to 0.5 g steel spheres in a silicon or PTFE matrix and observed trains of strongly nonlinear solitary waves with small amplitudes (corresponding to forces of about 0.3 N) which propagate at 317 m/s which is below the speed of sound in air [START_REF] Daraio | Strongly nonlinear waves in 3D phononic crystals[END_REF].

Then, in 2006, the same team showed that the nonclassical, strongly nonlinear wave behavior appears in granular materials if the system is weakly compressed, which means that the precompression force is very small with respect to the wave amplitude. Oppositely, strongly compressed chain behavior approaches linear wave behavior [START_REF] Daraio | Tunability of solitary wave properties in one dimensional strongly nonlinear phononic crystals[END_REF].

The derivation of the dispersion relation of 3D granular crystals made of hexagonal arrangements of spheres was performed by Merkel et al. in 2010. They predict the existence of translational, rotational and coupled translational/rotational modes and show that the longitudinal modes are not changed when the rotational degree of freedom is taken into account [START_REF] Merkel | Dispersion of elastic waves in three-dimensional noncohesive granular phononic crystals: Properties of rotational modes[END_REF].

In 2011, Boechler et al. used granular crystals similar to the one shown on Figure 15 with a defect to design a rectifier (device that allows the propagation of some frequencies only in one direction) and an acoustic switch with sharp transitions between states [START_REF] Boechler | Bifurcation-based acoustic switching and rectification[END_REF]. Nonlinear resonances in diatomic granular chains have also been observed by Cabaret et al. who use the amplitude-dependent behavior that result from the geometry of the structures, specially the Hertzian contacts between the particles [START_REF] Cabaret | Amplitude-dependent phononic processes in a diatomic granular chain in the weakly nonlinear regime[END_REF].

Finally, in 2012, Yang and Daraio studied the propagation of stress waves in granular crystals composed of diatomic unit cells in bent elastic guides. The cells were made of centimeter-long alternating spherical and cylindrical particles. The waves were generated by striking the particle on the top with a force of 50 to 1000 N. Those structures possess band gaps around 5 and 10 kHz which can be tuned by modifying the precompression of the chain. They show highly nonlinear behaviors and propagate solitary waves [START_REF] Yang | Frequency-and Amplitude-Dependent Transmission of Stress Waves in Curved One-Dimensional Granular Crystals Composed of Diatomic Particles[END_REF].

Tunability

One of the interests of granular crystals is the possibility to tune their behavior from near-linear to highly nonlinear by varying the precompression which changes the ratio of static to dynamic particle displacement [START_REF] Boechler | Bifurcation-based acoustic switching and rectification[END_REF].

A simulation by Göncü et al. proved that the band structure of a 2D granular crystal composed of silicon rubber and polytetrafluoroethylene (PTFE) cylinders could be tuned more efficiently using pattern transformation rather than changing the particles' mechanical properties, creating new gaps around 5 kHz [START_REF] Göncü | Exploiting pattern transformation to tune phononic band gaps in a two-dimensional granular crystal[END_REF].

A way of tuning the solitary waves speed by a factor of two by applying a precompression of 2.38 N with a magnetically-induced interaction or only the gravitational preload of 0.017 N was demonstrated by Daraio et al. [START_REF] Daraio | Tunability of solitary wave properties in one dimensional strongly nonlinear phononic crystals[END_REF].

A numerical model by Spadoni and Daraio shows an acoustic lens capable of generating acoustic pulses with a pressure amplitude of 675 Pa, corresponding to 57 dB, which is two orders of magnitude larger than what is reachable with linear lenses. The focal point position is tuned by varying the time and space distribution of the static precompression of the chains that compose the lens. Experiments using an array of steel spheres in Teflon sheets produced focused waves in a polycarbonate plate in accordance with the simulations. The wavelength is determined by the size of the particles [START_REF] Spadoni | Generation and control of sound bullets with a nonlinear acoustic lens[END_REF].

Applications

Nonlinear acoustic lenses created by Spadoni and Daraio with granular PnCs can achieve better focusing and allow higher focal power than linear ones. So, they are expected to improve the performances of current devices for biomedical imaging, nondestructive testing and sonars. They could also be used to constitute nonintrusive scalpels, for cancer treatment, for example [START_REF] Spadoni | Generation and control of sound bullets with a nonlinear acoustic lens[END_REF].

The ability of nonlinear granular PnCs to behave as rectifiers or switches suggest that they could be used to control the flow of energy in several applications, such as energy-harvesting materials with frequencydependent absorption and emission and thermal computers. Their nonlinear response allows them to change their state when small perturbations are applied, which makes them suited for sensing applications [START_REF] Boechler | Bifurcation-based acoustic switching and rectification[END_REF].

The control of nonlinear resonances in diatomic granular chains are expected to lead to the creation of new devices, such as passive amplitude-dependent amplifiers and attenuators [START_REF] Cabaret | Amplitude-dependent phononic processes in a diatomic granular chain in the weakly nonlinear regime[END_REF].

One of the possible future direction in research about those crystal could be the ability to engineer the dispersion relation and create gaps to create tunable vibration filtering devices and systems that could be insulated from noise vibrations [START_REF] Deymier | Acoustic Metamaterials and Phononic Crystals[END_REF].

Nanoscale Phononic Crystals

The first part provided a general description of PnCs. The second part explained the specific case of granular PnCs, whose nonlinear properties have been investigated. Another particular domain, for which a brief review is now provided, consists of nanoscale PnCs. Those crystals deal with very high frequencies (THz) and thermal effects. In this field, recent efforts have been made to engineer the band gaps. The most considered materials for those crystals are graphene and nanoporous silicon.

History

In 2009, Gillet et al. demonstrated that 3D arrays of germanium quantum dots in silicon reduce the thermal conductivity by several orders of magnitude compared to bulk silicon. This effect is produced by THz phonon confinement and reflection on layer interfaces which decrease the phonon group velocity [START_REF] Gillet | Atomic-Scale Three-Dimensional Phononic Crystals With a Very Low Thermal Conductivity to Design Crystalline Thermoelectric Devices[END_REF].

A similar effect was observed by Marconnet et al. in silicon nanobridges where the thermal conductivity can be reduced to 3 % of its value in bulk material, without affecting the electrical conductivity. They measured a thermal conductivity one order of magnitude lower than predicted by the model [START_REF] Marconnet | Phonon Conduction in Periodically Porous Silicon Nanobridges[END_REF].

Then, this effect was explained in 2012 by Dechaumphai and Chen who modeled the coherent (wave-like) and incoherent (particle-like) behavior of the phonons and computed the dispersion relation with a FDTD method. They found that the zone folding effect (formation of mini-bands in the superlattice because the Brillouin zone edge for two materials is smaller than the Brillouin zones of each material) has a major impact on the thermal conductivity [START_REF] Dechaumphai | Thermal transport in phononic crystals: The role of zone folding effect[END_REF].

A simulation by Sgouros et al. models PnCs made of graphene whose defects are constituted by adding carbon atoms, removing some or replacing some by silicon atoms. This study shows that only the substitution by silicon atoms creates band gaps [START_REF] Sgouros | Phononic band gap engineering in graphene[END_REF].

As explained by Maldovan in 2013, the control of thermal conduction can be achieved by transforming the heat flow to wave phonon transport by applying a 2D holes pattern to nanostructured alloys of silicon. By blocking high-frequency phonons (above 1 THz), preventing them from existing, the heat transfer can only be done by low-frequency phonons (below 1 THz) which can be guided by nanoscale crystals [START_REF] Maldovan | Narrow Low-Frequency Spectrum and Heat Management by Thermocrystals[END_REF].

Applications

The possible applications of nanoscale PnCs considered by Gillet et al. concern the design of new energyconversion devices, such as Peltier refrigerators for domestic use as well as highly efficient hybrid thermal/electric car engines. They also note that the thermoelectric devices should be CMOS-compatible (based on silicon, for instance) in order to improve the power management of computer microprocessors with high-density of transistors [START_REF] Gillet | Atomic-Scale Three-Dimensional Phononic Crystals With a Very Low Thermal Conductivity to Design Crystalline Thermoelectric Devices[END_REF].

The crystals designed by Maldovan open a new field called "thermocrystals" whose expected applications are heat waveguides, thermal lattices, heat imaging, thermo-optics, thermal diodes and thermal cloaking [START_REF] Maldovan | Narrow Low-Frequency Spectrum and Heat Management by Thermocrystals[END_REF].

State of the Art Summary

Based on the state of the art provided above, the following observations can be made.

• On one hand, nonlinearities have only been studied in granular chains, which use the geometrical nonlinearities not the nonlinear elastic constants of the materials.

• On the other hand, there is still no simulation effectively showing nonlinear behavior in a nongranular (bulk) PnC. The need for a more general understanding of nonlinear effects in PnCs have been stated by different authors [START_REF] Cabaret | Amplitude-dependent phononic processes in a diatomic granular chain in the weakly nonlinear regime[END_REF]. This is why this thesis will concentrate on exploring the nonlinear effects in PnCs.

This introduction chapter presented both the history and the state of the art concerning the study of PnCs and provided a literature review.

Chapter 1 will deal with a mathematical description of one-dimensional superlattices. It will introduce some preliminary results and equations that are useful for the modeling of 1D PnCs. It will also present the main numerical methods commonly used for this purpose.

Chapter 2 will present the study of wave dispersion in a bone structure composed of collagen and hydroxy apatite alternate constituent layers. Energy dissipation in bones will be related to their hydration through a chemical model and numerical simulations.

Introduction

Chapter 3 demonstrates the switching behavior of a nonlinear elastic resonator constituted of steel inclusions in a silica matrix. A PnC is modeled and studied with the spectral elements method. A theoretical model is provided to explain its hysteretic behavior.

Chapter 4 provides a full derivation of the mixing laws of second and third order elastic constants for any system composed of nonlinear inclusions in a linear isotropic matrix, for small concentrations, through a homogenization procedure. The amplification of the nonlinear parameters is studied.

Chapter 5 exploits the results of chapter 4 and shows that amplification of cubic nonlinearities and reduction of quadratic nonlinearities can be obtained in a heterogeneous propagative system.

Chapter1

Analytical and Numerical Study of 1D Superlattices

Introduction

This chapter provides some essential results for the study of one-dimensional phononic crystals. The first section will define 1D superlattices and their governing equations.

The second part will detail an analytical method to solve the elastic wave propagation in those structures, called the matrix transfer method. It will then explain how we can calculate analytically the dispersion curves -or band diagrams -of structures by resolving the wave equations for the frequency or wave vector.

In 1D, the wave equation has analytical solutions that will be computed using a transfer matrix method.

We will study the wave propagation in and transmission through bilayers as well as the analytical expressions of the band structure.

In a third section, we will present the numerical methods used for the study of 1D PnCs in this thesis, namely the Pseudospectral method, the Finite Difference Time Domain method and the Spectral Energy Density method.

One-Dimensional Superlattices

The structures studied in this chapter are 1D periodic crystals consisting in a succession of layers, as shown on Figure 1.1. There are two types of layers, 1 and 2, whose thickness are respectively d 1 and d 2 , and whose impedances are respectively Z 1 and Z 2 .

...

Z 1 Z 2 Z 1 Z 2 ... ← d 1 → ← d 2 → ← d 1 → ← d 2 →

9

We first write the fundamental equations of acoustics:

           ∂p ∂x = -ρ(x) ∂v ∂t ∂p ∂t = -Y (x) ∂v ∂x , c = Y ρ , (1.1)
where p is the pressure, v is the particle velocity, ρ is the mass density and Y is the Young's modulus.

Those equations are equivalent to:

             ∂τ ∂x = ρ(x) ∂ 2 u ∂t 2 τ = Y (x) ∂u ∂x , ( 1.2)
where τ = -p is the stress and u is the particle displacement.

For an harmonic wave, the generic solution in a layer with impedance Z i is an harmonic solution, i.e., a plane wave composed of an incident wave with amplitude A and a reflected wave with amplitude B. In this case, considering p = Pe jωt and v = V e jωt , we can write the solution as follows:

       P = Z i A(ω)e -jkx -Z i B(ω)e jkx V = A(ω)e -jkx + B(ω)e jkx .
(1.3)

Analytical Method: Transfer Matrix Method

The Transfer Matrix Method will use the solution described by equation (1.3) written as a product of a matrix and a vector. This formalism will allow us to write the solutions of wave propagation in more complex structures, such as sets of bilayers with or without defects. We will show that particular cases, such as quarter-wavelength layers have a more simple and interesting behavior.

Propagation in a Layer

From equation (1.3), the propagation in a layer, from position x = 0 to position x can be written as a product of matrices:

      P V       x =       Z i e -jk i x -Z i e jk i x e -jk i x e jk i x             A B       0 =       Z i -Z i 1 1             e -jk i x 0 0 e jk i x             A B       0 , ( 1.4)
which we can write:

      P V       x = F i H i       A B       0 , with F i =       Z i -Z i 1 1       and H i =       e -jk i x 0 0 e jk i x       . (1.5)
In equation (1.5), H i describes the wave propagation in the layer i at a distance x.

We can also note that:

      P V       0 = F i       A B       0 , (1.6) so:       A B       0 = F -1 i       P V       0 , with F -1 i = 1 2       1 Z i 1 -1 Z i 1       . (1.7)
This way, we obtain for the propagation through a layer of impedance Z i and a length x:

      P V       x = F i H i F -1 i       P V       0 = M i       P V       0 , with M i = F i H i F -1 i =       cos(k i x) -jZ i sin(k i x) - j Z i sin(k i x) cos(k i x)       .
(1.8)

Propagation in a Bilayer

To write the propagation in a bilayer, we use the continuity property of P and V at the interfaces:

      P V       x + =       P V       x - , ( 1.9)
where x + belongs to Z 2 and x -belongs to Z 1 .

Thus, we have:

      P V       d 1 +d 2 = M 2 M 1       P V       0 = M       P V       0 , ( 1.10)
with:

M =       cos(k 2 d 2 ) -jZ 2 sin(k 2 d 2 ) - j Z 2 sin(k 2 d 2 ) cos(k 2 d 2 )             cos(k 1 d 1 ) -jZ 1 sin(k 1 d 1 ) - j Z 1 sin(k 1 d 1 ) cos(k 1 d 1 )       =        cos(k 1 d 1 ) cos(k 2 d 2 ) -Z 2 Z 1 sin(k 1 d 1 ) sin(k 2 d 2 ) -jZ 1 cos(k 2 d 2 ) sin(k 1 d 1 ) -jZ 2 sin(k 2 d 2 ) cos(k 1 d 1 ) - j Z 2 cos(k 1 d 1 ) sin(k 2 d 2 ) - j Z 1 cos(k 2 d 2 ) sin(k 1 d 1 ) -Z 1 Z 2 sin(k 1 d 1 ) sin(k 2 d 2 ) + cos(k 1 d 1 ) cos(k 2 d 2 )        ,
or in a condensed form:

M = (F 2 H 2 F -1 2 )(F 1 H 1 F -1 1 ). (1.11)
In a periodic medium (d = d 1 + d 2 ), the propagating modes are the Bloch modes and the eigen values of M are given by: Y 1 = e jkd and Y 2 = e -jkd .

(1.12)

Since the trace of the matrix M is:

Tr(M) = Y 1 + Y 2 = 2 cos(kd), (1.13)
the dispersion law of the bilayer periodic infinite medium can be written as:

cos(kd) = 1 2 Tr(M) = cos(k 1 d 1 ) cos(k 2 d 2 ) - 1 2 Z 1 Z 2 + Z 2 Z 1 sin(k 1 d 1 ) sin(k 2 d 2 ). (1.14)
Using equation (1.14), we can calculate the dispersion curves of layered media such as the bone structure of chapter 2 and the periodic inclusions of steel in silica of chapter 3. Figure 1.2 shows those dispersion curves for such infinite media. On Figure 1.2b, the red dots represent the analytical solutions for the band edges, when a = 1 (see section 1.3.3); as Z 1 Z 2 , the dots don't exactly match the edges. For a 1D superlattice, the first Brillouin zone extends for wave vectors k between -π/d and +π/d.

Considering the term F -1 2 F 1 appearing in equation (1.11), we can now give a physical interpretation of the matrices F i :

F -1 2 F 1 = 1 2        1 Z 2 1 -1 Z 2 1              Z 1 -Z 1 1 1       = Z 1 + Z 2 2Z 2        1 Z 2 -Z 1 Z 1 +Z 2 Z 2 -Z 1 Z 1 +Z 2 1        = 1 t p       1 r p r p 1       = M t 12 , (1.15)
where .16) is the reflection coefficient between a medium of impedance Z 1 and one of impedance Z 2 and similarly

r p = Z 2 -Z 1 Z 1 + Z 2 , ( 1
t p = 2Z 2 Z 1 + Z 2 , (1.17)
is the transmission coefficient of the same media, so it appears that M t through an interface between medium 1 and medium 2.

Band Structure in a 1D Phononic Crystal

We now search for analytical solutions of the dispersion equation (1.14) rewritten as: cos(kd) = cos(X) cos(aX) -β sin(X) sin(aX) = f (X), (1.18) with:

X = k 1 d 1 , aX = k 2 d 2 and β = 1 2 Z 1 Z 2 + Z 2 Z 1 > 1. (1.19)
First, we are looking for analytical expressions of the position and width of the band gaps. At the band edge, kd = nπ so:

cos(kd) = f (X) = ±1.
(1.20)

We will consider separately the two edges, where kd = 2nπ, i.e. where f (X) = 1 and where kd = (2n + 1)π, i.e. where f (X) = -1.

First case: f (X) = -1

In this case, the dispersion equation (1.18) becomes: cos(X) cos(aX) -β sin(X) sin(aX) = -1,

which we can write:

1 -t 2 1 + t 2 × 1 -τ 2 1 + τ 2 -β × 2t 1 + t 2 × 2τ 1 + τ 2 = -1, (1.22)
where we have defined:

t = tan X 2 and τ = tan aX 2 . (1.23)
This equation is equivalent to:

1 + t 2 τ 2 -2βtτ = 0. (1.24)
Defining Y = tτ gives:

Y 2 -2βY + 1 = 0. (1.25)
The two solutions of equation (1.25) are:

Y = β ± β 2 -1, ( 1
α = 1 2 Z 1 Z 2 + Z 2 Z 1 + 1 4 Z 1 Z 2 + Z 2 Z 1 2 -1, (1.28)
or after a simplification to:

α = 1 2Z 1 Z 2 Z 2 1 + Z 2 2 ± Z 2 1 -Z 2 2 .
(1.29)

Finally, the two solutions of equation (1.25) are given by:

Y = Z 1 Z 2 and Y = Z 2 Z 1 .
(1.30)

Second case: f (X) = 1
In this case, the dispersion equation is:

cos(X) cos(aX) -β sin(X) sin(aX) = 1, (1.31)
which we can write as:

1 -t 2 1 + t 2 × 1 -τ 2 1 + τ 2 -β × 2t 1 + t 2 × 2τ 1 + τ 2 = 1.
(

This equation is equivalent to:

t 2 + τ 2 + 2βtτ = 0, (1.33)
or:

t 2 τ 2 + 2β t τ + 1 = 0 if τ 0. (1.34)
Defining Z = t τ gives:

Z 2 + 2βZ + 1 = 0. (1.35)
The two solutions of equation (1.35) are:

Z = -β ± β 2 -1, (1.36)
or introducing α as in the previous case:

α and -

1 α with α = β + β 2 -1, (1.37)
the two solutions of equation (1.35) are then:

Z = - Z 1 Z 2 and Z = - Z 2 Z 1 .
(1.38)

Remark

As we have defined:

X = k 1 d 1 and aX = k 2 d 2 , (1.39) then if d 1 = d 2 ,
we obtain:

a = k 2 d 2 k 1 d 1 = k 2 k 1 . (1.40)
So, if we impose ρ 1 = ρ 2 then:

Z 2 Z 1 = ρ 2 c 2 ρ 1 c 1 = c 2 c 1 = c 2 ω × ω c 1 = k 1 k 2 = 1 a .
(1.41)

Special case: a = 1
In this case, we get from the dispersion equation, for f (X) = 1:

cos 2 (X) -β sin 2 (X) = 1, (1.42)
which gives:

-(4 + 4β)t 2 = 0. (1.43)
The solution of equation (1.43) is:

t = 0, so X = 2nπ. (1.44)
Now, when f (X) = -1, we have:

cos 2 (X) -β sin 2 (X) = -1, (1.45)
which we can write as:

1 -t 2 1 + t 2 2 -β 2t 1 + t 2 2 = -1.
(1.46)

Equation (1.46) can be simplified, leading to:

t 4 -2βt 2 + 1 = 0. (1.47)
The four solutions of equation (1.47) are represented on Figure 1.3. They are given by:

tan X 2 ∈ Z 1 Z 2 , - Z 1 Z 2 , Z 2 Z 1 , - Z 2 Z 1 . (1.48)
We only consider solutions where:

X = ω c 1 d 1 0, (1.49)
because ω 0. The values that we obtain for X/2 between 0 and π/2 are met periodically every X/2 = nπ, or Defining

X = 2nπ.
θ 1 = arctan Z 2 Z 1 π 4 if Z 2 < Z 1 and θ 2 = π 2 -θ 1 , (1.50)
we get:

tan(θ 1 ) = sin(θ 1 ) cos(θ 1 ) = Z 2 Z 1 . (1.51)
Calculating now tan(θ 2 ) we get: .52) which means that the first gap extends from 2θ 1 to π -2θ 1 giving, in terms of frequency:

tan(θ 2 ) = sin π 2 -θ 1 cos π 2 -θ 1 = cos(θ 1 ) sin(θ 1 ) = Z 1 Z 2 , ( 1
f ∈ θ 1 c 1 πd 1 , c 1 2d 1 - θ 1 c 1 πd 1 . (1.53)
It appears that the middle of the first gap is met at:

f c = c 1 4d 1 , (1.54)
and the width of the first gap is:

∆f = c 1 2d 1 - 2θ 1 c 1 πd 1 = c 1 2d 1 1 - 4θ 1 π = 2f c 1 - 4θ 1 π . (1.55)
From equation (1.54) and (1.55), the relative width of the first gap can be expressed as:

∆f f c = 2 1 - 4θ 1 π . (1.56)
Following similar steps, we calculate the edges of the second gap. We start with:

tan(θ 3 ) = - Z 1 Z 2 = -tan(θ 2 ) = tan(-θ 2 ), (1.57)
which leads to:

θ 3 = π -θ 2 = π - π 2 -θ 1 = π 2 + θ 1 . (1.58)
The second edge of the second gap is given by θ 4 = πθ 1 . The second gap extends from 2θ 3 to 2θ 4 , so for:

f ∈ c 1 2d 1 + θ 1 c 1 πd 1 , c 1 d 1 - θ 1 c 1 πd 1 . (1.59)
The middle of the second gap is met at:

3f c = 3c 1 4d 1 , (1.60)
and the width of the second gap is given by:

∆f = c 1 2d 1 - 2θ 1 c 1 πd 1 = c 1 2d 1 1 - 4θ 1 π = 2f c 1 - 4θ 1 π , (1.61)
which is the same width as the first gap.

To sum up, we have shown that the central frequency and width of the gap, when a = 1, i.e.

k 1 d 1 = k 2 d 2 are: f cn = (2n -1) c 1 d 1 , for n 1 and ∆f = c 1 2d 1 1 - 4θ 1 π with θ 1 = arctan Z 2 Z 1 , Z 2 < Z 1 . (1.62)
In fact, the central frequencies of the gaps correspond to k 1 d 1 = π/2 so d 1 = (2n -1)λ 1 /4 (and as a = 1,

d 2 = (2n -1)λ 2 /4).
Moreover, the widths of the gap increase as θ 1 decreases, so when Z 1 is very large compared to Z 2 .

Finally, with a = 1, if a frequency f 0 is in the first gap, then, the higher odd harmonics 3f 0 , 5f 0 ... will also be in band gaps, whereas the even harmonics 2f 0 , 4f 0 ... will be between gaps, i.e. in passing bands.

Special case: a = 2

In this case, we have:

k 2 d 2 = 2k 1 d 1 and τ = tan(X) = 2t 1 -t 2 .
(1.63)

Introducing as previously

Y = tτ = 2t 2 1 -t 2 , (1.64)
we obtain from the dispersion equation, in the case where f (X) = 1:

Y = t 2 (2 + Y ) and t = ± Y 2 + Y . (1.65)
Equation ( 1.65) has the following four solutions:

t = tan X 2 ∈ Z 1 2Z 2 + Z 1 , - Z 1 2Z 2 + Z 1 , Z 2 2Z 1 + Z 2 , - Z 2 2Z 1 + Z 2 . (1.66)
Now, in the case where f (X) = -1, introducing

Z = t τ = t(1 -t 2 ) 2t = 1 -t 2 2 , ( 1.67)
we derive the following equation:

t 2 = 1 -2Z = 1 + 2Y or simply t = ± √ 1 + 2Y . (1.68)
Equation (1.68) also has four solutions:

t = tan X 2 ∈ 2Z 1 + Z 2 Z 2 , - 2Z 1 + Z 2 Z 2 , 2Z 2 + Z 1 Z 1 , - 2Z 2 + Z 1 Z 1 . (1.69)
We finally have eight solutions, represented on Figure 1.4, given by:

t = tan X 2 ∈ δ 1 , δ 2 , -δ 1 , -δ 2 , 1 δ 1 , 1 δ 2 , - 1 δ 1 , - 1 δ 2 , , (1.70)
where

δ 1 = 1 tan θ 1 and δ 2 = 1 tan θ 2 . (1.71)
If we define:

θ 1 = arctan 1 δ 1 and θ 2 = arctan 1 δ 2 , ( 1.72)
then:

θ 3 = π 2 -θ 2 and θ 4 = π 2 -θ 1 .
(1.73) In this case, in the frequency band 0, c 1 2d 1

, we have two gaps:

Gap1: f ∈ θ 1 c 1 πd 1 , θ 2 c 1 πd 1 k = π a , Gap2: f ∈ c 1 2d 1 -θ 2 c 1 πd 1 , c 1 2d 1 -θ 1 c 1 πd 1 , ( 1.74)
where the middle frequencies are given by:

Gap1: f c1 = θ 1 +θ 2 2 × c 1 πd 1 , Gap2: f c2 = c 1 2d 1 -θ 1 +θ 2 2 × c 1 πd 1 , ( 1.75)
and the gap widths by: Gap1:

∆f 1 = c 1 (θ 2 -θ 1 ) πd 1 , Gap2: ∆f 2 = ∆f 1 .
(1.76)

This resolution of the dispersion equation (1.18), to obtain analytical expressions of the central frequencies and widths of gaps, allows us to design crystals with this specific feature: if a frequency f 0 is in the first gap, then all its higher harmonics 2f 0 , 3f 0 , 4f 0 ... will also be in gaps. However, as we loose the symmetry that we had with a = 1, those harmonics will not necessarily be in the middle of the gaps and, consequently, their attenuation may not be optimal.

Propagation of Amplitudes

We can also link the amplitudes A and B at position d to amplitudes A and B at position 0 with:

      P V       d - 1 = F 1 H 1       A B       0 and       P V       d + 1 =       P V       d - 1 , (1.77) leading to:       P V       d + 1 = F 1 H 1       A B       0 and       A B       d + 1 = F -1 2       P V       d + 1 . (1.78)
After propagation across the Z 2 layer, we obtain:

      P V       d + =       P V       d - = (F 2 H 2 )(F -1 2 F 1 H 1 )       A B       0 , ( 1.79)
and finally:

      A B       d + = (F -1 1 F 2 H 2 )(F -1 2 F 1 H 1 )       A B       0 = C 1       A B       0 , (1.80)
where we have introduced the following notations:

C 1 =       α β β * α *       , with                      α = (cos(k 2 d 2 ) -jS sin(k 2 d 2 )) e -jk 2 d 2 β = jD 12 sin(k 2 d 2 )e jk 1 d 1 S = 1 2 Z 1 Z 2 + Z 2 Z 1 D 12 = 1 2 Z 1 Z 2 -Z 2 Z 1 . (1.81)
So, we have expressed the values of the amplitudes A and B at position d as a product between a matrix and the vector of the initial values of those amplitudes.

We note that:

det(C 1 ) = αα * -ββ * = 1, (1.82)
and:

1 2 Tr(C 1 ) = α + α * 2 = Re(α) = cos(k 1 d 1 ) cos(k 2 d 2 ) -S sin(k 1 d 1 ) sin(k 2 d 2 ) = 1 2
Tr(M).

(1.83)

Transmission Through a Bilayer

We now study the propagation of an elastic wave through the bilayer described on Figure 1.5 with the formalism of the preceding paragraphs. We want here to calculate the coefficients of transmission and reflection through such a bilayer structure.

Z 0 Z 1 Z 2 Z L 1 → → t r ← ← d 1 → ← d 2 → Figure 1
.5 -A 1D bilayer with impedances Z 1 and Z 2 with external layers with impedances Z 0 and Z L .

For arbitrary Z 0 and Z L , we derive using equation (1.80) the following set of equations:

      A B       L = (F -1 L F 1 )C 1 (F -1 1 F 0 )       A B       0 = M       A B       0 . (1.84)
Now, using the fact that no wave enters the bilayer from the output medium of impedance Z L , i.e. B L = 0, we obtain the coefficients of transmission t and reflection r which are given by t = A L /A 0 and r = B 0 /A 0 :

      t 0       =       M 11 M 12 M 21 M 22             1 r       0 .
(1.85)

Resolving the system of equations, we obtain:

r = - M 21 M 22 , and t = det(M) M 22 .
(1.86)

Transmission Through N Bilayers

We now study the propagation of a elastic wave through a system composed of N bilayers, as shown on

Z 0 Z 1 Z 2 Z 1 Z 2 Z 1 Z 2 Z L 1 → → t r ← 1 2 ... N ← d 1 → ← d 2 → Figure 1.6 -A 1D system composed of N bilayers.
According to the previous paragraphs, for a system of N bilayers, we can write:

      A B       LN = C 1       A B       L(N -1) = C N 1       A B       0 .
(1.87) So, considering an input layer with impedance Z 0 and an output layer with impedance Z L , the transmission is given by:

      t 0       = M T       1 r       , with M T = (F -1 L F 1 )C N 1 (F -1 1 F 0 ). (1.88)
Figure 1.7 represents the normalized transmission spectrum for 5 and 7 steels inclusions in silica, respectively. We can observe that the number of oscillations on the left and right sides of the gap is equal to the number of inclusions. We can also note that the attenuation in the gap increases with the number of inclusions.

Transmission Through N Bilayers with a Defect

We now study the propagation of a elastic wave through the following system composed of N bilayers and a central defect, as represented on Figure 1.8: In this case, the transmission is given by:

Z 0 Z 1 Z 2 Z 1 Z 2 Z C Z 2 Z 1 Z 2 Z 1 Z L 1 → → t r ← 1 ... N L 1 ... N R ← d 1 → ← d 2 → ← d C →
      t 0       = M T       1 r       , with M T = M T R H C M T L , ( 1.89)
with:

M T L = (F -1 C F 1 )C N L 1 (F -1 1 F 0 ), M T R = (F -1 L F 2 )C N R 2 (F -1 2 F C ), C 2 = (F -1 2 F 1 H 1 )(F -1 1 F 2 H 2 ).
(1.90) Figure 1.9 represents the normalized transmission spectrum for 5 and 7 steels inclusions in silica, respectively, with a central defect (d C = 2d 1 ). We can observe a peak in the band gap. Due to the symmetry of the system, the peak is exactly at the middle of the gap. Also, its width decreases with the number of inclusions (its quality factor increases). The central position of this peak is optimal as the attenuation is higher at the center of the band gap.

Reflected Impedance Through a Multilayer

We will now calculate equivalent impedance of impedance Z L seen through a system composed of N bilayers, as represented on Figure 1.10. The impedance at the relative position x in a layer with impedance Z C is defined by:

Z C Z 1 Z 2 Z 1 Z 2 Z 1 Z 2 Z L 1 2 ... N ← L → ← d 1 → ← d 2 →
Z(x) = P V x = Z C Ae -jk c x -Be jk c x Ae -jk c x + Be jk c x . (1.91)
In the semi-infinite medium with impedance Z L , there is no reflected wave, so:

∀x L + N d, Z(x) = Z L . (1.92)
Using the transfer matrices described in section 1.3.6, one can write:

      A N B N       = C N -1 1       A 1 B 1       , ( 1.93)
and:

      A Z L 0       = (F -1 Z L F 2 H 2 )(F -1 2 F 1 H 1 )       A N B N       = (F -1 Z L F 1 )(F -1 1 F 2 H 2 )(F -1 2 F 1 H 1 )       A N B N       . (1.94)
Introducing equation (1.93) into equation (1.94), we obtain:

      A Z L 0       = (F -1 Z L F 1 )C 1       A N B N       = N       A 1 B 1       , with N = (F -1 Z L F 1 )C N 1 . (1.95)
Inverting this equation, we get:

      A 1 B 1       = N -1       A Z L 0       , with N -1 = 1 N 22 N 11 -N 12 N 21       N 22 -N 12 -N 21 -N 11       , ( 1.96)
where N 11 , N 22 , N 12 and N 21 are the elements of the matrix N .

Using the definition of the acoustic impedance in the layer of impedance Z 1 at position L:

Z(L) = Z(x = L + ) = Z 1 A 1 -B 1 A 1 + B 1 = Z 1 (N -1 ) 11 -(N -1 ) 21 (N -1 ) 11 + (N -1 ) 21 , (1.97)
we finally obtain:

Z(L) = Z 1 N 22 + N 21 N 22 -N 21 . (1.98)
We note that Z(L) can be complex, thus introducing a phase shift when the wave is reflected on the interface between the cavity and the supercell.

Transmission Through N Bilayers with Quarter-Wavelength Layers

We want to obtain a simple analytical expression of the transmission coefficient of a wave through N bilayers with layer thickness equal to the quarter of the wavelength. This will enable us to find the link between the quality factor Q of the resonator described in section 1.3.7, i.e. a superlattice with a defect, and the number of bilayers and the impedance mismatch between the layers of the structure. We first present two preliminary results before performing the calculation.

Preliminary results

It is easy to check the following diagonalization:

      α β β α       =       1 1 1 -1             α + β 0 0 α -β             1 2 1 2 1 2 -1 2       , (1.99)
and the following commutation:

      α 1 β 1 β 1 α 1             α 2 β 2 β 2 α 2       =       α 2 β 2 β 2 α 2             α 1 β 1 β 1 α 1       =       α 1 α 2 + β 1 β 2 α 1 β 2 + α 2 β 1 α 2 β 1 + α 1 β 2 α 1 α 2 + β 1 β 2       . (1.100)

Calculation

Since we use quarter-wavelength layers, we rewrite the matrix C 1 of equation (1.81) using:

d 1 = λ 1 4 and d 2 = λ 2 4
, or equivalently:

k 1 d 1 = k 2 d 2 = π 2 . (1.101)
So, we obtain:

C 1 =       α β β α       , with                      α = -S β = -D 12 S = 1 2 Z 1 Z 2 + Z 2 Z 1 D 12 = 1 2 Z 1 Z 2 -Z 2 Z 1 . (1.102)
We note that in this case:

det(C 1 ) = α 2 -β 2 = 1, (1.103)
and:

α + β = Z 1 Z 2 , α -β = Z 2 Z 1 . (1.104)
Using a recurrence and the diagonalization given in equation (1.99), we can easily prove that:

C N 1 =       α N β N β N α N       , with          α N = (α+β) N +(α-β) N 2 β N = (α+β) N -(α-β) N 2
.

(1.105)

Moreover, we note that:

α 2 N -β 2 N = det C N 1 = (det (C 1 )) N = 1, (1.106)
and:

α N + β N = (α + β) N , α N -β N = (α -β) N . (1.107)
In the general case, with no hypothesis on the values of Z 0 and Z L , according to equation (1.88) the coefficients of reflection and transmission are given by:

r = - M T 21 M T

22

, and t = det(M T ) M T

22

.

(1.108)

Using expressions (1.5) and (1.7), we calculate:

F -1 1 F 0 = 1 2        1 Z 1 1 -1 Z 1 1              Z 0 -Z 0 1 1       = 1 2        1 + Z 0 Z 1 1 -Z 0 Z 1 1 -Z 0 Z 1 1 + Z 0 Z 1        .
(1.109)

So, given their shapes, according to equation (1.100), C N 1 and (F -1 1 F 0 ) commute. So:

M T = (F -1 L F 1 )C N 1 (F -1 1 F 0 ) = (F -1 L F 1 )(F -1 1 F 0 )C N 1 = (F -1 L F 0 )C N 1 , (1.110)
and we obtain:

F -1 L F 0 = 1 2       1 Z L 1 -1 Z L 1             Z 0 -Z 0 1 1       = 1 2       1 + Z 0 Z L 1 -Z 0 Z L 1 -Z 0 Z L 1 + Z 0 Z L       , (1.111)
leading to the following expression of the transmission matrix M T :

M T = 1 2       1 + Z 0 Z L 1 -Z 0 Z L 1 -Z 0 Z L 1 + Z 0 Z L             α N β N β N α N       = 1 2       α N + β N + (α N -β N ) Z 0 Z L α N + β N + (β N -α N ) Z 0 Z L α N + β N + (β N -α N ) Z 0 Z L α N + β N + (α N -β N ) Z 0 Z L       . (1.112)
Introducing theses expressions in equation (1.108), the transmission coefficient becomes:

t = det M T M T 22
where det

M T = det F -1 L det F 0 det C N 1 = Z 0 Z L and M T 22 = 1 2 (n) N + 1 n N Z 0 Z L , ( 1.113)
where we have introduced the following notations:

α + β = Z 1 Z 2 = n and α -β = Z 2 Z 1 = 1 n . (1.114)
These expressions lead to:

t = Z 0 Z L 2 (n) N + 1 n N Z 0 Z L = 2Z 0 Z L n N + Z 0 n N = 2n N Z L Z 0 n 2N + 1 . (1.115)
Then, the coefficient of transmission in energy can be written as:

T = Z L Z 0 t 2 = Z L Z 0 4n 2N Z L Z 0 n 2N + 1 2 . (1.116)
We can note that in the special case where Z 0 = Z L , when looking at the energy (or power), the coefficient of transmission is then:

T = t 2 = 4n 2N (1 + n 2N ) 2 , (1.117)
which corresponds to expression (1) given by Beaky et al. in [START_REF] Beaky | Two-Dimensional Photonic Crystal Fabry-Perot Resonators with Lossy Dielectrics[END_REF]. In their article, N is the number of layers.

Here, N is the number of bilayers, so 2N is the number of layers.

Now, if we consider the case where Z 0 = Z 2 and Z L = Z 1 then the coefficient of transmission in energy becomes:

T = 4n 2N +1 (1 + n 2N +1 ) 2 .
(1.118)

Exploitation

We now consider the resonator constituted of a finite width defect surrounded by two symmetrical sets of bilayers. In such structure, the wave going from the input medium to the defect doesn't see an output medium with an impedance equal to Z C = Z 1 . In fact, as shown in section 1.3.8, the equivalent impedance can be complex, leading to an expression of the transmission coefficient T that becomes:

T = n 2N 4 (1 + n 2N ) 2 .
(1.119)

We know that, at the resonance, full transmission is reached, so the amplification factor in power compensates exactly the attenuation factors

(T 1 T 2 Q 2 = 1), so T 1 = T 2 = 1/Q. This leads, for n = Z 1 /Z 2 ≫ 1, to: Q = 4 Z 1 Z 2 N or ln Q = N ln Z 1 Z 2 + ln 4. (1.120)
For the three couples of materials with constants provided in table 1.1, the curve of ln(Q) vs. N given by equation (1.120) is compared to values of Q obtained in numerical experiments and plotted on Figure 1.11.

The matching between analytical and numerical results is excellent. Table 1.1 -Density, linear elastic constant and longitudinal elastic wave velocity for steel and silica.

Numerical Methods

We have given detailed equations of the transfer matrix method in order to provide several analytical results

for the study of one-dimensional PnCs. Those considerations will be useful in chapters 2 and 3.

We will now provide some details about three numerical methods for the study of 1D superlattices as well as information about their implementations in our programs.

Pseudospectral Method

The Pseudospectral (PS) method that we will use in chapter 3 to solve numerically the propagation of an acoustic wave in a heterogeneous medium is described in [START_REF] Bou Matar | Two-dimensional axisymmetric numerical simulation of supercritical phase conjugation of ultrasound in active solid media[END_REF]. The wave propagates in a 1D medium, in the z direction.

The following equations are implemented in a Fortran code. First is Newton's second law:

∂v ∂t = 1 ρ ∂τ ∂z , ( 1.121)
where v(z, t) is the particle velocity, t is the time, ρ is the density, τ(z, t) is the stress and z the position.

The stress is given by:

τ = Y (z, ε)ε where ε = ∂u ∂z , (1.122)
where Y (z, ε) is the Young's modulus and u(z, t) is the particle displacement, which is linked to the particle velocity by: ∂u ∂t = v.

(1.123)

In order to solve equation (1.121), the Pseudospectral method calculates the spatial derivatives using a staggered grid, as described in [START_REF] Özdenvar | Causes and reduction of numerical artefacts in pseudo-spectral wavefield extrapolation[END_REF]. The derivative ∂τ/∂z is calculated by multiplying each point of the Fourier transform of τ(z) over z by jke jk∆z/2 and taking the inverse Fourier transform of the result:

∂τ ∂z = F -1 jke jk∆z/2 F (τ(z)) , (1.124)
where F and F -1 are the forward and inverse Fourier transform, respectively, k is the wave number and ∆z the step of the grid. This way of calculating the space derivatives by a wavenumber multiplication in the wavenumber domain gives the PS method a greater accuracy and a better stability than that of the finite-difference methods [START_REF] Özdenvar | Causes and reduction of numerical artefacts in pseudo-spectral wavefield extrapolation[END_REF].

Similarly, the derivative ∂u/∂z is calculated as:

∂u ∂z = F -1 jke -jk∆z/2 F (u(z)) . (1.125)
In our code, the time stepping is performed using a staggered fourth order Adams-Bashforth method [START_REF] Ghrist | Staggered time integrators for wave equations[END_REF] by which stress and particle velocity are updated at alternating half time steps to integrate forward in time, as shown on Figure 1.12.

Finite Difference Time Domain

The Finite Difference Time Domain method is a simple method that consists in discretizing the time and space and solving the differential equations in an iterative manner.

1D FDTD

Solving the 1D elastic wave propagation is done in the following way. The space is discretized into a grid with step ∆x and the equations are solved at each time step ∆t. First, the strain ε is calculated from the displacement u for each grid point i:

ε i (t) = u i+1 (t) -u i (t) ∆x . (1.126)
Then, the stress τ is deduced from the strain:

τ i (t) = Y i (t)ε i (t).
(1.127) Discretizing Newton's second law, given in equation (1.121), for each grid point, we get:

v i (t + ∆t) -v i (t) ∆t = 1 ρ τ i (t) -τ i-1 (t) ∆x . (1.128)
So the speed v is updated with:

v i (t + ∆t) = v i (t) + F i (t) ∆t ρ , (1.129)
where we have introduced the force F given by:

F i (t) = τ i (t) -τ i-1 (t) ∆x . (1.130)
Finally, the displacement is updated by:

u i (t + ∆t) = u i (t) + v i (t + ∆t)∆t. (1.131)
The values of ε i , u i and v i are saved for post-calculation analysis and the whole calculation of those values is performed again for the next time step.

2D FDTD

Solving the 2D elastic wave propagation is done in the following way. The space is discretized into a 2D grid with steps ∆x and ∆y and the equations are solved at each time step ∆t. First, the stress τ is calculated from the displacement u for each grid point i, j:

τ xx,i,j (t) = C 11,i,j (t) u x,i,j+1 (t)-u x,i,j (t) ∆y + C 12,i,j (t) u y,i,j (t)-u y,i-1,j ∆x τ yy,i,j (t) = C 11,i,j (t) u y,i,j (t)-u y,i-1,j (t) ∆x + C 12,i,j (t) u x,i,j+1 (t)-u x,i,j ∆y τ xy,i (t) = C 44,i,j (t) u x,i+1,j (t)-u x,i,j (t) ∆x + u y,i,j (t)-u y,i,j-1 ∆y . (1.132)
The force F is calculated by: and the speed v is updated with:

F x,i,j (t) = τ xx,i,j (t)-τ xx,i,
v x,i,j (t + ∆t) = v x,i,j (t) + F x,i,j (t) ∆t ρ v y,i,j (t + ∆t) = v y,i,j (t) + F y,i,j (t) ∆t ρ . (1.134)
Finally, the displacement is updated by:

u x,i,j (t + ∆t) = u x,i,j (t) + v x,i,j (t + ∆t)∆t u y,i,j (t + ∆t) = u y,i,j (t) + v x,i,j (t + ∆t)∆t . (1.135)

Spectral Energy Density

For a specified wave vector k, the frequency spectrum of the spectral energy density is obtained by adding the square of the Fourier transform of the discrete temporal signal for every mass in every unit cell. The value of the SED is the average kinetic energy per unit cell, for a given wave vector and a given frequency. When a peak is found in the SED spectrum, it means that a vibrational eigen mode exists for wave vector k at that frequency. This way, the SED calculation allows one to quantify the phonon modes in the 1D superlattice.

Figure 1.13 -Mass-spring system as it is used by the SED method. The structure is composed of N unit cells of M masses.

If, as represented on Figure 1.13, we designate the total number of unit cells as N , the number of masses per unit cell as M and the total amount of time during which the particle velocity is recorded as T , then the expression for the SED calculation is:

Φ( k, ω) = 1 4πT N M i=0 m i         T 0 N j=0 v i (j, t)e i k• r 0 -iωt dt         2 , ( 1.136)
where v i (j, t) is the velocity of mass number i in unit cell j. The shape of the frequency spread for eigen mode k is represented by the Lorentzian function:

Φ( k, ω) = I 1 + ((ω -ω c ) /γ) 2 , (1.137)
where I is the magnitude of the peak, ω c is the frequency of the center of the peak and γ is the half-width of the peak at half-maximum. The lifetime for phonon mode k is defined as τ = 1/(2γ). Non-degenerate wave vector modes are written k i = 2πn i /(aN ) where a is the lattice constant and n i is an integer ranging from -N + 1 to N . The non-degenerate wave vector modes depend on the simulation size. The Spectral Energy Density method relies on another method, for instance FDTD, to calculate the particle velocities v i used in equation (1.136).

Conclusion

In this chapter, we have provided definitions and a formalism for the study of wave propagation in 1D PnCs.

We have introduced equations and analytical results as well as numerical methods that will be useful in the following chapters, namely the study of wave dissipation in a 1D bone structure (chapter 2) and a nonlinear elastic resonator (chapter 3).

Chapter 2

Dispersion of Elastic Waves in a Bone Structure

Abstract

In this chapter, we study the scattering of elastic waves in a bone structure composed of alternate collagen and hydroxy-apatite constituent layers. This structure is considered a 1D PnC and we characterize it by obtaining its dispersion diagram and showing how its band structure evolves with the amplitude of the displacement. The dissipation mechanism of elastic waves due to the nonlinear behavior of the collagen with respect to the strain is explained by the mean of a spectral energy density calculation.

It is essential to understand how the composite structure of bone and teeth reacts to dynamics loads which are often responsible for fracture. This is why we study how the nonlinear behavior of collagen in equilibrium with water provides a means of filling vibrational bad gaps. This band-gap filling facilitates the propagation of vibrations in a larger range of frequencies, limiting the risk of failure. This work has been published in the Journal of the Mechanical Behavior of Biomedical Materials [START_REF] Guerder | Multi-phonon scattering processes in one-dimensional anharmonic biological superlattices: Understanding the dissipation of mechanical waves in mineralized tissues[END_REF].

Introduction

Mineralized biological tissues, such as bone and tooth, are hierarchical composite structures composed of a stiff hydroxyapatite (HAP) mineral phase, a compliant proteinaceous collagen phase, and water. At the nanoscale, bone and teeth are constituted of a periodic assembly of alternating regions of collagen and HAP in a hydrated environment with a repeat unit cell size of 67 nm [START_REF] Ten | Oral Histology: Development, Structure and Function[END_REF]. This periodic composite structure, forming a one dimensional (1D) superlattice, is believed to be responsible for the remarkable strength and toughness of these biological materials [START_REF] Jager | Mineralized Collagen Fibrils: A Mechanical Model with a Staggered Arrangement of Mineral Particles[END_REF][START_REF] Currey | Bones: Structure and Mechanics[END_REF][START_REF] Deymier-Black | Synchrotron X-Ray Diffraction Study of Load Partitioning During Elastic Deformation of Bovine Dentin[END_REF]. At the micrometer scale, mineralized tissues exhibit a large network of interconnected porosity, tubules in dentin and canaliculi and lacunae in bone, which allow for the transfer of nutrients, waste, and water throughout the tissue [START_REF] Boyde | An Electron Microscope Study of Fractured Dentinal Surfaces[END_REF][START_REF] Dillaman | Fluid Movement in Bone: Theoretical and Empirical[END_REF]. This porosity allows bone and teeth to remain in equilibrium with the water, maintaining hydration of the tissues. Water molecules exhibit a variety of different interactions with the HAP and collagen including the formation of water-bridges within the collagen helix, filling channels within the HAP, and surface hydration of the collagen and HAP phases [START_REF] Bertinetti | Surface Structure, Hydration, and Cationic Sites of Nanohydroxyapatite: Uhr-Tem, Ir, and Microgravimetric Studies[END_REF][START_REF] Simone | Role of Hydration in Collagen Triple Helix Stabilization[END_REF][START_REF] Ravikumar | Region-Specific Role of Water in Collagen Unwinding and Assembly[END_REF]. Three-point bend and notch testing indicate that hydration has a significant impact on the mechanical properties of mineralized tissues resulting in increased elastic moduli as well as decreased toughness and loss of plastic behavior [START_REF] Kahler | Fracture-Toughening Mechanisms Responsible for Differences in Work to Fracture of Hydrated and Dehydrated Dentine[END_REF][START_REF] Kruzic | Crack Blunting, Crack Bridging and Resistance-Curve Fracture Mechanics in Dentin: Effect of Hydration[END_REF][START_REF] Nyman | The Influence of Water Removal on the Strength and Toughness of Cortical Bone[END_REF]. Hydration increases the nonlinear behavior of collagen as well as its elastic modulus while increasing its toughness [START_REF] Grant | Effects of Hydration on the Mechanical Response of Individual Collagen Fibrils[END_REF][START_REF] Gevorkian | Stabilization and Anomalous Hydration of Collagen Fibril under Heating[END_REF]. It is theorized that these changes in the collagen mechanical behavior due to water content are the major cause for the changes in overall tissue mechanics.

Although these tests provide important information about how hydration affects material properties, they provide no information about the dynamic properties of these tissues.

The number one cause of fracture in both bones and teeth is trauma or impact. To avoid fracture, the energy from sudden impacts must be dissipated in order to limit the formation of stress concentrators that can lead to tissue failure and fracture. It is therefore essential to understand how the composite structure of bone and teeth reacts to these dynamics loads which are so often responsible for fracture. In this study, it is theorized that the nonlinear behavior of collagen in equilibrium with water provides a means of filling vibrational band gaps that arise from the periodicity of the HAP/collagen structure assembly. Under high load (high deformation) conditions, the nonlinearity of the mechanical response of the collagen-water system may open up multiphonon scattering channels leading to a filling of the band gaps in the vibrational band structure of the HAP/collagen superlattice. We model the collagen/water system within the context of the thermodynamics of stressed open systems. This model leads to nonlinear stress-strain relationships of the collagen. The nonlinear model of the open collagen/water systems is incorporated into a dynamic model of the HAP/collagen superlattice. The propagation of elastic waves through that superlattice is investigated using the finite difference time domain (FDTD) method in conjunction with the spectral energy density (SED) method. Vibrational (phonon) wave band structures of the HAP/collagen/water system show that at high amplitudes, vibrational waves can interact with each other through multiphonon scattering channels that can fill the band gaps inherent in the band structure of elastic superlattices. This band-gap filling facilitates the propagation of vibrations in a larger range of frequencies, providing an effective mechanism for dissipation of mechanical energy, thus limiting the risk of failure. This chapter is organized as follows. In section 2.3, we introduce the model of the open collagen/water system. This system is modeled within the context of the thermodynamics of stressed solid solutions. This model results in the formulation of the nonlinear stress-strain response of the collagen/water solid, i.e. a strain dependent elastic modulus. Section 2.4 presents the model of the HAP/collagen periodic structure in the form of a 1D superlattice as well as the methods of FDTD and SED that are employed to investigate the dynamic response of the superlattice. In particular we focus on the calculation of the vibrational/phonon band structure of the superlattice as a function of the energy (amplitude) of the propagating elastic waves.

In section 2.5, we report the results of the calculations and provide an analysis of the effect of the nonlinear elastic modulus. In particular, it is shown that high amplitude (energy) waves can interact with each other through multiple phonon scattering processes. These interactions lead to the opening of new channels for the dissipation of the elastic energy over a wider range of frequencies compared to the case of low energy waves.

Finally, conclusions are drawn in section 2.6 as to the relationship between the observed behavior and bone fragility due to decreased hydration.

Models

Thermodynamics of a Stressed Solid Solution

To address the problem of the mechanical behavior of bone material in the presence of water, we develop the chemo-mechanical equations of states of materials that can adsorb fluids under stress based on the work of Larché and Cahn [START_REF] Larché | Thermodynamics of stressed solids[END_REF][START_REF] Larché | The interactions of composition and stress in crystalline solids[END_REF]. The total internal energy of the material is obtained as an integral of an internal Helmholtz energy density f ′ :

E = V ′ f ′ dV ′ , (2.1)
where the energy density is given by:

f ′ = f ′ (T , ε, ..., c ′ i , ...). (2.2)
The prime indicates that all densities are relative to the reference state for measuring strain. T , ε, and c ′ i are the temperature, strain and molar density of chemical constituent i. We consider K variable chemical species in the chosen materials. The differential form of equation (2.2) takes the form:

df ′ (T , ε ij , c ′ I , ...) = s ′ (T , ε ij , c ′ I , ...)ds ′ + σ ij (T , ε ij , c ′ I , ...)dσ ij + I,K M I,K (T , ε ij , c ′ I , ...)dc ′ I .
(

The functions: s ′ (T , ε ij , c ′ I , ...), σ ij (T , ε ij , c ′ I , ...) and M I,K (T , ε ij , c ′ I , ...) are the density of entropy, stress and diffusion potential equations of state. The diffusion potential is used when considering a substitutional solid solution that constrains the molar densities according to: c ′ 1 + ... + c ′ I + c ′ K = c ′ 0 where c ′ 0 is the density of substitutional sites of the different chemical species. Introducing mole fractions X I = c ′ I /c ′ 0 the diffusion potentials are therefore defined as:

∂f ′ ∂X I T ,ε,X I K = c ′ 0 M I,K with I = 1, ..., K -1. (2.4)
By choosing K as a dependent chemical specie, one may treat the problem with only K -1, independent variables. Let us now simplify the problem to a binary solution i.e. I = 1 and K = 2. Note that the substitutional solid solution representation is equivalent to an interstitial solid solution if 1 represents the interstitial specie (e.g. water) and 2 the interstitial sites (i.e. available sites for water in a collagen matrix).

Denoting by X the composition in specie 1 and using 2 as the dependent specie, the diffusion potential becomes:

M 1,2 = 1 c ′ 0 ∂f ′ ∂X 1 T ,ε . (2.5)
It is convenient to define the free energy density Φ ′ by the Legendre transformation where strain is replaced by stress as variable:

Φ ′ = f ′ -σ ij ε ij . (2.6)
The differential form of the density Φ ′ is:

dΦ ′ = -ε ij dσ ij -s ′ dT + c ′ 0 M 1,2 dX 1 . (2.7)
From this relation we deduce the following Maxwell relation:

-c ′ 0 ∂M 1,2 ∂σ ij T ,X 1 = ∂ε ij ∂X 1 T ,σ kl . (2.8)
In the case of a binary solution which elastic properties depend on the composition, the right hand side of equation (2.8) has to be written as:

∂ε ij ∂X 1 T ,σ kl = ∂ε c ij ∂X 1 + ∂S ijkl ∂X 1 σ kl , ( 2.9)
where we have dropped the subscript in the differentials for the sake of simplifying the notation. ε c ij are the components of the chemical strain and S ijkl are the components of the compliance tensor. The chemical strain is stress free and is only associated with the expansion or contraction of the material upon a change in composition. To simplify the notation, we take η ij = ∂ε c ij /∂X 1 where the linear coefficients η ij are the components of the chemical expansion coefficient tensor. The simplest relationship between the change in composition (X -X 0 ) and the chemical strain is therefore:

ε c ij = (X -X 0 )η ij δ ij , (2.10)
where δ ij is the Kroenecker symbol and X 0 is the composition of the reference state for measuring strain. To obtain the second term in equation (2.9), we have used Hooke's law:

s ij = C ijkl (e kl -e c kl ), (2.11)
or:

ε ij -ε c ij = S ijkl σ kl .
(2.12)

C ijkl are the components of the stiffness tensor. The quantity ε m ij = ε ij -ε c ij is the mechanical strain. Inserting equation (2.9) into equation (2.8) and after integration, the diffusion potential becomes:

c ′ 0 M 1,2 = - ∂ε c ij ∂X σ ij - 1 2 ∂S ijkl ∂X σ kl σ ij + φ(X), (2.13)
where φ is some unknown function of composition. This unknown function is eliminated by choosing a hydrostatic state of pressure P as reference state. Inserting equation (2.10) into equation (2.13), the diffusion potential equation of state for the binary solution is now given as:

M 1,2 = µ 1 (T , P, X 2 ) -µ 2 (T , P, X 1 ) -V ′ 0 η ij σ ij -η kk P - 1 2 ∂S ijkl ∂X 1 σ kl σ ij + 1 2 P 2 ∂S ijkl ∂X 1 δ kl δ ij , (2.14)
where V ′ 0 = 1/c ′ 0 . In equation (2.14), we have introduced the chemical potentials of species 1 and 2. Under hydrostatic pressure, the diffusion potential is nothing but a difference in chemical potential. The chemical potentials are defined as:

µ i (T , P, X i ) = 1 c ′ 0 ∂f ′ ∂X 1 T ,ε,X i j .
(2.15)

The condition for chemical equilibrium of the binary solid solution in contact with a binary fluid solution is determined by the conservation of the diffusion potential,

M 1,2 (σ, x) = µ 1 -µ 2 and µ F 1 = µ 1 and µ F 2 = µ 2 .
(2.16)

We now make the temperature dependency implicit and drop T from the equations. Subtracting the diffusion potentials of a stressed solid solution, M 1,2 and a solid solution under hydrostatic pressure, M 1,2 yields:

M 1,2 -M 1,2 = µ F 1 -µ F 1 -(µ F 2 -µ F 2 ).
(2.17)

µ F 1 , µ F 1 , µ F 2 , µ F
2 are the chemical potentials of the species 1 and 2 in the fluid when the solid is subjected to a stress or to a hydrostatic pressure only, respectively. Provided that the fluid behaves like a chemical reservoir, this is the case if we consider the fluid to be a reservoir of water only, the difference on the right hand side of equation (2.16) becomes identically equal to 0 and the diffusion potentials of the stressed and unstressed system are equal. Note that the diffusion potential of the stressed solid is evaluated at the composition X while the hydrostatic diffusion potential is the difference in chemical potential of the species 1 and 2 at equilibrium hydrostatic composition X 0 . The diffusion potential of the stressed system can be calculated from:

M 1,2 = µ 1 (P, X 0 1 ) -µ 2 (P, X 0 2 ). (2.18)
Equation (2.18) is sufficient to solve for the change in composition, X -X 0 , of the solid to maintain equilibrium under stress with the fluid reservoir. For this, we use the equation of state 2.14.

µ 1 (T , P, X 0 1 ) -µ 2 (T , P, X 0 2 = 1 -X 0 1 ) = µ 1 (T , P, X 1 ) -µ 2 (T , P, X 2 = 1 -X 1 ) -V ′ 0 η ij σ ij -η kk P - 1 2 ∂S ijkl ∂X 1 σ kl σ ij + 1 2 P 2 ∂S ijkl ∂X 1 δ kl δ ij . ( 2
f m = RT (X 1 ln X 1 + (1 -X 1 ) ln(1 -X 1 )) + Ω(1 -X 1 )X 1 .
(2.20)

Ω represents the interaction energy between the species 1 (water) and 2 (water sites in collagen). In this case, the difference in chemical potential is given by:

µ 1 (T , P, X 1 ) -µ 2 (T , P, X 2 = 1 -X 1 ) = RT ln X 1 1 -X 1 + Ω(1 -2X 1 ). (2.21)
We now assume that the equilibrium mole fraction of water in collagen in absence of stress is X 0 1 = 0.5. Molecular dynamics simulations of the interaction between collagen and water indicate that the occupancy of water sites, both internal and external to the collagen helix, is approximately 40-60 % [START_REF] Simone | Role of Hydration in Collagen Triple Helix Stabilization[END_REF][START_REF] Ravikumar | Region-Specific Role of Water in Collagen Unwinding and Assembly[END_REF]. Under this condition, in absence of stress, the difference of chemical potential given by equation (2.21) is zero. We further assume that the elastic coefficients of collagen are independent of water content. This is justified based on the fact that the speed of sound and the density of polymer-based media and water are very similar.

Finally, reducing the problem to a 1D one, we rewrite equation (2.21) in the form:

RT ln X 1 1 -X 1 + Ω(1 -2X 1 ) = V ′ 0 ησ. (2.22)
To obtain equation (2.22) we have also neglected the hydrostatic pressure P compared to the stress. The chemo-mechanical problem of a collagen solid matrix in equilibrium with a water reservoir then can be solved by eliminating composition between equation (2.22) and equation (2.23):

ε = σ Y + η(X 1 -X 0 1 ). (2.23)
Equation (2.23) is the 1D version of equation (2.12) combined with equation (2.10). There, Y is Young's modulus of the collagen with the stress free water content, X 0 1 = 0.5.

Nonlinear Young's Modulus of Collagen

Bowman has measured the stress versus strain relationship in demineralized bovine bone under uniaxial tension and open conditions (i.e. in equilibrium with water), from which we extracted the nonlinear stressstrain curve for tensile stress [START_REF] Bowman | The tensile behavior of demineralized bovine cortical bone[END_REF]. The part of the curve for the negative values of the strain is the symmetrical of the positive part, as it is in the equations above. Figure 2.1 shows a fit to this data obtained by eliminating composition between equations (2.22) and (2.23). We found that best fit is achieved for the following conditions: η = 0.12, Y = 7.75 × 10 8 Pa, Ω = 5448.9 J mol -1 and V ′ 0 = 4.56 × 10 -4 m 3 mol -1 . These values are in agreement with values of η = 0.09 and Y = 0.5 -1 GPa, measured for type 1 collagen [START_REF] Iii | Heterogeneity of the Mechanical Properties of Demineralized Bone[END_REF][START_REF] Heim | Determination of the Elastic Modulus of Native Collagen Fibrils Via Radial Indentation[END_REF][START_REF] Fullerton | Micro-Ct Dilatometry Measures of Molecular Collagen Hydration Using Bovine Extensor Tendon[END_REF].

Equations (2.22) and (2.23) which describe the open-system stress-strain relationship by elimination of the compositional variable X, include five parameters. After fitting the model to the available experimental stress-strain curve, one obtains three of the parameters that are in excellent agreement with published data for collagen. The other two parameters take on values that are physically meaningful.

From the fit, the Young's modulus and the elastic energy are written as Taylor series up to the 8th order.

Y (ε) = a 0 + a 2 ε 2 + a 4 ε 4 + a 6 ε 6 + a 8 ε 8 .
(2.24)

From expression (2.25), we obtain the stress σ = Y (ε)ε and the elastic energy This formula is valid for |ε| 0.09. In our simulations, we will always ensure that the strain will not overpass this value. Collagen has a rubber-like behavior, which explains why values of strain as high as 0.12 can be reached [START_REF] Bowman | The tensile behavior of demineralized bovine cortical bone[END_REF].

E(ε) = ε 0 σ(ε)dε: E(ε) = a 0 2 ε 2 +
The Young's modulus, stress and elastic energy versus the strain obtained with this Taylor series are represented in Figure 2.1.

The experimental data from [START_REF] Bowman | The tensile behavior of demineralized bovine cortical bone[END_REF] only contains information for the positive values of strain. We make the assumption that the relation is symmetrical for the negative values of strain.

Methods

Matrix Transfer Method

As demonstrated in chapter 1, section 1.3.2, the elastic band structure of an infinite linear AB superlattice can be calculated analytically by [START_REF] Camley | Transverse elastic waves in periodically layered infinite and semi-infinite media[END_REF]:

cos(2kL) = cos(k A L) × cos(k B L) -0.5 × sin(k A L) × sin(k B L) × Z A Z B + Z B Z A . ( 2

.26)

Where L = 33.5 nm is the length of a layer, Z i is the impedance of material i, ω is the frequency, k i = ω/c i is the wave vector and c i is the particle velocity in material i.

Equation (2.26) provides an analytical solution for the band structure of the PnC in the linear case. We use Matlab to implement a this equation in a transfer matrix method. We evaluate the expression for 5000 frequencies between zero and 20 GHz.

SED-FDTD Method Description

We now consider the 1D model illustrated in Figure 2.2 as representative of the nanoscale, periodic assembly of alternating regions of collagen and HAP in bone and teeth. The finite super lattice is composed of N = 2560 masses of materials A and B connected by nonlinear springs. Material A is chemically inert and represents HAP. Material B represents the open collagen/water system i.e. collagen that can adsorb water from a reservoir. The behavior of elastic modes in this model is simulated using a modified Spectral Energy Density -Finite Difference Time Domain (SED-FDTD) method. The SED-FDTD method has been used recently with success to simulate the propagation of elastic wave propagation in PnCs [START_REF] Thomas | Predicting Phonon Dispersion Relations and Lifetimes from the Spectral Energy Density[END_REF]. In its original version, this method solves the springs and masses equation by discretizing time and space (x-axis) and by replacing derivatives by finite differences. We extend this approach to include chemo-mechanical effects by using the nonlinear stress-strain relation described in section 2.3, equation (2.24). Figure 2.2 -Schematic representation of the system. The unit cells (UC) are heterogeneous media composed 32 masses of a material A mimicking hydroxy apatite that does not change composition when subjected to a stress and 32 masses of material B. B is a collagen-like material that can interact with a chemical reservoir of water it is embedded into. Each mass is connected to its two neighbors by two springs. A spring is linear (respectively nonlinear) if the mass on its left is composed of material A (resp. material B). Periodic boundary conditions are imposed at the end of the system, which contains 40 UCs.

The nonlinear spring constant depends on the water content in accordance with the nonlinear stress-strain relation of the collagen/water open system. Owing to the 1D nature of the system, we assume that lateral diffusion of water into material B (collagen) is much faster than the acoustic wave and that material B is always at equilibrium with respect to chemical composition. Under this condition, as explained in chapter 1, section 1.4.2, the 1D mass-spring equation is given by:

m i ∂ 2 u i ∂t 2 = f i + f i+1 with m i = ρ i dx 3 , (2.27)
where t is time, ρ i is the mass density of mass i and m i its mass, u i (t) its displacement. The displacement is related to the mechanical strain by ε = ∂u/∂x and σ is the stress. In 1D, we describe the stress-strain relation by:

σ = Y (ε)ε, (2.28)
or the force-displacement relation by:

f i = k i du i = ρ i c 2 L,i ∆xdu i , (2.29)
where k i = ρ i c 2 L,i ∆x and c L,i = Y i (ε)/ρ i and ∆x is the distance between two adjacent masses. Equation (2.27) is solved discretely and takes the form:

ρ i ∆x 3 ∂ 2 u i ∂t 2 = ρ i c 2 L,i (u i+1 -u i ) -ρ i-1 c 2 L,i-1 (u i -u i-1 ) ∆x. (2.30)
In equation (2.30) we assume that the mass density is independent of composition. This is the case in material A. Since the mass density of water and collagen are similar, the mass density of discrete points in material B are also taken to be constant.

Structure

The structure is composed of N = 40 unit cells each of them composed of M = 64 masses (32 of them composing the material A and the 32 others composing the material B). Hence, the total number of masses is

M × N = 2560.
Each mass is connected to its two neighbors by two springs. A spring is linear (respectively nonlinear) if the mass on its left is composed of material A (resp. material B). This structure is represented on Figure 2.2. The nonlinear spring constant is calculated using equation (2.28) where Y (ε) is the Taylor expansion described in section 2.3.

Simulation parameters

We employ a spatial grid with mesh size ∆x = 1.04 nm. The time integration step, ∆t, is given by ∆t = ∆x/(40c A ) = 6.0 × 10 -15 s where c A is the speed of sound in medium A. The total simulation time is 2 22 ∆t. Periodic boundary conditions are applied at the free ends of the homogeneous regions to simulate an infinite superlattice. The thickness of the segments of material A and B are designated by L A and L B . For the sake of simplicity we take: L A = L B = 33.5 nm. Hence the length of a unit cell is a = 67 nm.

The density, Young's modulus and elastic wave velocity constants used in our model for hydroxy apatite and collagen are given in table 2.2. Table 2.2 -Density, Young's modulus and elastic wave velocity for hydroxy apatite and collagen.

Layer

We use the SED-FDTD method to calculate the presence of acoustic modes in the crystal for different wave vectors. In particular, the SED method, described in chapter 1, section 1.4.3 enables the projection of the vibrations of a structure onto a set of plane waves characterized by a wave vector (wave number in 1D) and a frequency. With this method, an initial random displacement is applied to each of the 2560 masses which imparts an initial potential energy to the structure. Then, the system is free to evolve during the total number of time steps (2 22 ) and the speed of each mass is recorded during the last 2 21 time steps. For each mass, this speed is projected on a considered wave number and a Fourier transform provides the frequency distribution of the energy in the solid. This energy is averaged over the 64 masses of each unit cell. This operation is repeated 200 times for which the results are averaged, in order to allow a sufficient variety of initial energy distribution. The SED-FDTD method is a way of knowing how often a phononic mode is visited.

The dimension chosen for the system, 40 unit cells, allows the whole operation to be repeated over

Time resolution and convergence

We compare the time resolution of our calculation with the diffusion time for water in collagen. If we approximate the one-dimensional system by a wire with a L = 0.5 nm cross section, and if we consider a diffusion coefficient for water in collagen of D = 2.5 × 10 -5 cm 2 s -1 [START_REF] Westover | Collagen hydration: Pulsed nuclear magnetic resonance studies of structural transitions[END_REF], the time for diffusion through the wire is on the order of t = L 2 /(6D) = 1.67 × 10 -11 s. If we compare this time to the maximum frequency that we study, namely 20 GHz, we obtain a characteristic time of 5 × 10 -11 s. This time is significantly larger than the time for diffusion through the model system. So, even at high frequencies, we can assume that water content of the system is at equilibrium at all time. This assumption is even more valid when considering low frequency modes.

As for the convergence of the method, the SED-FDTD calculation was performed using 16, 32 and 64 masses per unit cell in the limit of small displacements. We observed that the band diagram obtained with 16 masses per unit cell almost matched the linear results below 10 GHz but not above this frequency.

The band diagram obtained with 64 masses per unit cell matched almost perfectly the linear band diagram as shown on Figure 2.3. The band diagram obtained with 32 masses per unit cell was an intermediary result, with a match better than that for 16 but worse than that for 64 masses per unit cell. Hence, the value of 64 masses per unit cell was chosen for all the SED-FDTD calculations.

Results

Matrix Transfer Method

The band diagram of the structure calculated with the matrix transfer method is plotted on Figure 2.3. As one can see, five transmission modes are obtained, below 20 GHz. They are separated by absolute band gaps.

In the frequency ranges corresponding to these band gaps, no wave can propagate: only evanescent modes can appear.

SED-FDTD Method

The spectral energy density of the structure for an initial displacement of ∆x/150 (low amplitude) and that for an initial displacement of ∆x/45 (large amplitude) are provided on Figure 2.4.

At low amplitude, the spectrum shows very thin peaks at specific frequencies corresponding to the vibrational modes that are found when the transmission is solved analytically (Figure 2.3). The system behaves according to linear elasticity: it shows a harmonic behavior. In this regime, the elastic coefficient (Young's modulus) of the collagen (B) region retains its value at small strain. For all the frequency range (0 to 20 GHz), an excellent agreement between this nonlinear result and the linear result obtained with equation (2.26) is observed.

At large amplitude, the system becomes strongly nonlinear. The peaks broaden significantly. A filling of the band gaps is observed, as shown on Figure 2.6. This band gap filling is more visible at higher frequencies (10 to 20 GHz), which allows us to think that as frequency increases, the band gaps will be more and more Figure 2.3 -Band structure of the AB superlattice. The wavevector is expressed in reduced units of 1/d. This is calculated using the matrix transfer method. filled. At higher frequencies (hundreds of GHz), this would mean that the energy is dissipated as thermal phonons.

We now study the influence of each term of the Taylor series of the Young's modulus. Taking each term into account or not will result in different energy functions. The terms of the Taylor series of the energy function are represented on Figure 2.5. The peak broadening and frequency shift are interpreted as follows.

If we limit the Young's modulus function of strain to its constant term, which is exactly the linear case, then the Taylor series expansion of the elastic energy function is limited to its second order term. In that case, the spectrum from the SED-FDTD calculation only contains the "primary" frequencies which are that of a linear system. This result is shown on Figure 2.6. Now, if we limit the Young's modulus function of strain to its terms up to order 2, the energy function contains a term of order 2 and a term of order 4. As shown on Figure 2.6, the system simulated in the SED-FDTD program behaves as a nonlinear system with the frequency shift and band gap filling described above.

If we limit the Young's modulus function to its terms up to order 4 (resp. 6), the energy function contains terms of order 2, 4 and 6 (resp. 2, 4, 6, 8). Taking those terms into account (see Figure 2.6) does not change the behavior of the system, which proves that the second order term of the Young's modulus (fourth order term of the energy) plays the major role for the nonlinear behavior and the terms of higher order do not play 1, namely the integral of the Taylor series decomposition of the stress as a function of strain. The dotted and dashed lines are the terms of orders 2, 4, 6 and 8 of this Taylor series. The term of order 6 is particularly interesting as it is gives a concave down curve characteristic of an unstable system (its second derivative is negative). a significant role. Hence, we can assert that this nonlinear behavior associated with a fourth-order function in strain elastic energy is mainly due to four-wave (four-phonon) interactions. Indeed, considering that the dynamic strain is a linear superposition of plane waves with different wave number and frequencies and raising this superposition to the fourth power to estimate the energy, leads to four-wave interactions. These scattering interactions would conserve momentum and frequency and may involve a variety of processes such as the splitting of a single phonon into three others, the scattering of two phonons forming two others, etc. [START_REF] Manktelow | Multiple scales analysis of wave-wave interactions in a cubically nonlinear monoatomic chain[END_REF][START_REF] Swinteck | Phonon Scattering in One-Dimensional Anharmonic Crystals and Superlattices: Analytical and Numerical Study[END_REF]. Those new phonons interactions offer more channels for the dispersion of mechanical energy.

We can relate the behavior of our nonlinear superlattice to the behavior of a multiple-well mass-spring system. A multiple-well system is a system for which the representation of the elastic energy as a function of strain is not parabolic but is a superposition of infinity of parabolas. To each parabola corresponds a single value of the spring constant. Hence, representing a non-quadratic energy map by a multiplicity of parabolic wells is equivalent to introducing a continuum of values of spring constants. At high amplitudes, waves will sample wider ranges of strain values, effectively visiting the multiplicity of energy wells. The band diagram will show a continuity of modes above the "primary" frequencies of the linear system. Those bands broaden and finally fill the gaps. The passing bands of this system will also be shifted to higher frequencies. Figure 2.6 -SED spectrum for a specific wave vector. Spectral energy density of the phononic crystal structure with an initial displacement of ∆x/45 for wave vector π/a. Black line: with the Young's modulus truncated at the order zero (constant value, linear case); red dashes: with the Young's modulus truncated at the order 2; green dots: with the Young's modulus truncated at the order 4. Curves with Young's modulus truncated at orders 6 and 8 are not included since they do not carry significant changes compared to order 4. Inset: magnified highlight of the band gap filling in the frequency range between 0 and 8 GHz.

In our case, the reason why the spring constants increase as a function of the magnitude of the strain is that the model we chose for the Young's modulus is symmetrical and has only even orders in its Taylor series expansion. If an asymmetrical function was used for the Young's modulus, its Taylor series expansion would contain odd order terms. In that case, nonlinear modes with frequencies less than the "primary" frequencies would appear. We would also observe a shift of the passing bands to lower frequencies.

Conclusion

A 1D PnC composed of a superlattice of alternate collagen and hydroxy-apatite constituent layers is studied. This serves as a model of mineralized biological tissue. The collagen layers are treated as an open system that can absorb or desorb water. This collagen/water open system results in a nonlinear Young's modulus.

This open system model is fitted to experimental data. We study the dynamic response of the nonlinear superlattice to the propagation of elastic waves. A FDTD method coupled to a SED calculation is performed to obtain the wave band structure. In the linear limit, the band structure is composed of passing bands separated by band gaps that forbid the propagation of elastic waves over some ranges of frequency. The nonlinearity of the collagen layers in the superlattice gives rise to multi-wave (phonon) scattering processes that lead to partial band gap filling. Multi-phonon scattering processes constitute ways of opening new channels for the dissipation of mechanical energy. This mechanism for mechanical energy dissipation is the direct consequence of the hydration of the collagen.

The risk of dentin and bone fracture has been shown to increase with age. During aging, there is also a marked decrease in the level of interconnected porosity in both tissues [START_REF] Kinney | Age-Related Transparent Root Dentin: Mineral Concentrations, Crystallite Size and Mechanical Properties[END_REF][START_REF] Weinstein | Endogenous Glucocorticoids Decrease Skeletal Angiogenesis, Vascularity, Hydration, and Strength in Aged Mice[END_REF]. This filling of the porosity is associated with decreased fluid flow and therefore decreased hydration of the tissues. It has previously been theorized that this decrease in hydration may in part be responsible for the increased fragility of bone and dentin with age; however, the associated mechanism remained unclear [START_REF] Manolagas | Estrogen-Centric to Aging and Oxidative Stress: A Revised Perspective of the Pathogenesis of Osteoporosis[END_REF]. The results presented here provide a possible explanation for the relationship between bone fragility and decreased hydration.

The presence of interstitial water in the collagen phase increases its plastic behavior, but more importantly allows it to act nonlinearly. When mineralized tissues undergo a traumatic event, the nonlinearity of the collagen allows for band gap filling, thus allowing larger range of frequencies to propagate mechanical waves through the material. This ease of propagation diminishes the formation of stress concentrations by vibrational modes that would otherwise be banned from propagation, reducing the risk of fracture. In turn, a decrease in hydration will lead the collagen to act more linearly, which is done here by limiting the Young's modulus to its constant term, and as a result limit the number of channels that can dissipate elastic energy and frequencies which are allowed to propagate, thus increasing the development of stress concentrations and the possibility of fracture. This project is an application of the study of PnCs outside its usual fields (acoustics, materials science).

The knowledge and techniques developed for solving the equations of propagation of the elastic waves can have direct applications to fields such as biomechanics and biology.

Chapter 3

A Nonlinear Elastic Resonator

Abstract

In 2002, Soljačić et al. described a resonator based on cubic nonlinearities in a photonic crystal with a bistable behavior [START_REF] Soljačić | Optimal bistable switching in nonlinear photonic crystals[END_REF]. The existence and behavior of this system follow the generic rules described by Landau and Lifshitz [START_REF] Landau | Physique théorique: Mécanique. Physique théorique[END_REF]. The purpose of the study presented in this chapter is to design a similar resonator for quadratic and cubic nonlinearities in a PnC.

We present a nonlinear PnC structure capable of elastic bistable switching. The structure is comprised of a one-dimensional elastic PnC composed of alternate steel and silica constituent layers and a nonlinear steel defect located at its middle. Numerical simulations and analytical modeling of the transmission of a longitudinal elastic wave through the crystal have been performed. These calculations show that the central nonlinear inclusion constitutes a resonating cavity which is able to commute the elastic wave depending on the particle velocity. The key parameters for this effect are the nonlinear elastic constants of the material, the difference between the excitation frequency and the resonance frequency of the cavity as well as its quality factor. The PnC structure reproduces for elastic waves what has previously been observed for electromagnetic waves in a nonlinear photonic crystal.

Introduction

In recent years, thermal and elastic switches, consisting of granular PnCs, which are arrays of packed elastic spheres that interact through Hertzian contacts, have been studied and proposed for different applications, such as controlling and sensing the flux of energy. A bifurcation-based switch and rectificator has been described [START_REF] Boechler | Bifurcation-based acoustic switching and rectification[END_REF] and nonlinear resonances in diatomic granular chains have been observed [START_REF] Cabaret | Amplitude-dependent phononic processes in a diatomic granular chain in the weakly nonlinear regime[END_REF]. These models and experiments are expected to lead to the creation of new devices, such as amplifiers and attenuators. These new functionalities are associated with the amplitude-dependent behavior that results from the geometry of the structures, especially the Hertzian contacts between the particles, but do not depend from the nonlinear elastic constants of the constitutive materials. Very few studies have been performed on the use of the 49 nonlinearities of the materials to control the behavior of PnCs. More work is needed on this aspect of PnCs to extend the realm of their applications.

Recently, the use of constitutive materials with nonlinear elastic properties as well as the magnetoacoustic effect have been proposed to tune the band-gaps of two-dimensional PnCs [START_REF] Robillard | Tunable magnetoelastic phononic crystals[END_REF][START_REF] Bou Matar | Band gap tunability of magneto-elastic phononic crystal[END_REF][START_REF] Gei | Band-gap shift and defect-induced annihilation in prestressed elastic structures[END_REF]. Modeling and simulations of systems using these properties and effect have been reported [START_REF] Vasseur | Band structures tunability of bulk 2D phononic crystals made of magneto-elastic materials[END_REF]. In optics, a nonlinear photonic crystal, consisting of a resonating cavity with a nonlinear optical index, capable of performing optimal bistable switching has been described by an analytical model accompanied with numerical experiments [START_REF] Soljačić | Optimal bistable switching in nonlinear photonic crystals[END_REF]. The purpose of this chapter is to determine whether we can obtain in elastodynamics similar results with those observed in optics, namely a bistable behavior relatively to the particle velocity, due to nonlinearities of the elastic constants of the material constituting the resonating cavity. We consider here a structure composed of a one-dimensional (1D) PnC containing a nonlinear cavity. We observe and demonstrate without ambiguity the switch behavior in transmission of this structure. This chapter is organized as follows: first, section 3.3 draws the basis of the study of nonlinear oscillators, starting from the equations of Landau and Lifshitz [START_REF] Landau | Physique théorique: Mécanique. Physique théorique[END_REF]. Then, section 3.4 provides the analytical model for the nonlinear resonator. Section 3.5 describes the structure of the PnC that we propose and the numerical experiment that was performed. Section 3.6 presents these numerical results. Then, in section 3.7, a complete analytical interpretation and a theoretical model of the switching phenomenon is provided. In section 3.7.1, the effects of the different physical and geometrical parameters on the switching capability are analyzed. Finally, conclusions are drawn concerning the possibilities opened by this study, such as the tunability of PnCs through the magnetoacoustic effect and the design of two-dimensional elastic resonators.

Nonlinear Oscillators

This section describes the equations of quadratic and cubic nonlinear oscillators. For more details, an excellent reference for this subject is Landau's book [START_REF] Landau | Physique théorique: Mécanique. Physique théorique[END_REF].

Base Equation

The base equation of nonlinear oscillators is that of a Duffing oscillator:

ẍ + 2λ ẋ + ω 2 res x = F cos(ω 0 t) -αx 2 -βx 3 , (3.1)
where x is the displacement, α and β the quadratic and cubic nonlinear constants respectively, λ is the damping coefficient, F is homogeneous to an acceleration, ω 0 is the pulsation of the forced oscillation and ω res is the resonance frequency of the oscillator. When we use a source a = F/ω 2 res , homogeneous to a displacement, equation (3.1) is replaced by: ẍ + 2λ ẋ + ω 2 res x = aω 2 res cos(ω 0 t) -αx 2 -βx 3 .

(3.2)

Linear case

In the linear case, i.e. α = 0 and β = 0, the relation between amplitude of the forced oscillation b and the external excitation (F, ω 0 ) or (a, ω 0 ), near the resonance, is:

b 2 ξ 2 + λ 2 = F 2 4ω 2 res = a 2 ω 2 res 4 , ( 3.3)
where ξ = ω 0 -ω res .

Nonlinear case

For such a resonator, it is well known [START_REF] Landau | Physique théorique: Mécanique. Physique théorique[END_REF][START_REF] Nayfeh | Nonlinear Oscillations[END_REF] that the nonlinear behavior of the oscillations leads to a frequency shift of the resonance proportional to the square of the amplitude of vibration b:

ω ′ res = ω res + χb 2 , ( 3.4)
where:

χ = 3β 8ω res - 5α 2 12ω 3 res . (3.5)
It is easy to see on equation (3.5) that the nonlinearity factor χ can be positive or negative in the case of purely cubic nonlinearities, depending on the sign of β. However, in the case of purely quadratic nonlinearities, χ will always be negative, independently of the sign of α.

In the nonlinear case, we replace ξ, defined in the linear case, by ξ ′ :

ξ ′ = ω 0 -ω ′ res = ω 0 -ω res -χb 2 = ξ -χb 2 . (3.6)
In such way, the relation between amplitude of the forced oscillation b and the external excitation (a, ω), near the resonance, can be written as:

b 2 a 2 = ω 2 res 4 (ξ ′2 + λ 2 ) = ω 2 res 4 (ξ -χb 2 ) 2 + λ 2 . ( 3.7)
If the damping coefficient λ is low enough, the equation ξ(b) for the nonlinear resonance is then:

ξ = χb 2 ± aω res 2b 2 -λ 2 .
(3.8) For χ = 0 (linear case, Figure 3.1a), the resonance curve has the shape of a Lorentzian function centered at ξ = 0. For χ 0, for instance χ = 0.2 s 2 m -1 , the behavior of the resonance curve depends on the value of the excitation a.

Study of the Resonance Curves

For values of the excitation around a = 5 m (Figure 3.1b), the shape of the curve is similar to that in the

Position of points C and D

To give further insight of the behavior of the considered nonlinear oscillator, analytical expressions of the characteristic points C and D of the nonlinear resonance curve, shown on Figure 3.1d, will now be derived.

We define C and D as the two points where db/dξ → ∞. The derivation of equation (3.7) gives:

b 2 2 ξ -χb 2 1 -χ2b db dξ + 2b db dξ ξ -χb 2 2 + λ 2 = 0. (3.9)
Regrouping the terms depending on db/dξ gives:

2b 2 ξ -χb 2 + db dξ -4 ξ -χb 2 χb 3 + 2b ξ -χb 2 2 + λ 2 = 0, (3.10)
or:

db dξ -4ξχb 3 + 4χ 2 b 5 + 2bξ 2 + 2χ 2 b 5 -4ξχb 3 + 2bλ 2 = -2b ξb -χb 3 . (3.11)
Rearranging this equation leads to the following expression of db/dξ:

db dξ = -(ξb -χb 3 ) 3χ 2 b 4 + ξ 2 -4ξχb 2 + λ 2 .
(3.12)

As we are interested in the points where db/dξ → ∞, we obtain the following equation:

3χ 2 b 4 + ξ 2 -4ξχb 2 + λ 2 = 0. (3.13)
Equation (3.13) is the characteristic equation for the position of points C and D. We suppose that the damping coefficient λ is low enough so that χ 2 b 4 -λ 2 > 0. The solution of this equation is then given by:

ξ = 2χb 2 ± χ 2 b 4 -λ 2 . (3.14)

Maximum value

The maximum value of the amplitude is reached for db/dξ = 0.

According to equation (3.12), when db/dξ = 0, this is equivalent to:

bξ + χb 3 = 0 or ξ = χb 2 .

(3.15)

In the other hand, according to equation (3.8), we have:

ξ = χb 2 ± aω res 2b 2 -λ 2 . (3.16)
Combining those two equations, we obtain:

b 2 max λ 2 = a 2 ω 2 res 4 , ( 3.17)
leading finally to:

b max = aω res 2λ = Qa with Q = ω res 2λ . (3.18)
So, independently of the nonlinearities, the maximum amplitude is always the quality factor Q times the excitation amplitude a.

Critical value

We define a k , the critical value of a which corresponds to the case where points C and D are the same. This case is shown on Figure 3.1c. This time, we resolve equation (3.13) according to b 2 . The discriminant of this equation is:

∆ ′ = 4 χ 2 ξ 2 -3χ 2 λ 2 . (3.19)
Putting ∆ ′ = 0, we get:

ξ 2 = 3λ 2 , ( 3.20)
which leads to the following solution:

b 2 = 4χξ ± √ ∆ ′ 6χ 2 = 2ξ 3χ . (3.21)
Injecting this solution in equation (3.7), we obtain:

2 √ 3λ 3χ        √ 3λ - 2 √ 3λ 3 2 + λ 2        = a 2 k ω 2 res 4 , (3.22)
giving finally:

a 2 k = 4ω res 3 √ 3Q 3 χ . ( 3.23)
The system is bistable for a > a k . The threshold a k can be decreased by increasing the nonlinearity factor χ or the quality factor Q.

According to the expression of χ given in equation (3.5), we obtain, for the quadratic-only nonlinear resonator:

a 2 k = 48ω 4 res 15 √ 3Q 3 α 2 , ( 3.24)
and for the cubic-only nonlinear resonator:

a 2 k = 32ω 2 res 9 √ 3Q 3 β . (3.25)
The quadratic case always corresponds to a softening of the material (χ < 0). The cubic case corresponds to a softening or a hardening of the material, depending on the sign of β.

Analytical Model

After those considerations about nonlinear oscillators, we introduce three preliminary calculations which will be useful for our study of a nonlinear acoustic resonator.

Derivation of the Resonator Equation

In this section, we write the Lagrange equation for an elastic resonator, as in [START_REF] Preobrazhenskii | Nonlinear self-action of sound waves in an antiferromagnet with easy-plane anisotropy[END_REF], and derive it in order to obtain the expression of the Duffing equation describing the resonance behavior of the system, and link its parameters with the physical and geometrical parameters of this system.

We expand elastic deformation into normal vibration modes, writing the displacement as:

u(z, t) = n A n (t)U n (z). (3.26)
The Lagrangian of the elastic resonator is [START_REF] Klimov | Ferromagnetic resonance and magnetoelastic demodulation in thin active films with an uniaxial anisotropy[END_REF]:

L = 1 2 n m n Ȧ2 n -U L -U N L , ( 3.27)
where

m n = S ρU 2 n dS. (3.28)
The linear U L n and nonlinear U N L n part of the elastic energy are given by:

U L n = n 1 2 µ n A 2 n , (3.29)
U N L n = mnq 1 3! Φ mnq A n A m A q + mnqq ′ 1 4! Ψ mnqq ′ A n A m A q A ′ q , ( 3.30)
with

µ n = Ĉ2 ε n (z) 2 dr, ( 3.31)
Φ mnq = Ĉ3 ε m (z)ε n (z)ε q (z)dr, (3.32)
Ψ mnqq ′ = Ĉ4 ε m (z)ε n (z)ε q (z)ε ′ q (z)dr. (3.33)
In nonlinear elasticity, the strain ε is linked, to second order in 1D, to the displacement by:

ε(z) = ∂u ∂z + 1 2 ∂u ∂z 2 . (3.34)
Considering only the mode n, not interacting with other modes of the resonator, the equation of motion for A n (t) can be derived from Lagrange's differential equation:

∂ ∂t ∂L n ∂ Ȧn - ∂L n ∂A n + ∂D ∂ Ȧn = 0, (3.35)
where the expression of the Rayleigh dissipation is given by:

D = δ n m n Ȧ2 n . (3.36)
Considering in a first approximation (i.e. neglecting the geometrical nonlinearities term) that the strain is:

ε(z, t) ≈ ∂u ∂z = n A n (t) ∂U n ∂z , ( 3.37)
we rewrite the Lagrangian of the mode n as:

L n = 1 2 m n Ȧ2 n -U L n -U N L n , (3.38)
with

U L n = Ĉ2 2 A 2 n (t) z ∂U n (z) ∂z 2 dz, ( 3.39)
and

U N L n = Ĉ3 3! A 3 n (t) z ∂U n (z) ∂z 3 dz + Ĉ4 4! A 4 n (t) z ∂U n (z) ∂z 4 dz. (3.40)
Using equation (3.35), we obtain the following equation of the nonlinear resonator:

Än + 2δ n Ȧn + A n ω 2 res = -Ĉ3 2ρ z ∂U n (z) ∂z 3 dz z U 2 n dz A 2 n -Ĉ4 6ρ z ∂U n (z) ∂z 4 dz z U 2 n dz A 3 n , (3.41)
with ω res the resonance frequency. We recognize the expression of the Duffing resonator [START_REF] Landau | Physique théorique: Mécanique. Physique théorique[END_REF][START_REF] Nayfeh | Nonlinear Oscillations[END_REF] considered in section 3.3.1:

ẍ + 2λ ẋ + ω 2 res x = -αx 2 -βx 3 , (3.42)
where x is the displacement and:

ω 2 res = Ĉ2 ρ z ∂U n (z) ∂z 2 dz z U 2 n dz , (3.43) α = Ĉ3 2ρ z ∂U n (z) ∂z 3 dz z U 2 n dz , ( 3.44
)

β = Ĉ4 6ρ z ∂U n (z) ∂z 4 dz z U 2 n dz . (3.45)
In the specific case considered in this thesis, of a purely longitudinal elastic wave, we have: We now study the transmission through the PnC composed of N bi-layers and a central defect, similar to the one shown on Figure 3.2. The incident power P in is transferred from the material on the left to the defect by an evanescent wave. Let the elastic power be defined locally as

ω 2 res = C 11 ρ z ∂U n (z) ∂z 2 dz z U 2 n dz . ( 3

Complete Model of the Transmission

Z 2 Z 1 Z 2 Z 1 Z 2 Z 1 Z 2 Z 1 Z 2 Z 1 Z 2 a → a ′ → b → b ′ P in T 1 Q T 2 P out
P = Za 2 ω 2 /2 = Zv 2 /2 (in W .m -2
), where Z is the local impedance, a is the local displacement, v is the particle velocity and ω is the angular excitation frequency. Writing T 1 the transmission factor in power through the section of the crystal on the left of the defect, the elastic power P ′ in entering the defect is given by:

P ′ in = T 1 P in = Z 1 2 a ′ 2 ω 2 , (3.51)
with a ′ the displacement amplitude that enters the defect, ω the angular excitation frequency of the wave, and Z 1 the impedance of the defect.

In the nonlinear case, according to equation (3.7), the wave is amplified in the defect, and its amplitude becomes

b: b 2 a ′ 2 = ω 2 res 4 (ξ -χb 2 ) 2 + 4λ 2 = 1 2(ω-ω res ) ω res -2χb 2 ω res 2 + 1 Q 2 = Q 2 1 + ∆ + 2χQb 2 ω res 2 , (3.52) with ∆ = - ξ λ = 2Q (ω res -ω) ω res . (3.53)
Here, Q is the quality factor of the resonator, and ω res = 2πf res . As for the incident power, the power inside the defect is transferred to the material on the right by an evanescent wave. Writing T 2 as the transmission factor in power through the section of the crystal on the right of the defect, the output power P out is given by:

P out = Z 2 2 b ′2 ω 2 = Z 1 2 b 2 T 2 ω 2 , (3.54)
where b ′ is the displacement amplitude in the matrix on the right side of the crystal and Z 2 is the impedance of the matrix. Using equation (3.52), we can write:

P out P in = Z 1 T 1 T 2 b 2 Z 1 a ′ 2 = T 1 T 2 Q 2 1 + ∆ + 4χQP out ω res Z 1 T 2 ω 2 2 .
(3.55)

In our case, the PnCs on each side of the cavity are symmetrical, so, T 1 = T 2 . Also, the resonator has an amplification factor Q (its quality factor) for the input wave in displacement (or particle velocity), which is an amplification factor Q 2 in energy (or power). At the resonance, full transmission is reached, so the amplification factor in power compensates exactly the attenuation factors, so T 1 T 2 Q 2 = 1 which gives us:

T 1 = T 2 = 1 Q . (3.56)
Using this last expression, the transmission coefficient can be written as:

P out P in = 1 1 + P out P 0 -∆ 2 , ( 3.57)
with

P 0 = - Z 1 ω res ω 2 4χQ 2 . (3.58)
As we will see now, P 0 is a reference power, which characterizes the behavior of the system. We have defined P 0 with a minus sign, so that, as χ < 0 (softening of the material), P 0 > 0. We would define it with a plus sign if we had χ > 0 (hardening of the material).

Study of the "S" Curve

Equation (3.57) gives an implicit relation between P in and P out , which we study here. An example of such curve (that we will call the "S" curve due to its fold) is shown on Figure 3.3. Point A is the point described in equation (3.59), which corresponds to a full transmission (P in = P out ).

Point A corresponds to a full transmission of energy through the resonator, which, according to equation (3.57), is reached when:

P out P in = 1 ⇐⇒ P out = ∆P 0 = P in , ( 3.59)
so:

A (∆P 0 , ∆P 0 ) . (3.60)
Points B and C correspond to a critical values where the system will "jump" from an unstable to a stable state. We find their coordinates by differentiating equation (3.57) written as:

P in = P out P 0 -∆ 2 P out + P out . (3.61)
The differentiation gives:

dP in dP out = 3 P 2 out P 2 0 -4∆ P out P 0 + 1 + ∆ 2 . (3.62)
Looking at the point where dP in /dP out = 0, we obtain:

dP in dP out = 0 ⇐⇒ P out P 0 = 2∆ 3 ± ∆ 2 9 - 1 3 . (3.63)
Finally, their coordinates are:

B        P 0 27 2∆(9 + ∆ 2 ) -2(∆ 2 -3) 3/2 , 2∆P 0 3 + P 0 ∆ 2 9 - 1 3        , C        P 0 27 2∆(9 + ∆ 2 ) + 2(∆ 2 -3) 3/2 , 2∆P 0 3 -P 0 ∆ 2 9 - 1 3        .
(3.64)

The width of the "S" curve is:

P C in -P B in = 4P 0 27 (∆ 2 -3) 3/2 .
(3.65)

Models and Methods

The initial structure consists in a 1D PnC constituted of five steel slabs inserted inside a SiO 2 (silica) matrix.

We consider here only longitudinal elastic waves. The physical characteristics (density ρ, linear elastic modulus C 11 and longitudinal elastic wave velocity c L ) of steel and silica are given in Table 3. We perform a time resolution of the equations of elastodynamics, using a pseudo-spectral method described in chapter 1 [START_REF] Bou Matar | Two-dimensional axisymmetric numerical simulation of supercritical phase conjugation of ultrasound in active solid media[END_REF]. The reason for performing a time domain resolution is that the nonlinear effects cannot be taken into account by a "static" method such as the transfer matrix method as described in [START_REF] Markoš | Wave Propagation: From Electrons to Photonic Crystals and Left-Handed Materials[END_REF].

The pseudo-spectral method is appropriate for nonlinear simulations as it allows a high-order resolution [START_REF] Fornberg | A practical guide to pseudospectral methods[END_REF], which is useful in our case where the second and third harmonics will play an important role. With this method, simulating on a space grid of n points is equivalent to using a FDTD resolution at order n (n = 256 in our case).

The following equations are implemented in the pseudo-spectral method:

∂v ∂t = 1 ρ ∂τ ∂z , ( 3.66)
and

∂τ ∂t = k(t) ∂v ∂z , ( 3.67)
where v is the particle velocity, and τ is the stress. In equation (3.67), k(t) is defined as [START_REF] Naugolnykh | Nonlinear Wave Processes in Acoustics[END_REF]:

k(t) = ∂τ ∂ε = C 11 1 -2Γε -3δε 2 , ( 3.68)
where ε = ∂u/∂z is the strain and Γ and δ are dimensionless quantities involving the nonlinear characteristics of the PnC components. Γ and δ are related to the longitudinal second C 11 , third C 111 , and fourth C 1111 elastic constants by:

Γ = - 3 2 + C 111 2C 11 , ( 3.69)
δ = - 1 2 + C 111 C 11 + C 1111 6C 11 . ( 3.70)
Here, to simplify the calculation, only the defect is made of a nonlinear medium. All the materials of the other layers are considered as linear materials, i.e. the nonlinear constants Γ and δ of these materials vanish.

Our simulations with all the steel inclusions being nonlinear or only the defect show that this does not change the observations, as the nonlinear interaction between the wave and the medium is negligible where no resonance happens. The quadratic and cubic nonlinear elastic constants C 111 and C 1111 (according to [START_REF] Yang | Effect of aging on the third-order elastic moduli of 18N i maraging steel[END_REF]), and the corresponding dimensionless quantities Γ and δ of the defect are given in table The incident longitudinal wave is generated by a particle velocity, imposed at the position z = 4.6 µm and the output particle velocity is measured at z = 73.9 µm. Two Nearly-Perfectly Matched Layers (NPML), with thickness 3.0665 µm, defined as explained in [START_REF] Bou Matar | A nodal discontinuous Galerkin finite element method for nonlinear elastic wave propagation[END_REF] are used, one placed at each end of the structure.

Results

In this section, the results obtained using the pseudo-spectral method are detailed. First, a time resolution is performed using a large bandwidth source and without taking the nonlinearities into account, which allows to draw the transmission spectrum of the structure. Then, another time resolution at a specific excitation frequency is performed to observe an hysteretic behavior due to the nonlinear effects.

Linear Resonator

The transmission spectra for longitudinal waves at normal incidence of the two linear PnCs (perfect and defected) were calculated between 0 and 1 GHz with the transfer matrix method as described in chapter 1 [START_REF] Markoš | Wave Propagation: From Electrons to Photonic Crystals and Left-Handed Materials[END_REF]. The spectrum of the perfect crystal is represented on Figure 3.5. For this perfect structure, a wide band gap is observed between 300 MHz and 700 MHz. The spectrum of the defected structure, calculated for a mean of validation, confirms the results obtained with the pseudo-spectral method. The transmission spectra for the defected structure obtained with the transfer matrix method and the pseudo-spectral method are not distinguishable.

The simulation is first performed using an excitation consisting in a brief (large bandwidth) Ricker impulsion [START_REF] Bou Matar | A nodal discontinuous Galerkin finite element method for nonlinear elastic wave propagation[END_REF]. The transmission spectrum is calculated and normalized with the transmission of an homogeneous SiO 2 bar. This way, the resonance peak of the exact simulated structure is obtained and represented on the insert of Figure 3.5. By fitting it with a Lorentzian curve, its characteristics are obtained: its central frequency is f res = 491.95 MHz and its quality factor, defined as the ratio between f res and its full width at half maximum (FWHM), is Q = 608.36. The inset of Figure 3.5 presents a zoom of this peak. The central location of the defect in the structure ensures a maximal amplitude (close to 1) of the defect peak.

Nonlinear Resonator

The source is a sine function with frequency f c = 490.1 MHz, which is below the resonance frequency f res = 491.95 MHz of the nonlinear cavity. This sine function imposes a particle velocity that is varied incrementally in 60 levels, from 0.76 to 45.7 m s -1 . Each level is maintained for 4 millions dt, corresponding to 10 µs, since the timestep is dt = 2.5 ps. Calculations run first through the 60 levels in ascending order and then in descending order.

The particle velocity is recorded at the receiver after a period of 3 millions dt for each level. An hysteretic curve, represented on Figure 3.6, is observed. For both ascending and descending variations, the values reached by the output particle velocity are not proportional to the input particle velocity: for both curves, a positive step is observed, but the positions of these steps are different. The step for the increasing levels corresponds to a higher input particle velocity than the step corresponding to the decreasing levels.

As a short hand-waving interpretation, we anticipate that due to the nonlinearities in the resonating cavity, a frequency shift appears as the particle velocity increases. Due to the softening of the nonlinear material with increasing strain, the resonance frequency of the cavity, initially f res , decreases and reaches the excitation frequency f c . So, a resonance happens in the cavity and the elastic wave is transmitted through the crystal. Figure 3.7 shows such a frequency shift, for a structure with 7 inclusions.

Discussion

This section provides an interpretation of the results described above. For this purpose, an analytical model is proposed and a parametric study of its parameters is provided.

The complete model of the transmission obtained in section 3.4 allows us to draw a theoretical curve of the output elastic power versus the input elastic power and compare it with the results of the numerical simulations. This theoretical curve has the shape of an "S" and will therefore be referred as the "S" curve in the following parametric study. Figure 3.8 shows a reasonably good agreement between the hysteresis curve derived from the numerical experiment and the "S" curve. On this figure, we can see that the system possesses two stable states. The first one, corresponding to low particle velocities, is displayed on the bottom left corner of the figure. The second one, corresponding to high particle velocities, is displayed on the top right corner of the figure. Moreover, in a range of powers of P in the system becomes bi-stable. In this range, the state chosen by the system depends on its previous state: high output powers when high input powers have already been reached and low output powers in the other case. The portion of the "S" curve with a negative slope corresponds to unstable states Figure 3.6 -Output particle velocity versus the input particle velocity. The arrows show the paths followed by the output particle velocity for increasing and decreasing source levels. and therefore cannot be reached physically.

Parametric Study

The previous analytical result allows one to predict the characteristics of the hysteretic behavior of the system for a wide range of variation of the parameters. A parametric study of this "S" curve will now be performed to show how the shape of the hysteresis depends on the resonance parameters, namely the quality factor Q, the nonlinear parameter β and the dimensionless quantity, ∆.

Then, the dependence of the quality factor Q on the number of bilayers of the system and the impedance break of those bilayers will be studied. Finally, the respective effects of the quadratic and cubic nonlinearities will be compared.

Study of the "S" curve

The theoretical curve, which analytical formula is given in equation (3.57), has been plotted on Figure 3.9 for different varying parameters, all other parameters being maintained constant, with values close to those used in Figure 3.8. On these figures, the straight line corresponds to the full transmission (P in = P out ). As one Figure 3.7 -Resonance frequency of the defected structure for different input particle velocities. A frequency shift to lower frequencies is clearly visible when the input particle velocity increases. This study was made for a structure with 7 inclusions. can see, the full transmission is only reached for one point of each "S" curve and this point is located just before the step on the path followed when the input power decreases. Figure 3.9a shows the variation of the shape of the "S" curve when β is varied from -4 × 10 31 m -2 s -2 to -7.5 × 10 31 m -2 s -2 while ∆ = 4.51 and Q = 600 are kept constant. As one can see, the curve appears contracted in both x-and y-directions when the absolute value of β is increased . Thus, the shape itself is not modified, only the magnification factor. Notably, the negative slope in the middle of the curve remains the same. This behavior indicates that when the nonlinearity of the material increases, the step will be reached at a lower excitation power and the hysteresis curve will reach a saturation value that is lower than that for smaller nonlinearities.

On Figure 3.9b, the value of ∆ is varied from 1.34 to 4.76 while β = -4.5 × 10 31 m -2 s -2 and Q = 600 are kept constant. As previously a dilatation of the curve is observed in both x-and y-directions, but an enlargement of the curve also occurs. One can interpret this as the fact that when ∆ (i.e. the difference between the resonance frequency and the excitation frequency) is increased, the width of the hysteresis is increased as well as the required power to reach the step.

On Figure 3.9c, the value of Q is varied from 500 to 800 while ∆ = 4.51 and β = -4.5 × 10 31 m -2 s -2 are kept constant. This variation leads to a dilatation of the curve in the x-direction and to a shrinkage of the curve along the y-direction. Then, when the quality factor Q of the cavity is increased, a higher input power is required to reach the step, and the height of the step is decreased. Two important parameters of the "S" curve, namely the threshold power P C in and the width of the hysteresis P C in -P B in are now studied. The evolution of P C in and P C in -P B in versus β, ∆ and Q are plotted on Figure 3.10. Figure 3.10a shows the variations of P C in and P C in -P B in as a function of β when ∆ and Q are kept fixed. The decreasing behavior of these quantities with increasing β (in absolute value) indicates that a higher nonlinear constant will result in a lower threshold power and a smaller hysteresis.

In Figure 3.10b the variations of P C in and P C in -P B in versus ∆ when β and Q are kept fixed are reported and show that both P C in and P C in -P B in increase with ∆. This means that a larger difference between the excitation frequency and the resonance frequency will result in a higher threshold power and a larger hysteresis.

Finally Figure 3.10c, shows that both P C in increases proportionally and P C in -P B in increases almost proportionally with Q when ∆ and β are fixed. This means that a resonating cavity with a higher quality factor will require a larger energy to resonate, and that its hysteresis will be wider.

We observed that when the nonlinear parameter δ is increased by a factor X, then, the maximum value for the strain reached in the default is decreased by a factor √ X. We observed that when the quality factor is increased by a factor X, then, the maximum value for the strain reached in the default is decreased by a factor √ X. Another effect of increasing the quality factor is an increase of the number of wave periods needed for exciting the resonator. The order of magnitude of the time that is necessary to observe a resonance in a cavity with quality factor Q is 5QT where T is the wave period. In our case, T /dt ≈ 800, so we need about 2.4 million timesteps to reach the resonance.

Finally, we tried to include quadratic nonlinearities in our simulations. In this case, we couldn't observe a step in the output particle velocity, unless we increased the cubic nonlinearities. So, we observed that quadratic nonlinearities seem to break the hysteretic phenomenon, but, at the time of writing this manuscript, we cannot propose a satisfactory explanation. We can only make the hypothesis that the second harmonic will escape the resonating cavity and, consequently, the quality factor of the defect will decrease and the behavior of the system will move away from the theoretical model.

Conclusion

A one-dimensional PnC composed of steel inclusions in a SiO 2 matrix, with a nonlinear, half-wavelength central inclusion has been designed. Numerical simulations of the transmission of an elastic wave through the crystal have been performed. It has been shown that the central nonlinear inclusion constitutes a resonating cavity which is able to commute the elastic wave depending on the particle velocity. An analytical model of this phenomenon has been provided, based on the derivation of the Lagrange equation. A parametric study has been performed for the parameters of the system, namely the nonlinear elastic constants, the quality factor of the cavity and the difference between the excitation frequency and resonance frequency of the cavity. The relations between those parameters and the behavior of the system have been explored and the hysteresis phenomena has been explained with a complete model of the transmission through the bilayers.

This realization is considered to be the first step to the simulation of tunable PnCs. Because of the considerable dependence of the hysteresis phenomena on the nonlinear elastic constants, one being able to dynamically tailor these constants would create a controllable switch. So, the dynamic variation of the nonlinear elastic constants induced by the magneto-elastic effect will be analyzed in a future study, aiming the creation of magnetically controlled PnC commutators, with tailored properties [START_REF] Bou Matar | Band gap tunability of magneto-elastic phononic crystal[END_REF].

Also, this result in a one-dimensional structure will lead to the design and modeling of a two-dimensional elastic resonator. This structure will be simulated numerically using efficient tools that have been developed

for numerical simulation of nonlinear elastic wave propagation.

Chapter 4

Mixing Laws for the Quadratic and Cubic Elastic Constants

Introduction

In chapter 3, we have designed a nonlinear elastic resonator. We have performed numerical studies and explained their results with a complete analytical model taking into account second-and third-order nonlinearities. We have shown that the switching behavior of this resonator relies on cubic nonlinearities and can be improved by increasing those cubic nonlinearities and decreasing quadratic nonlinearities.

In this chapter, we will derive the mixing laws for the quadratic and cubic elastic constants. For nonlinear inclusions of a material inside a linear matrix, we will calculate the values of the linear as well as quadratic and cubic nonlinear constants of the effective material, through an homogenization procedure. We will demonstrate that both quadratic and cubic nonlinearities can be amplified in the effective material, for specific concentrations of the nonlinear inclusions.

Even though the inclusions would have "reasonable" nonlinear parameters, we hope we will obtain an effective material with very high cubic and low quadratic elastic nonlinearities. This effective material could then compose the defect of the resonator of chapter 3.

In a first section, we will derive the equations of a one-dimensional mass-spring system to show that it is actually possible to amplify both quadratic and cubic nonlinearities. Then, we will explain the 3D fluid model and calculate the 2D and 3D mixing laws for the second-and third-order elastic parameters.

We will finally provide equations for the prediction of the optimal concentrations and the corresponding amplification of quadratic and cubic parameters, according to the derivation order and the elastic constants of the constituent materials.

1D Mass-Spring System Model

This section is adapted from [START_REF] Zaitsev | A model of anomalous elastic nonlinearity of microinhomogeneous media[END_REF]. We study the mass-spring system model presented on Figure 4.1. This system is composed of N = N 1 + N 2 springs among which N 1 have an elastic coefficient K 1 and N 2 have an elastic coefficient K 2 . The initial length of the springs is L and L ≪ λ where λ is the wavelength of the excitation. We write the expression of the stress σ, for the material of the rigid (K 1 ) and soft (K 2 ) springs, respectively:

σ = M 1 ε (1 + f 1 (ε)) and σ = M 2 ε (1 + f 2 (ε)) , (4.1)
where f 1 (ε) and f 2 (ε) are nonlinear functions of ε which verify: f 1 (ε) ≈ f 2 (ε) ≪ 1, and M 1 = K 1 /L and M 2 = K 2 /L. M 1 and M 2 are the springs constants per unit length of the springs.

We consider the case of "soft" inclusions:

K 1 K 2 = Q = M 1 M 2 ≫ 1. (4.2)
In the quasi-static case (λ ≫ L), we have:

σ = K 1 X 1 = K 2 X 2
where X 1 and X 2 are the elongations of the springs with nonlinear constants K 1 and K 2 respectively.

The total elongation of a chain with length N × L is:

X t = X 1 N 1 + X 2 N 2 . (4.3)
The strain in that chain is given by:

ε = X t N L = X 1 L N 1 N + X 2 L N 2 N = X 1 L (1 -c) + X 2 L c = X 1 L (1 -c + Qc) , (4.4)
since X 2 = QX 1 . From now, c = N 2 /N is the concentration of soft springs.

The average energy density W is linked to the elastic energy stored in the soft and hard springs by:

W = N 1 N W h + N 2 N W s , ( 4.5)
where:

W h = X 1 0 M 1 L ξ L 1 + f 1 ξ L dξ, and W s = X 2 0 M 2 L ξ L 1 + f 2 ξ L dξ. (4.6)
So, the average energy density can be rewritten as:

W = N 1 N L X 1 0 M 1 ξ L 1 + f 1 ξ L dξ + N 2 N L X 2 0 M 2 ξ L 1 + f 2 ξ L dξ. (4.7)
And from equation (4.4), we get:

X 1 = Lε 1 + c(Q -1)
and

X 2 = QX 1 . (4.8)
The stress is then given by:

σ = ∂W ∂ε = ∂W ∂X 1 ∂X 1 ∂ε where ∂X 1 ∂ε = L 1 + c(Q -1) , ( 4.9)
and

∂W ∂X 1 = (1 -c) M 1 L X 1 L 1 + f 1 X 1 L + c M 2 L QX 1 L 1 + f 2 QX 1 L ∂X 2 ∂X 1 . (4.10)
We define:

α = X 1 L = ε 1 + c(Q -1) , (4.11)
and we use QM 2 = M 1 and ∂X 2 /∂X 1 = Q in equation (4.10), leading to:

∂W ∂X 1 = (1 -c) M 1 L α (1 + f 1 (α)) + c M 1 L Qα (1 + f 2 (Qα)) = M 1 L α (1 + c(Q -1) + (1 -c)f 1 (α) + cQf 2 (Qα)) .
(4.12)

Finally, we obtain the following stress-strain relation for the inhomogeneous string of springs:

σ = ∂W ∂ε = M 1 1 + c(Q -1) × α (1 + c(Q -1) + (1 -c)f 1 (α) + cQf 2 (Qα)) = M 1 α 1 + 1 -c 1 + c(Q -1) f 1 (α) + cQ 1 + c(Q -1) f 2 (Qα) , (4.13)
which corresponds to equation ( 5) in [START_REF] Belyaeva | Nonlinear elastic properties of microinhomogeneous hierarchically structured media[END_REF].

Low-concentration Case

In the low-concentration case, i.e. when c ≪ 1/Q, then:

α = ε 1 + c(Q -1)
≈ ε so Qα ≈ Qε. (4.14) In this case, the inclusions undergo a strain that is Q ≫ 1 times bigger than the average strain the medium undergoes. This fact explains why the nonlinear parameters of the heterogeneous medium can become quite large.

High-concentration Case

We now define the nonlinear functions f 1 and f 2 with a quadratic term Γ i and a cubic term

δ i : f 1 (ε) = Γ 1 ε+δ 1 ε 2 and f 2 (ε) = Γ 2 ε + δ 2 ε 2 .
We replace f 1 and f 2 by their expressions in equation ( 4.13) and we obtain:

σ = M 1 α 1 + 1 -c 1 + c(Q -1) Γ 1 α + δ 1 α 2 + cQ 1 + c(Q -1) Γ 2 Qα + δ 2 (Qα) 2 = M 1 ε 1 + c(Q -1) 1 + (1 -c)Γ 1 ε (1 + c(Q -1)) 2 + (1 -c)δ 1 ε 2 (1 + c(Q -1)) 3 + cΓ 2 Q 2 ε (1 + c(Q -1)) 2 + cδ 2 Q 3 ε 2 (1 + c(Q -1)) 3 = M 1 ε 1 + c(Q -1)         1 + Γ 1 ε 1 -c + c Γ 2 Γ 1 Q 2 (1 + c(Q -1)) 2 + δ 1 ε 2 1 -c + c δ 2 δ 1 Q 3 (1 + c(Q -1)) 3         . (4.15)
Finally, we rewrite the stress in the following form:

σ = M 1 A K ε 1 + Γ 1 εB K + δ 1 ε 2 C K , ( 4

.16)

with:

A K = 1 1 + c(Q -1)
and

B K = 1 -c + c Γ 2 Γ 1 Q 2 (1 + c(Q -1)) 2 and C K = 1 -c + c δ 2 δ 1 Q 3 (1 + c(Q -1)) 3 . (4.17)
In the specific case where Γ 1 = Γ 2 , the expressions of A K and B K correspond to equations ( 7) and ( 8) in [START_REF] Belyaeva | Nonlinear elastic properties of microinhomogeneous hierarchically structured media[END_REF].

We will now solve for the optimal concentration that maximizes the nonlinear parameters of order n.

Here, X n represents B K (n = 2) or C K (n = 3), x i represents Γ i (n = 2) or δ i (n = 3). Taking the derivative of the nonlinear parameter X n by the concentration c, we obtain:

∂X n ∂c = ∂ ∂c       1 -c + c x 2 x 1 Q n (1 + c(Q -1)) n       = x 2 x 1 Q n -1 (1 + c(Q -1)) n + 1 + c x 2 x 1 Q n -1 (-n(Q -1)) (1 + c(Q -1)) n-1 , (4.18)
which is null if and only if:

x 2 x 1 Q n -1 (1 + c(Q -1)) -n(Q -1) 1 + c x 2 x 1 Q n -1 = 0. (4.19)
The previous expression gives us the optimal concentration c opt (n) as:

c opt (n) = 1 (n -1)(Q -1) - n (n -1) x 2 x 1 Q n -1 ≈ 1 (n -1)Q if Q ≫ 1. (4.20)
For n = 2 and n = 3, the optimal concentrations are given by:

c opt (n = 2) = 1 Q -1 - 2 Γ 2 Γ 1 Q 2 -1 ≈ 1 Q and c opt (n = 3) = 1 2(Q -1) - 3 2 δ 2 δ 1 Q 3 -1 ≈ 1 2Q . (4.21)
These expressions correspond to that obtained by the mixing law of the quadratic and cubic nonlinear coefficients of fluids, as it will be shown in the next section.

We can notice that the optimal concentrations for the quadratic and cubic nonlinearities are different: for n = 2 and n = 3, c opt ≈ 1/((n -1)Q). This allows us to amplify the cubic nonlinearity without changing too much the quadratic one. This will be useful for the nonlinear oscillator described in chapter 3.

When c = c opt (n), the linear elastic parameter A K becomes:

A K (c opt (n)) = 1 1 + Q-1 (n-1)Q = (n -1)Q (n -1)Q + Q -1 = (n -1)Q nQ -1 ≈ n -1 n , ( 4.22)
so, the linear parameter does not vary much.

On the other hand, the nonlinear parameters B K and C K , at c = c opt (n) ≈ 1/((n -1)Q), are:

B K opt , C K opt ≈ (n -1)Q -1 + x 2 x 1 Q n ((n -1)Q + Q -1) n × (n -1) n Q n (n -1)Q ≈ x 2 x 1 Q n (n -1) n-1 Q n-1 Q n n n ≈ x 2 x 1 × (n -1) n-1 Q n-1 n n (4.23)
At the optimal concentration for quadratic nonlinearities, the values of B K and C K are:

B K opt = B K (c opt (n = 2)) ≈ Γ 2 Γ 1 × Q 4 et C K (c opt (n = 2)) ≈ δ 2 δ 1 × Q 2 8 . (4.24)
At the optimal concentration for cubic nonlinearities, the values of B K and C K are:

B K (c opt (n = 3)) ≈ Γ 2 Γ 1 × 2Q 9 et C K opt = C K (c opt (n = 3)) ≈ δ 2 δ 1 × 4Q 2 27 . (4.25)
The important points to notice here are the ratios of amplification: K 1 /(4K 2 ) for n = 2 and 4K 2 1 /(27K 2 2 ) for n = 3, and the optimal concentrations: K 2 /K 1 for n = 2 and K 2 /(2K 1 ) for n = 3. In the next section, we will see that they are the same for fluids.

We notice that we can increase the cubic nonlinearity by a factor 100 by choosing:

4 27 Q 2 = 100 which is Q = K 1 K 2 ≈ 26 and c opt = 1 2Q = 1 52 (4.26)
In this case:

B K 1 2Q = 1 -1 2Q + 1 2Q Γ 2 Γ 1 Q 2 1 + 1 2Q (Q -1) 2 ≈ Γ 2 Γ 1 Q 2 1 + Q 2Q 2 = Γ 2 Γ 1 2Q 9 = 52 9 Γ 2 Γ 1 (4.27)
We notice that Γ 2 /Γ 1 is increased by a factor 5.8 when δ 2 /δ 1 is increased by a factor 100, which is 17.3 times bigger. This favors the cubic nonlinearities.

Mixing Law for a Fluid in 3D

For a fluid, we define the quadratic and cubic nonlinear parameters with the usual notations:

β = B l 2A l + 1 and γ = C l 6A l -1, (4.28)
where A l , B l , and C l are the first three coefficients of the Taylor's expansion of the pressure as a function of the density.

It has been demonstrated [START_REF] Everbach | A corrected mixture law for B/A[END_REF] that the mixing law for β and γ are:

1

K 2 m β m = N i=1 β i c i K 2 i and 1
K 3 m γ -β B l A l m = N i=1 γ -β B l A l i c i K 3 i , ( 4.29)
where N is the number of constituents of the mixture, K is the compressibility factor and c i is the volume fraction of each constituent N i=1 c i = 1 . The mixing law for K is [START_REF] Everbach | A corrected mixture law for B/A[END_REF]:

1 K m = N i=1 c i K i . (4.30)
For a mixture with two constituents, we use the notations c 2 = c and c 1 = 1 -c. If the second constituent is more compressible than the first one (K 2 < K 1 ), then:

β m = c K 2 m K 2 2 β 2 + (1 -c) K 2 m K 2 1 β 1 , (4.31)
and 1 K m = c K 2 + 1 -c K 1 so K m = K 1 K 2 cK 1 + (1 -c)K 2 . ( 4.32)
We solve for the concentration c that will optimize the nonlinear parameter β m of the mixture:

∂β m ∂c = β 2 K 2 2 - β 1 K 2 1 K 2 m + cβ 2 K 2 2 + (1 -c) β 1 K 2 1 ∂K 2 m ∂c , ( 4.33)
which is null if and only if:

cβ 2 K 2 2 + (1 -c) β 1 K 2 1 = - K 2 m ∂K 2 m ∂c β 2 K 2 2 - β 1 K 2 1 . ( 4.34)
We know that:

K 2 m = K 2 1 K 2 2 (cK 1 + (1 -c)K 2 ) 2 . ( 4.35)
So, taking the derivative of this expression as a function of c leads to:

∂K 2 m ∂c = -2K 2 1 K 2 2 (K 1 -K 2 ) (K 1 c + (1 -c) K 2 ) 3 . ( 4
K 2 m ∂K 2 m ∂c = K 2 1 K 2 2 (cK 1 + (1 -c)K 2 ) 2 × (K 1 c + (1 -c) K 2 ) 3 -2K 2 2 K 2 1 (K 1 -K 2 ) = K 1 c + (1 -c) K 1 2(K 1 -K 2 ) . ( 4.37)
Finally, introducing this latest result in equation (4.34) leads to the equation that we have to solve:

cβ 2 K 2 2 + (1 -c) β 1 K 2 1 = K 1 c + (1 -c) K 2 2(K 1 -K 2 )
β 2 K 2 2 - β 1 K 2 1 ⇐⇒ c = 2(K 2 -K 1 )
β 1 K 2 1 + K 2 β 2 K 2 2 - β 1 K 2 1 (K 1 -K 2 ) β 2 K 2 2 - β 1 K 2 1 . ( 4.38)
So, the concentration which maximizes the nonlinear quadratic (n = 2) parameter β m is:

c opt (n = 2) = 1 K 1 K 2 -1 - 2 β 2 β 1 K 2 1 K 2 2 -1 . ( 4.39)
The same calculation for optimizing the cubic (n = 3) parameter δ m = γ m -β m B m /A m gives:

c opt (n = 3) = 1 2 K 1 K 2 -1 - 3 2 δ 2 δ 1 K 3 1 K 3 2 -1 . ( 4.40)
If the second medium is much more compressible than the first one (K 2 ≪ K 1 ), then:

c opt (n = 2) ≈ K 2 K 1 -2 β 1 β 2 K 2 2 K 2 1 and c opt (n = 3) ≈ K 2 2K 1 - 3 2
β 1 β 2 K 3 2 K 3 1 . ( 4.41)
Now, if moreover, β 1 is of the same order as β 2 and δ 1 is of the same order as δ 2 , then:

c opt (n = 2) ≈ K 2 K 1 and c opt (n = 3) ≈ K 2 2K 1 which we write c opt (n) = K 2 (n -1)K 1 . ( 4.42)
Medium Compressibility factor (Pa) Quadratic nonlinear parameter Cubic nonlinear parameter

Water K 1 = 2.2 × 10 9 β 1 = 3.5 γ 1 = 5.25, δ 1 = -12.25 Air K 2 = 1.42 × 10 5 β 2 = 1.2 γ 2 = -1.04, δ 2 = -1.52
Table 4.1 -Compressibility factor and quadratic and cubic nonlinear constants of air and water.

The values of the nonlinear parameters that we obtain are:

β m c opt (n = 2) = K 1 K 2 β 2 4 + β 1 4 ≈ K 1 K 2 β 2 4 and δ m c opt (n = 3) = K 2 1 K 2 2 4δ 2 27 + 8δ 1 27 ≈ K 2 1 K 2 2 4δ 2 27 . ( 4.43)
As K 2 ≪ K 1 , the optimal concentrations can be very small. In such case, it appears that a small concentration of compressible inclusions leads to a huge (i.e. 4K 2 1 /(27K 2 2 )) increase of the global cubic nonlinearity of the medium, when the quadratic global nonlinearity is "only" increased by a factor K 1 /(4K 2 ).

We can observe that we obtain the same optimal concentrations c opt (n = 2) ≈ 1/Q = K 2 /K 1 and c opt (n = 3) ≈ 1/(2Q) = K 2 /(2K 1 ) are the same for the fluid and for the spring system of section 4.2. The amplification factors K 1 /(4K 2 ) and 4K 2 1 /(27K 2 2 ) are also the same. As an example, we can consider the case of air bubbles in water. The compressibility factors and nonlinear constants of air and water are given in table 4.1. For this system, the optimal concentration is:

c opt (n = 2) ≈ K 2 K 1 ≈ 6 × 10 -5 and c opt (n = 3) ≈ K 2 2K 1 ≈ 3 × 10 -5 , (4.44)
and the corresponding amplifications of the global nonlinearities in the medium are:

β m c opt (n = 2) ≈ K 1 K 2 β 2 4 ≈ 4650

Relations Between Constants for Isotropic Solids

Before extending the previous mixing laws for the nonlinear elastic constants in isotropic solids, we recall, in this section, the relations between:

• the Lamé constants: λ, µ,

• the compression modulus: K 2D in 2 dimensions or K 3D in 3 dimensions,

• the tensors C IJ for the linear and C IJK for the quadratic nonlinear elastic constants.

Compression Modulus

In 2D, we have:

K 2D = λ + µ = (λ + 2µ) -µ = ρ 0 c 2 l -c 2 t , (4.46)
with the following condition on the Poisson coefficient:

-1 < ν 2D = K 2D -µ 2K 2D < 1 2 , ( 4.47)
which corresponds to:

0 < µ < 3K 2D . ( 4.48)
In 3D, we have:

K 3D = λ + 2 3 µ = (λ + 2µ) - 4 3 µ = ρ 0 c 2 l - 4 3 c 2 t , (4.49)
with the following condition on the Poisson coefficient:

-1 < ν 3D = 3K 3D -2µ 2(3K 3D + µ) < 1 2 . (4.50)

Tensor for the Linear Constants

The C IJ tensor of an isotropic material is symmetric. Its components are linked to the Lamé constants by:

C 11 = C 22 = C 33 =λ + 2µ, C 12 = C 13 = C 23 =λ, C 44 = C 55 = C 66 =µ. (4.51)
Moreover, the components can be obtained from the wave speeds, using:

Longitudinal/plane: c l = c p = λ+2µ ρ 0 = C 11 ρ 0 , Transversal/shear: c t = c s = µ ρ 0 = C 44 ρ 0 , ( 4.52)

Notation by Notation by Notation by Notation by Murnaghan, Hughes and Kelly

Taupin and Bernstein Bland Eringen and Suhubi

l = B + C m = A 2 + B n = A A = n B = m -n 2 C = l -m + n 2 ν 1 = 2C ν 2 = B ν 3 = A 4 α = C 3 β = B γ = A 3 l E = A 3 + B + C 3 m E = -A -2B n E = A
µ = ρ 0 c 2 s , C 11 = λ + 2µ = ρ 0 c 2 p , C 12 = λ = ρ 0 (c 2 p -2c 2 s ). (4.53)

Tensor for the Nonlinear Quadratic Constants

The C IJK tensor of an isotropic material is symmetric. There are 20 non-null constants in the tensor. They can be given using three constants : A, B and C, called the Landau constants.

The quadratic elastic constants can be expressed with the three constants of Murnaghan, Hughes and Kelly Finally, the link between the Landau constants and the C IJK tensor is:

A = C 111 2 -3C 112 2 + C 123 , B = C 112 2 -C 123 2 , C = C 123 2 .
(4.54)

Units

Those constants are usually expressed in GPa or dyn/cm 2 . 1 GPa = 1 × 10 10 dyn/cm 2 .

Landau Coefficients in a Heterogeneous Medium

In the previous sections, we have observed that the 1D mass-spring system and the inclusion of air bubbles in water both possess quadratic and cubic nonlinear constants that are amplified compared to that of the inclusions. For both systems, we obtain the same expressions for the quadratic amplification and the same expressions for the cubic amplification.

The optimal concentrations are the same for both systems, and they are different for quadratic and for cubic nonlinearities. This will allow one to favor the cubic constants over the quadratic ones, by choosing a concentration equal to the optimal concentration for the cubic amplification.

In this section, we will extend the 1D model to 2D and 3D systems of nonlinear inclusions in a linear matrix. Our starting point is the model developed in [START_REF] Giordano | Nonlinear elastic Landau coefficients in heterogeneous materials[END_REF] that we rederive.

Derivation to the Second Order

We consider nonlinear inclusions of a material noted 2 in a linear matrix noted 1. The inclusions are ellipsoids.

In 3D

We write the elastic energy U stored inside a deformed material as a function of the strain tensor ε:

U = U( ε), (4.55)
and the stress can therefore be derived as follows:

T ij = ∂U ∂ε ij . (4.56)
If the material is isotropic, then the energy is an isotropic tensor and it can be written as [START_REF] Giordano | Nonlinear elasticity of composite materials[END_REF] :

U = µ Tr( ε2 ) + λ 2 Tr 2 ( ε) + A 3 Tr( ε3 ) + B Tr( ε) Tr( ε2 ) + C 3 Tr 3 ( ε), (4.57)
where, Tr is the trace operator. Tr( ε), Tr( ε2 ) and Tr( ε3 ) are called the "invariants" of the stress tensor ε.

We rewrite this expression up to the third order using the Einstein notation:

U = µε ij ε ji + λ 2 (ε ii ) 2 + A 3 ε ij ε jk ε ki + Bε ii ε kh ε hk + C 3 (ε ii ) 3 . ( 4.58)
We differentiate this expression to obtain T nm :

T nm = ∂U ∂ε nm = µδ in δ jm ε ji + µε ij δ jn δ im + 2λ 2 (ε jj )δ nm + A 3 δ in δ jm ε jk ε ki + A 3 ε ij δ jn δ km ε ki + A 3 ε ij ε jk δ kn δ im + Bδ nm ε kh ε hk + Bε ii (δ kn δ hm ε hk + ε kh δ hn δ km ) + 3C 3 (ε ii ) 2 δ nm . (4.59)
Or simply:

T nm = 2µε nm + λε ii δ nm + Aε mk ε kn + Bδ nm ε kh ε hk + 2Bε ii ε nm + C(ε ii ) 2 δ nm . ( 4.60)
So, finally, the tensor T is written:

T = 2µ ε + λ Tr( ε) Î + A ε2 + B Tr( ε2 ) Î + 2B ε Tr( ε) + C Tr 2 ( ε) Î . ( 4.61)
Here, λ and µ are the Lamé coefficients of the linear elastic theory. A, B and C are the second-order nonlinear moduli i.e. the Landau coefficients. We adopt the physical nonlinearity standpoint, with the geometrical nonlinearity neglected everywhere. The balance equations are based on the small-strain tensor and on the geometric Cauchy stress.

Restriction to 2D

In the 2D case, the tensor ε only has two invariants. We use the Cayley-Hamilton theorem and the following trick to express Tr( ε3 ) as a function of Tr( ε) and Tr( ε2 ). We define the characteristic polynomial of ε as:

P ε(λ) = det( ε -λ Î). (4.62)
According to the Cayley-Hamilton theorem, ε satisfies its own characteristic equation, so P ε( ε) = 0 ( 0 is the null matrix). Also, it is easy to check that for every 2 × 2 matrix A:

P A (λ) = det(A) -λ Tr(A) + λ 2 , ( 4.63)
so, with P ε applied to ε:

det( ε) Î -ε Tr( ε) + ε2 = 0. ( 4
ε2 = ε Tr( ε) - 1 2 Tr 2 ( ε) -Tr( ε2 ) Î , ( 4.66)
from which we obtain:

Tr( ε3 ) = 3 2 Tr( ε) Tr( ε2 ) - 1 2 Tr 3 ( ε). (4.67)
We replace ε2 and Tr( ε3 ) in equation ( 4.61) and we get:

T = 2µ ε + λ Tr( ε) Î + A ε Tr( ε) - 1 2 Tr 2 ( ε) -Tr( ε2 ) Î + B Tr(ε 2 ) Î + 2B ε Tr( ε) + C Tr 2 ( ε) Î . ( 4.68)
We regroup the terms of the previous expression and we define:

e = B + A 2 , 3f = C - A 2 .
(4.69)

So, the final expression for T in 2D is:

T nm = 2µε nm + λε ii δ nm + eε hk ε kh + 2eε ii ε nm + 3f (ε ii ) 2 δ nm , (4.70)
or written as a tensor:

T = 2µ ε + λ Tr( ε) Î + e Tr( ε2 ) Î + 2e ε Tr( ε) + 3f Tr 2 ( ε) Î . ( 4.71)

Comment

In fact, if a 3D system is studied without energy considerations, four coefficients, A, B, C and D, must be considered. In that case, the expression for the stress is:

T = 2µ ε + λ Tr( ε) Î + A ε2 + B Tr( ε2 ) Î + C Tr 2 ( ε) Î + D ε Tr( ε). (4.72)
If the thermodynamics of the system is considered, as in the Green formulation of the elasticity theory [START_REF] Giordano | Nonlinear elastic Landau coefficients in heterogeneous materials[END_REF][START_REF] Giordano | Nonlinear elasticity of composite materials[END_REF][START_REF] Landau | Theory of Elasticity[END_REF], those coefficients reduce to three independent numbers in 3D : A, B and C, as D = 2B and two independent numbers in 2D : e and f : A = 0, B = e, C = 3f and D = 2B = 2e.

In that case, we obtain the expressions (4.61) and (4.71) for the stress. For the study below, we will use the Green elasticity theory and its three independent coefficients.

Expressions for the Quadratic Nonlinearities

We write the bulk modulus as:

K 3D = λ + 2 3 µ, ( 4.73)
and the Poisson coefficient as:

ν = 3K 3D -2µ 2(3K 3D + µ) . ( 4.74)
The stress in a single ellipsoid inclusion is:

T = Ĉ2 ( ε) ε. (4.75)
and the stress in the matrix is:

T ∞ = Ĉ1 ε∞ , (4.76)
where Ĉ2 is the stiffness tensor of the inclusion and Ĉ1 the stiffness tensor of the matrix. The notation with ∞ means that the stress loading is located far from the inclusion.

If we write Ŝ the Eshelby tensor, which depends only on the geometry of the inclusion and on Ĉ1 , then the total stress in the ellipsoid inclusion is uniform and its value is [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF]83]:

εd = Î -Ŝ Î -Ĉ-1 1 Ĉ2 Ŝ -1 ε∞ . ( 4.77)
If this equation has a solution for a specific ε, called εd , then the nonlinear inclusion can be replaced by a linear one: Ĉ2 = Ĉ2 ( εd ) . (

The equation (4.77) is equivalent to: εd -Ŝ εd + Ŝ Ĉ-1 1 Td = ε∞ where Td = Ĉ2 ( εd ) εd , (

and according to the Green elasticity theory:

Td = 2µ 2 εd + K 2 - 2 3 µ 2 Tr ( εd ) Î + A ε2 d + B Tr ε2 Î + 2 ε Tr ( ε) + C Tr 2 ( ε) Î . ( 4.80)
Equation (4.80) applies to a single inclusion. Our study will concern a dispersion of inclusions, limited to the low-concentration case, where there is no interaction between the particles. So, we consider the average stress in the heterogeneous medium:

T = 1 V V T dv = 1 V Ĉ1 V m εdv + V e T dv = 1 V Ĉ1 V m εdv + V e T dv + 1 V Ĉ1 V m εdv - 1 V Ĉ1 V m εdv = Ĉ1 ε + c Td -Ĉ1 εd , (4.81)
where V m is the volume of the matrix, V e the volume of the inclusions and the total volume is V = V m + V e .

So:

T = 2µ 1 ε + K 1 - 2 3 µ 1 Tr ε Î +c (µ 2 -µ 1 ) εd + K 2 -K 1 - 2 3 (µ 2 -µ 1 ) Tr( εd ) Î + A ε2 d + B Tr ε2 d Î + 2 εd Tr( εd ) + C Tr 2 ( εd ) Î . ( 4.82)
Expression (4.82) is an exact one, for which we had to make no approximation. Now, in the case of a strongly dilute medium (c ≪ 1), we make the approximation that the elastic field inside the particles can be considered uniform, i.e. Td ≈ Td and εd ≈ εd , so:

T = c Td + Ĉ1 ( ε -c εd ) , (4.83)
and we make the approximation:

ε = c εd + (1 -c) ε∞ . (4.84)

Expressions for spherical inclusions

The Eshelby tensor for a spherical inclusion is [83]:

S ijkl = 1 15(1 -ν 1 )
δ ik δ jh + δ ih δ jk (4 -5ν 1 ) + δ kh δij (5ν 1 -1) , (

where ν 1 is the Poisson ratio of the matrix, so:

S ijkl ε d,kh = 2(4 -5ν 1 ) 15 (1 -ν 1 ) ε d,ij + 5(ν 1 -1) 15 (1 -ν 1 ) ε d,kk δ ij , (4.86)
where:

ν 1 = 3K 1 -2µ 1 2 (3K 1 + µ 1 ) , ( 4.87)
which gives us: Ŝ εd = 6 5

K 1 + 2µ 1 3K 1 + 4µ 1 εd + 1 5 3K 1 -4µ 1 3K 1 + 4µ 1 Tr ( εd ) Î . ( 4.88)
In this case, we can rewrite equation (4.79) as:

L εd + M Tr( εd ) Î + N ε2 d + O εd Tr( εd ) + P Tr ε2 d Î + Q Tr 2 ( εd ) Î = ε∞ , ( 4.89)
with:

L = 1 + 6 5 K 1 + 2µ 1 3K 1 + 4µ 1 µ 2 µ 1 -1 , (4.90)
M = 5K 2 -K 1 3 + 2 µ 2 µ 1 -4(µ 2 -µ 1 ) 5(3K 1 + 4µ 1 ) , ( 4.91
)

N = 3 5 A µ 1 K 1 + 2µ 1 3K 1 + 4µ 1 , ( 4.92)
O = 6 5 B µ 1 K 1 + 2µ 1 3K 1 + 4µ 1 , ( 4.93)
P = 1 15(3K 1 + 4µ 1 ) 15B -A 1 + 3 K 1 µ 1 , ( 4.94)
Q = 1 15(3K 1 + 4µ 1 ) 15C -2B 1 + 3 K 1 µ 1 . ( 4.95)
Using equation (4.84) and the expression of ε∞ given in equation (4.89), the expression of ε becomes:

ε = (c + (1 -c)L) εd + (1 -c) M Tr( εd ) Î + N ε2 d + O εd Tr( εd ) + P Tr ε2 d Î + Q Tr 2 ( εd ) Î , ( 4.96)
which we rewrite as:

ε = L ′ εd + M ′ Tr( εd ) Î + N ′ ε2 d + O ′ εd Tr( εd ) + P ′ Tr ε2 d Î + Q ′ Tr 2 ( εd ) Î , ( 4.97)
where:

L ′ = c + (1 -c)L M ′ = (1 -c)M N ′ = (1 -c)N O ′ = (1 -c)O P ′ = (1 -c)P Q ′ = (1 -c)Q. (4.98)
Now, we can calculate the five following expressions:

Tr ε = (L ′ + 3M ′ ) Tr( εd ) + (N ′ + 3P ′ ) Tr ε2 d + (O ′ + 3Q ′ ) Tr 2 ( εd ), ε 2 = L ′2 ε2 d + 2L ′ M ′ εd Tr( εd ) + M ′2 ε2 d Î , ε Tr ε = L ′ (L ′ + 3M ′ ) εd Tr( εd ) + M ′ (L ′ + 3M ′ ) Tr 2 ( εd ) Î , Tr ε 2 = L ′2 Tr ε2 d + M ′ (2L ′ + 3M ′ ) Tr 2 ( εd ), Tr 2 ε = (L ′ + 3M ′ ) 2 Tr 2 ( εd ). (4.99)
Inverting this system of equations, we obtain:

Tr 2 ( εd ) = Tr 2 ε (L ′ + 3M ′ ) 2 , Tr ε2 d = Tr ε 2 L ′2 - M ′ (2L ′ + 3M ′ ) L ′2 (L ′ + 3M ′ ) 2 Tr 2 ε , Tr( εd ) = Tr ε L ′ + 3M ′ - N ′ + 3P ′ L ′ + 3M ′         Tr ε 2 L ′2 - M ′ (2L ′ + 3M ′ ) L ′2 (L ′ + 3M ′ ) 2 Tr 2 ε         - O ′ + 3Q ′ (L ′ + 3M ′ ) 3 Tr 2 ε , εd Tr( εd ) = ε Tr ε L ′ (L ′ + 3M ′ ) - M ′ (L ′ + 3M ′ ) L ′ (L ′ + 3M ′ ) Tr 2 ε (L ′ + 3M ′ ) 2 Î . (4.100)
Now, we replace those expressions in that of T , i.e. equation (4.82). We obtain a constitutive equation for the homogenized material:

T = 2µ eff ε + K eff - 2 3 µ eff Tr ε Î + A eff ε 2 + B eff Tr ε 2 Î + 2 ε Tr εd + C eff Tr 2 ε Î , (4.101)
where:

µ eff = µ 1 + c µ 2 -µ 1 L ′ K eff = K 1 + c K 2 -K 1 L ′ + 3M ′ A eff = c A L ′2 -2c N ′ (µ 2 -µ 1 ) L ′3 B eff = 2c (N ′ M ′ -L ′ P ′ )(µ 2 -µ 1 ) L ′3 (L ′ + 3M ′ ) -c (N ′ + 3P ′ ) K 2 -K 1 -2 3 (µ 2 -µ 1 ) L ′2 (L ′ + 3M ′ ) + c B L ′2 C eff = c 9 9C + 9B + A (L ′ + 3M ′ ) 2 + c 9 A -3B L ′2 + c 9 (4N ′ + 6O ′ )(µ 2 -µ 1 ) L ′2 (L ′ + 3M ′ ) - 2c 9 3B -A L ′ (L ′ + 3M ′ ) + c 9 (3N ′ + 9P ′ )(K 2 -K 1 ) L ′2 (L ′ + 3M ′ ) - 4c 9 N ′ (µ 2 -µ 1 ) L ′3 - c 3 (9Q ′ + 3O ′ + 3P ′ + N ′ )(K 2 -K 1 ) (L ′ + 3M ′ ) 3 . (4.102)

Expressions for cylindrical inclusions

For parallel cylindrical inclusions, in a linear matrix, with the strain in a plane perpendicular to the cylinders, the expression of the Eshelby tensor gives:

Ŝ εd = 1 2 K 1 + 2µ 1 K 1 + µ 1 εd + 1 4 K 1 -2µ 1 K 1 + µ 1 Tr ( εd ) Î . ( 4.103)
Using the same procedure, for a 2D system, as the one used above for spherical inclusions, we obtain the following effective linear and quadratic Landau coefficients [START_REF] Giordano | Nonlinear elasticity of composite materials[END_REF]:

µ eff = µ 1 + c µ 2 -µ 1 L ′ K eff = K 1 + c K 2 -K 1 L ′ + 2M ′ A eff = c A L ′2 -2c N ′ (µ 2 -µ 1 ) L ′3 B eff = c [N ′ (µ 2 -µ 1 ) + BL ′ ] L ′3 - c(2P ′ + N ′ )(K 2 -K 1 ) L ′2 (L ′ + 2M ′ ) C eff = c 4C + 6B + A 4(L ′ + 2M ′ ) 2 + c A -2B 4L ′2 + c 2(N ′ + O ′ )(µ 2 -µ 1 ) + (2P ′ + N ′ )(K 2 -K 1 ) (2L ′2 (L ′ + 2M ′ ) -c 2P ′ + N ′ + 4Q ′ + 2O ′ )(K 2 -K 1 ) 2(L ′ + 2M ′ ) 2 -c N ′ (µ 2 -µ 1 ) L ′3 -c A + 2B 2L ′ (L ′ + 2M ′ ) , (4.104)
where:

L ′ = c + (1 -c) 1 + 1 2 K 1 + 2µ 1 K 1 + µ 1 µ 2 µ 1 - 1
M ′ = (1 -c) 2K 2 -K 1 1 + µ 2 µ 1 -2(µ 2 -µ 1 )
4(K 1 + µ 1 )

N ′ = (1 -c) A 4µ 1 K 1 + 2µ 1 K 1 + µ 1 O ′ = (1 -c) 2B 4µ 1 K 1 + 2µ 1 K 1 + µ 1 P ′ = (1 -c) 4B -A K 1 µ 1 8(K 1 + µ 1 ) Q ′ = (1 -c) 4C -2B K 1 µ 1 8(K 1 + µ 1 ) .
(4.105)

Derivation to the Third Order

Because of the complexity of the equations, the expressions for the quadratic and cubic nonlinearities will be given for all the terms for 2D systems, but only for the leading terms (corresponding to the longitudinal compression waves) for 3D systems.

We consider nonlinear inclusions of a nonlinear material noted 2 in a linear matrix noted 1. In 3D

If the material is isotropic, then the energy is an isotropic tensor and it can be written as [START_REF] Hamilton | Separation of compressibility and shear deformation in the elastic energy density (L)[END_REF]:

U = µ Tr( ε2 ) + λ 2 Tr 2 ( ε) + A 3 Tr( ε3 ) + B Tr( ε) Tr( ε2 ) + C 3 Tr 3 ( ε) + E Tr( ε) Tr( ε3 ) + F Tr 2 ( ε) Tr( ε2 ) + G Tr 2 ( ε2 ) + H Tr 4 ( ε). (4.106)
We rewrite this expression using the Einstein notation:

U = µε ij ε ji + λ 2 (ε ii ) 2 + A 3 ε ij ε jk ε ki + Bε ii ε kh ε hk + C 3 (ε ii ) 3 + Eε ii ε jk ε kq ε qj + F(ε ii ) 2 ε jk ε kj + G(ε ij ε ji ) 2 + H(ε ii ) 4 . ( 4.107)
We differentiate this expression to obtain T nm :

T nm = ∂U ∂ε nm = µδ in δ jm ε ji + µε ij δ jn δ im + 2λ 2 (ε jj )δ nm + A 3 δ in δ jm ε jk ε ki + A 3 ε ij δ jn δ km ε ki + A 3 ε ij ε jk δ kn δ im + Bδ nm ε kh ε hk + Bε ii (δ kn δ hm ε hk + ε kh δ hn δ km ) + 3C 3 (ε ii ) 2 δ nm
+ Eδ nm ε jk ε kq ε qj + Eε ii δ jn δ km ε kq ε qj + ε jk δ kn δ qm ε qj + ε jk ε kq δ qn δ jm

+ 2F(ε ii )δ nm ε jk ε kj + F(ε ii ) 2 δ jn δ km ε kj + ε jk δ kn δ jm + 2Gε pk ε kp δ in δ jm ε ji + ε ij δjnδ im + 4H(ε ii ) 3 δ in δ jm , (4.108)
which we write more simply: We regroup the terms of the previous expression and we define:

T nm = 2µε nm + λε ii δ nm + Aε mk ε kn + Bδ nm ε kh ε hk + 2Bε ii ε nm + C(ε ii ) 2 δ nm + Eδ nm ε jk ε kq ε qj + 3Eε ii ε nq ε qm + 2Fε ii δ nm ε jk ε kj + 2F(ε ii ) 2 ε nm + 4Gε pk ε kp ε nm + 4H(ε ii ) 3 δ nm .
e = B + A 2 3f = C - A 2 g = G h = 4H -2E l = 2F + 3E. (4.112)
The final expression for T in 2D is:

T nm = 2µε nm + λε ii δ nm + eε hk ε kh + 2eε ii ε nm + 3f (ε ii ) 2 δ nm + 4gε pk ε kp ε nm + h(ε ii ) 3 δ nm + l(ε ii ) 2 ε nm + lε ii δ nm ε jk ε kj , ( 4

Expressions for the Quadratic and Cubic Nonlinearities

To take into account both quadratic and cubic nonlinearities, we perform the same derivation as in section 4.5.1, but this time, we write expressions up to order three.

The full expressions with all the terms up to order three are extremely long to write in 3D. So, we perform the derivation up to the third order in 2D for all the terms and in 3D for the leading terms (corresponding to the longitudinal compression waves) only. Those derivations are provided in the appendix.

As in section 4.5.1, the homogenization procedure consists in writing the strain in the composite system ε∞ as a function of the strain inside the nonlinear inclusions εd and the Eshelby's tensor Ŝ:

ε∞ = εd + Ŝ Ĉ-1 1 Td -εd (4.115)
Then, ε∞ is expressed as a linear combination of the strain inside the nonlinear inclusions εd and its trace: For a low concentration c ≪ 1 of nonlinear inclusions, we make the following approximations: Td ≈ Td , εd ≈ εd , which mean that the elastic fields inside the particles is considered uniform. We also consider that Finally, we know all the necessary expressions to write:

Û                                            
Tr 3 ( εd ) Î                                             =                                             ε Tr ε Î ε Tr ε Tr ε 2 Î Tr 2 ε Î ε Tr ε 2 Tr ε Tr ε 2 Î ε Tr 2 ε Tr 3 ε Î                                             , ( 4.118)
with an explicit expression of the matrix Û . So, the effective parameters are given by the following expression:

                                            2µ eff K eff -µ eff 2e eff e eff 3f eff 4g eff l eff l eff h eff                                             =                                             2µ 1 K 1 -µ 1 0 0 0 0 0 0 0                                             T + c                                             2µ 2 -2µ 1 K 2 -µ 2 -K 1 + µ 1 2e e 3f 4g l l h                                             T Û-1 , (4.119)
where Û is known. The matrix Û can be inverted analytically in Maple, but the expressions of the effective parameters are too big to be written here. For simulation purposes, Û is inverted numerically in Octave or numpy. For the following study of the amplification of the nonlinear effective parameters, we use simplified expressions of the effective parameters: for each effective parameter that we study, we consider that the corresponding nonlinear parameter is non-null but that all other nonlinear parameters are null.

Amplification of the Nonlinear Effective Parameters

We consider a dispersion of nonlinear particles embedded in a linear matrix. We work in plain strain conditions, with a 2D system. Both particles and matrix are isotropic. Therefore, the matrix is described by: T1 = 2µ In the case of small concentration c ≪ 1, if we only consider longitudinal compression wave, then the stress-strain relation for the effective material obtained by an homogenization procedure is:

T = 2µ eff ε + (K eff -µ eff ) Tr( ε) Î + 3f eff Tr 2 ( ε) Î + h eff Tr 3 ( ε) Î . (4.122)

First observations

We are now trying to choose appropriate materials for the matrix and the inclusions as well as a specific concentration in order to observe an amplification of the nonlinear parameters, namely obtain an effective material with higher nonlinearities than that of the nonlinear inclusions. As we are more interested in longitudinal compression waves, we will concentrate on parameters f and h in 2D, or F and H in 3D. We want to obtain f eff /K eff > f /K 2D , h eff /K eff > h/K 2D and so on...

To obtain such amplifications, we will see that a high compressibility contrast is needed between the matrix and the inclusions

K 1 ≫ K 2 and µ 1 ≈ µ 2 ≪ K 2 . If, for instance: K 1 = 100K 2 and µ 1 ≈ µ 2 ≈ K 2 /10,
then, the inclusions are 100 times more compressible than the matrix, which can be qualified as hard and incompressible and the speeds of sound in the matrix verify c L ≫ c t .

Study of the quadratic amplifications

We suppose that f 0 inside the particles and that all other nonlinearities are null. We obtain the following linear and nonlinear effective parameters:

K eff = K 1 + c K 2 -K 1 c + (1 -c) µ 1 +K 2 µ 1 +K 1 µ eff = µ 1 + c µ 2 -µ 1 c + (1 -c) 1 + 1 2 µ 2 µ 1 -1 K 1 +2µ 1 K 1 +µ 1 f eff = cf (K 1 + µ 1 ) 3 [c(K 1 -K 2 ) + K 2 + µ 1 ] 3 . (4.123)
We are interested in determining the maximum value of the ratio f eff /K eff compared to f /K 2 . Those values represent the nonlinear parameter of the effective material and that of the inclusions, respectively. If we solve the equation d dc f eff K eff = 0, we can find the optimal volume fraction:

c f = K 1 4µ 1 µ 1 + K 2 K 1 -K 2 1 -1 -8 µ 1 K 1 . (4.124)
We note that c f exists only if µ 1 < K 1 /8. In general, we know that 0 < µ 1 < 3K 1 and therefore we are working on a very small region of the admitted zone (µ 1 varies from 0 to 3K 1 and the Poisson's ratio ν 1 is negative for µ 1 between K 1 and 3K 1 ).

We now determine the amplification of the quadratic nonlinearity, defined by: (f eff /K eff )/(f /K 2 ) and we prove that it is a function of x = µ 1 /K 1 and y = K 1 /K 2 . We finally obtain: where the leading term is y/4 -xy 2 /2. We have therefore an amplification proportional to y = K 1 /K 2 . Under the hypothesis x → 0 and y → ∞, we can simplify expression (4.124):

A f = f eff K eff f K 2 = 16(1 + x) 3 (1 - √ 1 -8x)xy 2 (1 + xy) 2 (y -1)(3 + √ 1 -8x)( √ 1 -8x -1 -4x) 2 . ( 4
c f ≈ x + 1 y . (4.127) Finally, if µ 1 ≪ K 1 , K 2 ≪ K 1 and µ 1 ≪ K 2 ,
the optimal concentration is:

c f = x + 1 y = µ 1 + K 2 K 1 ≈ K 2 K 1 , ( 4.128)
which corresponds to an amplification of:

A f = y 2 1 2 -xy = K 1 2K 2 1 2 - µ 1 K 2 ≈ K 1 4K 2 . (4.129)
We notice that the expressions of A f and c f are similar to that obtained for the quadratic terms of the mass-spring system and the fluid system.

Study of the cubic amplifications

We now consider h 0 within the particles and all other nonlinear constants are supposed null. We have:

h eff = ch(K 1 + µ 1 ) 4 [c (K 1 -K 2 ) + K 2 + µ 1 ] 4 . (4.130)
The optimal volume ratio, obtained by differentiating the previous expression according to c is:

c h = K 1 3µ 1 µ 1 + K 2 K 1 -K 2 1 -1 -3 µ 1 K 1 . (4.131) c h exists only if µ 1 < K 1 /3.
The ratio (h eff /K eff )/(h/K 2 ) can be written in terms of x and y as:

A h = h eff K eff f K 2 = 27(1 + x) 4 √ 1 -3x -1 x 2 y 3 (xy + 1) 3 (y -1) 2 + √ 1 -3x √ 1 -3x -1 -3x 3 . (4.132)
We develop its Taylor series expansion for x → 0 and y → ∞: where the leading term is 4y 2 /27 -4xy 3 /9. Similarly, we can simplify the expression for c h :

A h =
c h ≈ 1 2 x + 1 y = c f 2 . (4.134) Finally, if µ 1 ≪ K 1 , K 2 ≪ K 1 and µ 1 ≪ K 2 ,
the optimal volume fraction is:

c h = 1 2 x + 1 y = µ 1 + K 2 2K 1 ≈ K 2 2K 1 , ( 4.135)
which corresponds to an amplification of:

A h = 4y 2 9 1 3 -xy = 4K 2 1 9K 2 2 1 3 - µ 1 K 2 ≈ 4K 2 1 27K 2 2 . ( 4.136)
Once again, the expressions of A h and c h are similar to that obtained for the cubic terms of the mass-spring system and the fluid system.

Amplifications for 3D systems

We perform the same studies on the nonlinear effective constants for a 3D system in order to determine the amplification factors (C eff /K eff )/(C/K 2 ) and (H eff /K eff )/(H/K 2 ) as a function of x = µ 1 /K 1 and y = K 1 /K 2 .

For the second order, the optimal concentration that we obtain is:

c C = 1 y 1 + 4 3 xy ≈ 1 y , ( 4.137)
and the amplification factor is:

A C = C eff K eff C K 2 = y 4 1 - 8 3 xy ≈ y 4 . (4.138)
For the third order, the optimal concentration that we obtain is:

c H = 1 2y 1 + 4 3 xy ≈ 1 2y , (4.139) d = 2 d = 3 g = 2 c f ≈ x + 1 y ≈ 1 y c C ≈ 1 y 1 + 4 3 xy ≈ 1 y A f ≈ y 2 1 2 -xy ≈ y 4 A C ≈ y 4 1 -8 3 xy ≈ y 4 g = 3 c h ≈ 1 2 x + 1 y = c f 2 ≈ 1 2y c H = 1 2y 1 + 4 3 xy ≈ 1 2y A h ≈ 4y 2 9 1
3 -xy ≈ 4y 2 27

A H = 4y 2
27 (1 -4xy) ≈ 4y 2 27 Table 4.4 -Optimal concentrations and amplification factors for the nonlinear parameters.

and the amplification factor is:

A H = H eff K eff H K 2 = 4y 2 27 (1 -4xy) ≈ 4y 2 27 . (4.140)
Those results are summed up in table 4.4. As for 2D systems, we obtain the same optimal concentrations and amplification coefficients as for the mass-spring and fluid systems.

Summary

To summarize the expressions of the nonlinear amplifications and optimal concentrations given above, for 2D and 3D systems and for second-and third-order terms, we have defined generic expressions that take the dimensionality d (d = 2 for 2D systems, d = 3 for 3D systems) and the degree g (g = 2 for second-order leading term, g = 3 for third-order terms) as parameters:

c opt = 1 + bxy (6 -g)bx(y -1) 1 - √ 1 -ax , ( 4.141)
for the optimal concentration and: Amp = 4g 9-2g (1 + bx) g+1 (1 + bxy) g × So, those expressions will give the optimal concentrations and corresponding nonlinear amplifications:

y g y -1 b 2 × x 1 - √ 1 -ax g-1 × 1 5 -g + √ 1 -ax g+1 , ( 4
• (f eff /K eff )/(f /K 2 ) for d = 2 and g = 2,

• (h eff /K eff )/(h/K 2 ) for d = 2 and g = 3,

• (C eff /K eff )/(C/K 2 ) for d = 3 and g = 2,

• (H eff /K eff )/(H/K 2 ) for d = 3 and g = 3.

If, moreover, we make the hypothesis x → 0, x > 0 and y → ∞, the optimal ratio x * is given by: dAmp(x) dy = 0. (4.144)

We obtain:

x * → g -1 y × 2d -3 3d -5 , ( 4.145)
and the corresponding amplification factor is:

Amp(x * ) ≈ 4g 9-3g y g-1 (18 -5g) g-1 (6 -g) g+1 , (4.146) and so the optimal concentration becomes: We can interpret the amplification of the nonlinear parameters as follows: because the matrix is composed of an incompressible medium, and the inclusions of a less-incompressible medium, the stress applied to the whole material will concentrate around the inclusions.

c opt (x * ) ≈ g 2y × 18 -5g 6 -g . ( 4
So, for a level of stress τ 0 applied to the whole material, the inclusions will experiment a higher level of stress τ 1 > τ 0 , because this stress will be concentrated on smaller volumes. Consequently, the inclusions will show a nonlinear behavior corresponding to this level of stress τ 1 . The effective material appears highly nonlinear because under a stress τ 0 , it shows a behavior corresponding to that of the inclusions at stress level τ 1 > τ 0 . Finally, the behavior of the effective material appears more nonlinear than that of the nonlinear inclusions, for the same stress level τ 0 .

Values of the Constants

Values of constants of soft and hard materials are given in tables 4.5 and 4.6 respectively.

Moreover, table 4.7 provides a list of publications giving the linear and nonlinear elastic constants of some materials.

Exploitation

We consider a 2D system composed of inclusions of a (virtual) porous polymer in PDMS. The density, longitudinal and transverse speeds and compressibility constants for those materials are indicated in PDMS was chosen because of its high contrast between longitudinal and transverse speeds:

K 1 ≫ µ 1 .
The porous polymer has been chosen because of its low longitudinal speed compared to that of PDMS:

K 1 ≫ K 2 .
Such a porous polymer can be produced with the High Internal Phase Emulsion (HIPE) technique. Therefore, its longitudinal and transverse speeds are realistic and are deduced from [START_REF] Gokmen | Fabrication of porous "clickable" polymer beads and rods through generation of High Internal Phase Emulsion (HIPE) droplets in a simple microfluidic device[END_REF][START_REF] Brunet | Soft Acoustic Metamaterials[END_REF]. Its density is that of typical polymers (see table 4.5) but does not take into account the porosity, which value is unknown.

A study of nonlinear amplifications with a more realistic porous polymer has been performed and published after this thesis was written [START_REF] Guerder | Tuning the elastic nonlinearities in composite nanomaterials[END_REF].

The nonlinear parameters of the porous polymer are defined as follows: f = 10 × K 2D = 1.28 GPa and h = 100 × K 2D = 12.8 GPa.

The effective linear (µ eff and K eff ) and nonlinear (f eff and h eff ) parameters of the composite system are calculated. We observe that the optimal concentration for cubic amplification is 0.0067. For this concentration, the cubic amplification reaches 1042.6 and the quadratic amplification at this concentration is 4.69, as shown on

Conclusion

In this chapter, we have provided equations and interpretation for the amplification of quadratic and cubic nonlinearities in systems composed of nonlinear inclusions in a linear matrix.

We have shown that those parameters are amplified relatively to the compressibility ratio of the constituent materials and the characteristics of the matrix, and we have provided equations to predict those amplifications as well as the optimal inclusion concentration to reach this amplification.

Since it is possible to obtain a greater amplification for the cubic elastic constant compared to that of the quadratic elastic constant, our goal is to use a composite material as the defect in the nonlinear elastic resonator described in chapter 3.

In the present chapter, we have derived static equations, independent of any time-dependent phenomena such as wave propagation. In our dynamic systems like the bone structure and nonlinear resonator, the stress is applied to the constituent by the means of longitudinal elastic wave propagation. This is why the next chapter will deal with the dynamic behavior of a composite nonlinear material designed using the equations of this chapter. We will try to verify that the amplification of the nonlinear constants can also be obtained in the context of wave propagation through the composite material.

Chapter 5

Numerical Study of 2D Nonlinear Phononic Crystals

Introduction

In chapter 4, we have performed a derivation of the mixing laws of the quadratic and cubic nonlinear parameters in 2D and 3D and proved that for some concentrations of nonlinear inclusions in a linear matrix, those nonlinear parameters could be considerably amplified. We have provided an example of a system with cubic amplification, consisting in inclusions of a nonlinear porous polymer in a linear PDMS matrix.

The homogenization procedure used in the previous chapter and its results for the porous polymer inclusions in PDMS are valid for a static system and do not take into account any dynamic phenomenon such as wave propagation or the application of a periodic stress function to the system.

We will now verify, through numerical simulation of elastic wave propagation, that the dynamic behavior of a composite nonlinear system still follows those rules, i.e. that elastic wave propagation through a finite section of the heterogeneous system described in section 4.6.3 is close enough to that in the effective medium obtained with the homogenization procedure.

We will first describe Hedge, the framework which we used to perform the numerical simulations, and the Elastodynamics Operator that we implemented in this framework for the purpose of simulating elastic wave propagation in heterogeneous nonlinear media. Then, we will present the models and the results of the simulations and conclude.

The Discontinuous Galerkin Finite Elements Method

In order to perform a numerical study of 2D nonlinear PnCs, one needs to solve the elastic wave equation numerically, including for nonlinear problems and for materials with a complex shape. The number of available methods to solve this problem is relatively small, which makes the Discontinuous Galerkin (DG) method particularly interesting. First the principle and improvements of the DG method will be described, then its alternatives and finally its implementation on graphic processors.

Overview of the Discontinuous Galerkin method

In 2D, the space is divided into a mesh of contiguous triangular or quadrilateral elements (in 3D, those are tetrahedral volumes). The Finite Elements Method (FEM) approximates the unknown solution of the differential equation with polynomials on each element of the grid but does not allow a local statement on each element [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods[END_REF].

In the Finite Volume Method (FVM), the equation is solved by calculating the local average on each of those volumes. It is an iterative process and the value at t n only depends on the value at t n-1 and on the fluxes that enter and leave the volume. This method is rather simple to implement but does not allow high-order accuracy [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods[END_REF].

The DG method combines the advantages of those two methods: using polynomials of arbitrary order of the FEM and using numerical fluxes at the interfaces of the FVM. The main difference between DG and the FEM is that polynomials must be continuous at the interface of the elements in the FEM whereas the solution can be discontinuous across the element interfaces in DG, which allows the incorporation the well-established ideas of numerical flux functions from the FVM. So, by choosing carefully the values of those fluxes, it is possible to minimize the error introduced by the approximation [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods[END_REF].

The DG method was first introduced by Reed and Hill in 1973, to solve hyperbolic partial differential equations in the context of the neutron transport equation. It is also applied to elliptic problems, but those don't appear in elastodynamics (only in elastostatics where there is no time dependency) [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF].

Solving the equation with DG can be done in two ways: modal (decomposition of the local solution in a linear combination of several resonating modes) or nodal (resolution at chosen points before interpolating by decomposing the local solution in a linear combination of orthogonal polynomials). Moreover, the method can be p-adaptive (the order of the polynomials can vary from one element to another), h-adaptive (the element size can vary locally) or hp-adaptive (both can vary) [START_REF] Etienne | An hp-adaptive discontinuous Galerkin finiteelement method for 3-D elastic wave modelling[END_REF].

Application of the DG-FEM to elastodynamics

The 2D elastic wave equations can be written as a conservation equation which links the time derivative of a state vector Q to the space derivatives of the flux vectors F and G and to the source vector S:

∂Q(t, x) ∂t = ∂F(t, x) ∂x + ∂G(t, x) ∂y + S(t, x), (5.1)
where F ij are the components of the deformation gradient, P ij are the components of the Piola-Kirchoff stress tensor, v i the components of the particle velocity vector and x = [x, y] is the direction vector. This form makes the DG resolution easier, especially for the nodal method described by Hesthaven and Warburton and also facilitates the inclusion of nonlinear effects [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods[END_REF].

In computational seismology, Käser and Dumbser proposed to use the DG-FEM combined with the ADER (Arbitrary high-order DERivatives) method, in order to solve the equation with high orders of accuracy (practically up to order 10 for polynomials), both in time and in space. They explain that a recurrent issue is about simulating the absorbing (where the wave is not reflected) and the free surface boundaries and show that it can be done by imposing the numerical fluxes [START_REF] Käser | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -I. The two-dimensional isotropic case with external source terms[END_REF]. Then, de la Puente et al. extended this method in 2007 to 3D anisotropic materials and coupled it with viscoelastic effects [START_REF] De La Puente | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -IV. Anisotropy[END_REF].

In 2010, the ADER-DG method was improved by Castro et al. who added space-variable coefficients that allow to deal with heterogeneous materials, including inside the same cell of the mesh. However, this improvement only concerns 2D for now and does not deal well with discontinuities inside a cell [START_REF] Castro | Seismic waves in heterogeneous material: subcell resolution of the discontinuous Galerkin method[END_REF].

Bou Matar et al. included in 2012 the quadratic and cubic nonlinearities and the quadratic hysteretic nonlinearities by writing the equations of nonlinear elastodynamics in a conservative form. They also introduced equations for attenuation and absorbing boundaries, called Nearly Perfectly Matched Layers (NPML), and validated the method by comparing it to analytical results [START_REF] Bou Matar | A nodal discontinuous Galerkin finite element method for nonlinear elastic wave propagation[END_REF].

The application of DG-FEM to the 3D elastodynamics equation was analyzed by Delcourte and Glinsky-Olivier in 2013. They demonstrated its numerical stability and studied its convergence [START_REF] Delcourte | Analysis of a discontinuous Galerkin method for elastodynamic equations. application to 3D wave propagation[END_REF].

DG-FEM GPU implementation

Graphical Processor Units (GPUs) were initially dedicated to image processing for the graphical output of computers. They are now more and more used to general purpose computation, and some constructors sell computation-dedicated graphic cards with no video output, such as nVidia Tesla, which can perform 10 12 floating-point operations per second (1 TFLOPS).

Implementation of an Elastodynamics Operator

In order to solve for the propagation of elastic waves in a nonlinear heterogeneous medium, we implemented an Elastodynamics Operator for the framework Hedge.

First, we implemented the resolution of the equations of acoustics in a homogeneous linear material.

In the linear Elastodynamics Operator, the Piola-Kirchoff stress tensor P is calculated with:

P ij = C ijkl ∂u k ∂x l (5.3)
and the deformation gradient is:

F kl = ∂u k ∂x l + δ kl .
(5.4)

Then, we added a possibility to create heterogeneities in the medium and we added the quadratic nonlinear behavior of the materials composing the matrix and the inclusions. This way, the quadratic Elastodynamics Operator can solve wave propagation in heterogeneous anisotropic media in 1D, 2D and 3D.

In the nonlinear quadratic Elastodynamics Operator, the Piola-Kirchoff stress tensor P is calculated with:

P ij = C ′ ijkl ∂u k ∂x l , ( 5.5)
where:

C ′ ijkl = C ijkl + 1 2 M ijklmn ∂u m ∂x n and M ijklmn = C ijklmn + C ijlm δ kn + C ilmn δ jk + C iklm δ jn . ( 5.6)
Finally, we implemented a cubic nonlinear Elastodynamics Operator, which deals with both quadratic and cubic elastic nonlinearities. For the sake of simplicity, we only implemented this cubic Elastodynamic Operator for heterogeneous isotropic materials in 2D only. For this operator, following equation (4.113), the Piola-Firchoff stress tensor P is defined as (in Einstein notation):

P nm = 2µε nm + λε ii δ nm + eε hk ε kh + 2eε ii ε nm + 3f (ε ii ) 2 δ nm + 4gε pk ε kp ε nm + h(ε ii ) 3 δ nm + l(ε ii ) 2 ε nm + lε ii δ nm ε jk ε kj , (5.7)
where λ and µ are the linear elastic constants, e and f the quadratic elastic constants and g, h and l the cubic elastic constants. The components of the strain tensor are given by:

ε ij = 1 2 ∂u i ∂x j + ∂u j ∂x i . ( 5.8)
For our simulations, we only take into account the constants responsible for longitudinal compression waves, namely f and h and consequently, the expression is simplified (in Einstein notation):

P nm = µ ∂u n ∂x m + ∂u m ∂x n + λ ∂u i ∂x i δ nm + 3f ∂u i ∂x i 2 δ nm + h ∂u i ∂x i 3 δ nm .
(5.9)

The Nearly Perfectly Matched Layers

The PML (Perfectly Matched Layers) are perfectly matched absorbing layers. In theory, no wave can be reflected on it, whatever its incidence, its polarization or its frequency. In practice, some reflection can appear in a discretized structure. Such absorbing conditions are particularly useful when simulating infinite media (the wave seems to cross the border and does not come back).

In the context of the DG-FEM, a special type of PML, namely NPML (Nearly Perfectly Matched Layers), was proposed by Hu and Cummer [START_REF] Hu | The nearly perfectly matched layer is a perfectly matched layer[END_REF][START_REF] Cummer | A simple, nearly perfectly matched layer for general electromagnetic media[END_REF] for electromagnetism and acoustics [START_REF] Hu | Application of the nearly perfectly matched layer in acoustic wave modeling[END_REF]. A publication by Bou Matar et al. [START_REF] Bou Matar | A nodal discontinuous Galerkin finite element method for nonlinear elastic wave propagation[END_REF] proposes an improvement of NPML. It is that improvement that I used to implement the NPML in the Elastodynamics Operator.

F x + F ′′ x κ x - ∂ ∂y F y + F ′′ y κ y - ∂ ∂z F z + F ′′ z κ z = S(t), ( 5.11)
where the auxiliary variables F ′′ i are calculated with the following ordinary differential equations:

∂F ′′ i ∂t = -α i F ′′ i - σ i κ i F i + F ′′ i for all i ∈ {x, y, z} , (5.12)
and:

α i = 1.4π d i -i + i 0 d i and κ i = 1 and σ i = 6 C 11 ρ d i i -i 0 d i 2 , ( 5.13)
where d i is the width of the NPML and i 0 ∈ {x 0 , y 0 , z 0 } is its position.

Validation of the Elastodynamics Operator

We have simulated the propagation of an elastic wave in a 2D structure with a rectangular shape, represented on Figure 5.1 using the Elastodynamics Operator and compared it with the analytical result computed by the Fortran program EX2DDIR. This program by Berg and Flemming implements an algorithm given by Berg et al. [START_REF] Berg | Modeling the earth for oil exploration[END_REF] which uses the Cagniard-de Hoop method [START_REF] Hoop | A modification of Cagniard's method for solving seismic pulse problems[END_REF] to calculate the exact (analytical) response from a vertical directional point source in an elastic half-space with a free surface.

We use two receivers that both record the horizontal and vertical displacements. On

s(t) = (0.5 + a 1 (t -t D )
2 )e a 1 (t-t D ) 2 .

(5.14) The figures show a excellent agreement between the analytical solution and the numerical results obtained with the Elastodynamics Operator at order five.

We have evaluated the agreement between the analytical and numerical solutions for all orders between 1 and 9, which gives the convergence of the numerical method. This study, represented on Figure 5.4 shows that the difference between the numerical and analytical methods decrease from order 1 to order 4 and then 5.1 -Constants for the effective medium. Properties of the effective medium described in section 4.6.3 used to validate the measurement of the quadratic and cubic nonlinearities in Hedge. This structure is simulated using the CUDA version of the operator. The excitation is produced by a modulated sine function which expression is:

Numerical Studies on 2D Structures

s(t) = A × sin (2πf c t) × exp -((t-24/f c )f c /8) 2 , ( 5.15)
where A is adapted to produce a particle velocity of ~2 mm s -1 and f c = 40 Hz. The particle velocity is recorded at twenty positions. A Fourier transform of these time domain signals is made to extract the amplitudes of the fundamental (at 40 Hz), second and third hamonics (at 80 Hz and 120 Hz, respectively)

of the emitted wave. The obtained evolution of the second harmonic component as a function of distance is plotted on Figure 5.6 (stars). This evolution can be compared to the expected one. Indeed, in the case of a longitudinal plane wave in an isotropic solid, the generation of the second harmonic component P 2 is

Here, we would like to determine the nonlinear parameters f and h, defined in equation (4.112), as a function of β and η in the case where A l ≪ B l ≪ C l . From equation (4.112) we have:

3f = C - A 2 ≈ C ≈ - B l 2 , (5.18) h = 4H -2E = 4 24 (A l + C l ) - 8 3 A l ≈ C l 6 , ( 5.19)
leading to:

β = 1 + B l 2A l ≈ - C K ≈ - 3f ρ 0 c 2 l , ( 5.20
)

η = 1 + C l 2A l ≈ - C l 2K ≈ 3h ρ 0 c 2 l .
(5.21)

Plotting the evolution of the second harmonic component as a function of distance x calculated with equation (5.16) (blue line), with the value of β given by equation (5.20), on Figure 5.6, we find that the numerical and theoretical results are in excellent agreement. Now, measuring the amplitude of the fundamental component on the source and the amplitude of the second harmonic component, it is possible to determine the quadratic nonlinear parameter:

β = ρ 0 c 3 0 πf c × P 2 (x) xP 2 1 (0) . ( 5.22)
The calculated value for the nonlinear quadratic parameter β from the numerical results obtained with Hedge is in perfect agreement with the theoretical value, i.e. β = 139.

The same procedure can be used to measure the third order nonlinear parameter (h or η). The generation of the third harmonic component P 3 is governed by a more complex expression, but has been obtained by Liu et al. as [START_REF] Liu | Theoretical and experimental study of the third-order nonlinearity parameter C/A for biological media[END_REF]:

P 3 (x) = πf c 2ρ 2 0 c 5 0 P 3 1 (0) 3πf c c 0 β 2 x 2 + 3(η -1)x , ( 5.23)
where η = 1 + C l / (2A l ) is the cubic nonlinear coefficient of fluid. Now, when the quadratic nonlinear parameter is null (β = 0), then the amplitude of the third harmonic component increases linearly with the distance x:

P 3 (x) = πf c 2ρ 2 0 c 5 0 P 3 1 (0)3(η -1)x. ( 5.24)
Thus, as we did for the quadratic nonlinear parameter, we determine the cubic nonlinear parameter from the measurement of the evolution of the third harmonic component as a function of the distance x:

η = 4ρ 2 0 c 5 0 3ω × P 3 (x) xP 3 1 (0) - ωβ 2 x 2c 0 , ( 5.25)
where ω = 2πf c is the pulsation. Nevertheless, it is important to notice that the contribution of the cubic nonlinear parameter to the third harmonic component is in general quite small. In numerical experiment it is not a big issue, as we can suppress the quadratic nonlinear parameter contribution, i.e. setting it equal to zero, but it is a strong limitation for an experimental measurement of this parameter. Thus it is clear that the quadratic nonlinearity plays an important role in the generation of the third harmonic. This is consistent with the simulations performed in chapter 4.

First, we compare the evolution of the third harmonic component, when the quadratic nonlinear parameter is null (β = 0), as a function of distance x calculated with equation (5.24) (blue line) and with Hedge (stars), as shown on Figure 5.7. In this case, with β = 0, the third harmonic component increases linearly with the distance as it was the case for the second harmonic component. Once again, the numerical and theoretical results are in excellent agreement, and the calculated third order nonlinear parameter η = -60024 is the equal to the theoretical value.

Finally, plotting the evolution of the third harmonic component as a function of distance x calculated with equation (5.23) (blue line) and with Hedge (stars), considering both second and third order nonlinearities, we obtain excellent agreement between theoretical and numerical results, leading to the same value of the third order nonlinear parameter. We need to work with low frequencies, in order to have λ ≫ a ≫ r. This way, the inclusions will be "invisible" to wave propagation. To select this low frequency, we calculate with a Plane Wave Expansion method [START_REF] Deymier | Acoustic Metamaterials and Phononic Crystals[END_REF] the band structure of the considered PnC made of 0.04618 m radius inclusions of porous polymer in a PDMS matrix, and positioned on a square lattice with lattice parameter 1 m. The results displayed on Figure 5.10 show that we need to select a frequency below 50 Hz if we want to stay in the first Brillouin zone for the shear elastic waves propagating in the structure. Thus, the following simulations were performed with a center frequency of 20 Hz. This simulation is performed using the CUDA version of the operator. The particle velocity is recorded at eighteen positions spaced by 5 m. An example of the particle velocity obtained 90 m away from the source is represented on Figure 5.11. The excitation is produced by a modulated sine function which expression is given by equation (5.15), where A is adapted to produce a particle velocity of ~5 mm s -1 and f c = 20 Hz. The particle velocities calculated in the heterogeneous and in the effective homogeneous media overlap. The propagation in the considered heterogeneous medium seems to be similar to the one in a homogeneous medium, at the considered frequency, i.e. 20 Hz. However, in the inset of Figure 5.11, small time delay and maximal particle velocity can be seen. We will now evaluate the values of f eff and h eff and compare them with what was expected. To estimate these values of the quadratic and cubic nonlinear parameters, we use as in the previous section the value P 2 (x) and P 3 (x) of the second and third harmonics, that we measure on the transmission spectrum from Figure 5.12, where x is the length of the transmission. In order to verify the amplification of the nonlinear quadratic and cubic parameters in the composite system modeled in our simulations, we need to estimate the values of those nonlinear parameters, calculate the amplification ratio and compare it to the expected amplification calculated in chapter 4, section 4.6.

Conclusion

In this chapter we have presented a numerical scheme based on the Discontinuous Galerkin Method that is able to simulate nonlinear elastic wave propagation up to third order in heterogeneous solids of complex shape. We implemented this numerical scheme in the framework Hedge.

First, we validated this implementation on benchmarks for a linear case by comparing the results of the simulations with an analytical result. Then, we validated it for nonlinear cases [START_REF] Bou Matar | A nodal discontinuous Galerkin finite element method for nonlinear elastic wave propagation[END_REF]. For this validation, we compared the evolution of the second and third harmonic components recorded in Hedge with the prediction of their evolution according to expressions from the literature.

Those observations prove that the mixing law that we developed in chapter 4 is still valid in the scheme of a propagative system, even though we developed those equations using a static external stress.

A first consequence is that we proved the possibility for one to design composite materials with enhanced quadratic and cubic nonlinearities, of with a magnified ratio between cubic and quadratic nonlinearities.

Such a material could be used in the resonating cavity of the elastic resonator we designed in chapter 3 to improve its figure of merit, i.e. the relation between the width of the hysteresis and the minimal power needed to obtain that hysteresis.

The developed numerical scheme will also enable us to extend the study of nonlinear PnCs to higher dimensions, i.e. in two or three-dimension systems, due to the ability of this implementation to benefit from GPU and cluster-computing techniques.

General Conclusion

Conclusion

In this thesis, we have investigated the field of nonlinear phononic crystals and observed new phenomenon such as dispersion of elastic waves through four-phonon interaction as well as bistable switching of an elastic wave by a nonlinear resonator. Those phenomenon are due to the nonlinearities of the materials that constitute the crystals. Some properties of the second-and third-order nonlinear constants of materials were studied, such as how their values evolve when a nonlinear material is included in a linear matrix.

First, the PnCs were introduced and the state of the art was presented. We noted that very little work was done concerning the study of nonlinear nongranular PnCs. Notably, the influence of the nonlinear elastic parameters on the behavior of PnCs had never been studied.

Then, chapter 1 introduced preliminary calculations for the study of superlattices, i.e. 1D PnCs. Analytical results were demonstrated for the calculation of the dispersion curves of linear layered structures. This chapter also introduced numerical methods which are frequently used for simulating elastic wave propagation and calculating the energy density in 1D and 2D structures.

The results from chapter 1 were used in chapter 2 in which we studied elastic wave dispersion in a bone structure. We demonstrated that the stress-strain relation for a hydroxy apatite-collagen superlattice is highly nonlinear due to the ability of collagen to absorb and desorb water. We showed that this nonlinear behavior was responsible for energy dissipation through four-phonon interaction in correctly hydrated bones.

In chapter 3, we designed a nonlinear bistable switch based on a nonlinear elastic resonator. We used the third-order (cubic) elastic constants of a steel inclusion to shift the resonance frequency of a resonating cavity.

The low or high transmission through the structure depends on the particle velocity in the incident wave.

As the quadratic nonlinearities break this process, we studied, in chapter 4, the ways to lower the secondorder and increase the third-order elastic constants in a composite material. Using mixing laws obtained by a homogenization procedure, we calculated the values of those parameters for the material equivalent to a system of nonlinear inclusions in a linear matrix. We showed that for specific concentrations, that we calculated, the main cubic nonlinear elastic constant can be increased by a factor ~1000 while the quadratic constant is only increased by a factor ~5.

Finally, in chapter 5, we exploited the results of chapter 4 in 2D numerical simulations of wave propagation through a composite material constituted of porous polymer inclusions in a PDMS matrix. We demonstrated that the static mixing laws we obtained are also valid for propagating systems and, conse-quently, that they could be applied to the resonator designed in chapter 3.

Prospect

As prospects for future work, the first research work to do would be performing new simulations to exploit the mixing law in the 1D nonlinear elastic resonator presented in chapter 3. For this task, a study should be performed on the materials, to find the most appropriate constituents for the matrix and nonlinear inclusions of the central resonating cavity.

Then, this model of elastic resonator should be extended to 2D (an array of cylindrical inclusions in a matrix, infinite in the z-direction) and 3D (an array of cylindrical inclusions in a plate with finite thickness).

The later case is the closest to what would be an elastic wave switching device. When a convenient structure is proved by simulations to have the expected switching behavior, it could be actually created through micromachining, tested and characterized.

Finally, the use of the magnetoelastic effect to tune the nonlinear elastic constants of the material present in the defect should be studied. This would open the way for nonlinear dynamically tunable PnCs. The influence of a magnetic wave on the nonlinear elastic constants of materials through the magnetoelastic effect has already been studied and its use to tune the behavior of PnCs has been proposed. The study of the nonlinear elastic switch in chapter 3 is the first step for creating an magnetically controlled elastic switch.

The designed structure along with its analytical model and parametric study provide a good base for future studies that would include the magnetoacoustic effect.

Once the magnetoacoustic effect is successfully included in the model, the use of the spin reorientation transition effect would improve the tunability of the system, since its instability induces a large variation of the parameters from a small variation of the magnetic field. (A.17)

These expressions are finally rewritten as:

Û                                   
                                   =                                    ε Tr ε Î ε Tr ε Tr ε 2 Î Tr 2 ε Î ε Tr ε 2 Tr ε Tr ε 2 Î ε Tr 2 ε Tr 3 ε Î                                    , (A.18)
where Û is the upper-triangle matrix:

                                                         L ′ M ′ N ′ O ′ P ′ Q ′ R ′ S ′ T ′ 0 L' +2M' 0 2O ′ 2P ′ + N ′ 0 Q ′ + 2R ′ 0 S ′ + 2T ′ 0 0 L'(L' +2M') 0 M ′ (L ′ + 2M ′ ) 2O ′ L ′ 4O ′ M ′ + O ′ L ′ 2
L ′3 M ′ L ′2 2L ′ (L ′ M ′ + M ′2 ) 2M ′2 (L ′ + M ′ ) 0 0 0 0 0 0 L ′2 (L ′ + 2M ′ ) 0 2M'(L'+M') (L'+2M') 0 0 0 0 0 0 0 L ′ (L ′ + 2M ′ ) 2 M ′ (L ′ + 2M ′ ) 2 0 0 0 0 0 0 0 0 (L ′ + 2M ′ ) 3                                                          . (A.19)
We invert equation (A.18) to obtain εd , Tr( εd ), etc. and we inject this result in the expression of T : .20) to get:

T = c                                   2µ 2 -2µ 1 K 2 -µ 2 -K 1 + µ 1 2e e 3f 4g l l h                                   T                                  
                                   +                                   2µ 1 K 1 -µ 1 0 0 0 0 0 0 0                                   T                                    ε Tr ε Î ε Tr ε Tr ε 2 Î Tr 2 ε Î ε Tr ε 2 Tr ε Tr ε 2 Î ε Tr ε 2 Tr 3 ε Î                                    , ( A
T =                                                                      2µ 1 K 1 -µ 1 0 0 0 0 0 0 0                                   T + c                                   2µ 2 -2µ 1 K 2 -µ 2 -K 1 + µ 1 2e e 3f 4g l l h                                   T Û -1                                                                       ε Tr ε Î ε Tr ε Tr ε 2 Î Tr 2 ε Î ε Tr ε 2 Tr ε Tr ε 2 Î ε Tr ε 2 Tr 3 ε Î                                    , (A.21)
and therefore we can identify all the effective parameters as follows:

                                  2µ eff K eff -µ eff 2e eff e eff 3f eff 4g eff l eff l eff h eff                                   =                                   2µ 1 K 1 -µ 1 0 0 0 0 0 0 0                                   T + c                                   2µ 2 -2µ 1 K 2 -µ 2 -K 1 + µ 1 2e e 3f 4g l l h                                   T Û -1 (A.22)
The exact full expressions of the effective parameters are given by Maple. They are too big to be written here. For numerical simulations, the matrix Û is inverted numerically in Octave or with numpy in Python programs.

Derivation to the third order in 3D for the leading terms In 3D, for the sake of simplicity, we only consider the case of longitudinal compression waves.

So, we use the following expression for the stress-strain relation: which we can write in a matrix form:

            L ′ M ′ Q ′ W ′ 0 (L ′ + 3M ′ ) 3Q ′ 3W ′ 0 0 (L ′ + 3M ′ ) 2 6Q ′ (L ′ + 3M ′ ) 0 0 0 (L ′ + 2M ′ ) 3                        
εd Tr( εd ) Tr 2 ( εd ) Tr 3 ( εd )

            = Û            
εd Tr( εd ) Tr 2 ( εd ) Tr 3 ( εd )

            =             ε Tr ε Tr 2 ε Tr 3 ε             . (A.38)
The final expression for T is: Here again, the upper diagonal matrix Û is inverted to obtain T in terms of ε and of the properties of the matrix and the inclusions and concentration c. Finally, the effective parameters can be identified in:

T = c 2µ 2 εd + K 2 -
            2µ eff K eff -2 3 µ eff C eff 4H eff             T =             2µ 1 K 1 -2 3 µ 1 0 0             T + c             2µ 2 -2µ 1 K 2 -2 3 µ 2 -K 1 + 2 3 µ 1 C 4H             T Û -1 (A.40)
The exact full expressions of the effective parameters are given by Maple. They are too big to be written here. For numerical simulations, Û is inverted numerically in Octave or with numpy in Python programs.

List of Figures

Introduction 1 1

 1 Analytical and Numerical Study of 1D Superlattices 9 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2 One-Dimensional Superlattices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3 Analytical Method: Transfer Matrix Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3.1 Propagation in a Layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3.2 Propagation in a Bilayer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.3.3 Band Structure in a 1D Phononic Crystal . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.3.4 Propagation of Amplitudes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.3.5 Transmission Through a Bilayer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.3.6 Transmission Through N Bilayers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.3.7 Transmission Through N Bilayers with a Defect . . . . . . . . . . . . . . . . . . . . . 22 1.3.8 Reflected Impedance Through a Multilayer . . . . . . . . . . . . . . . . . . . . . . . 23 1.3.9 Transmission Through N Bilayers with Quarter-Wavelength Layers . . . . . . . . . . 24 1.4 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.4.1 Pseudospectral Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.4.2 Finite Difference Time Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 1.4.3 Spectral Energy Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 1.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 ix

Figure 1 -

 1 Figure 1 -Exemple de cristal phononique 2D. Cet exemple montre le cristal et ses paramètres de maille sur la gauche et une structure de bande sur la droite, ainsi que les directions de la zone de Brillouin.

Figure 2 -

 2 Figure 2 -Un super-réseau 1D avec des couches d'impédances Z 1 et Z 2 .Nous écrivons les équations fondamentales de l'acoustique :∂p ∂x = -ρ(x) ∂v ∂t et ∂p ∂t = -Y (x) ∂v ∂x (1) où x est la position, t le temps, p est la pression, v la vitesse particulaire, ρ la masse volumique Y le module d'Young. Pour une onde harmonique, la solution générique de ce système d'équations pour une couche d'impédance Z i est une onde plane, composée d'une onde incidente d'amplitude A et d'une onde réfléchie d'amplitude B. Dans ce cas, en considérant p = Pe jωt et v = V e jωt , nous pouvons écrire la solution ainsi :

Figure 3 -

 3 Figure 3 -Courbes de dispersion calculées avec la loi de dispersion. (a) Courbe de dispersion pour des inclusions d'hydroxy apatite dans du collagène d 1 = d 2 = 33.5 nm, Z 1 ≈ 13.8 × 10 6 kg m -2 s -1 et Z 2 ≈ 0.353 × 10 6 kg m -2 s -1 . (b) Courbe de dispersion d'inclusions d'acier dans de la silice d 1 = d 2 = 3.0665 µm, Z 1 ≈ 46.8 × 10 6 kg m -2 s -1 et Z 2 ≈ 13.1 × 10 6 kg m -2 s -1 .

Figure 4 -

 4 Figure 4 -Spectre en transmission pour une structure avec et sans défaut. Spectre en transmission pour cinq inclusions d'acier dans de la silice d 1 = d 2 = 3.0665 µm, Z 1 ≈ 46.8 × 10 6 kg m -2 s -1 et Z 2 ≈ 13.1 × 10 6 kg m -2 s -1 . Les milieux d'entrée et de sortie sont tous deux de la silice. (a) Structure parfaite. (b) Structure avec un défaut central en acier de largeur double : d C = 2 × d 1 .

Figure 5 -Tableau 1 -

 51 Figure 5 -Représentation schématique du modèle de la structure osseuse simulée.

  Pa) 9, 565 × 10 7 1, 543 × 10 11 -1, 571 × 10 13 6, 958 × 10 14 -1, 650 × 10 16

Figure 6 -

 6 Figure 6 -Graphique de densité d'énergie spectrale (DES) pour le vecteur d'ondes π/a. Ligne noire : avec le module d'Young tronqué à l'ordre zéro (Y constant, cas linéaire) ; tirets rouges : avec Y tronqué à l'ordre 2 ; points verts : avec Y tronqué à l'ordre 4. Les courbes avec Y tronqué à l'ordre 6 ou 8 ne sont pas représentées car elles se superposent presque parfaitement avec celle à l'ordre 4. Insert : agrandissement de la zone hachurée montrant le remplissage des bandes interdites entre 0 et 8 GHz.

Figure 7 -

 7 Figure 7 -Représentation du cristal phononique 1D considéré. (a) Structure « parfaite » (sans défaut). (b) Structure avec un défaut central. Zones rouges : inclusions d'acier ; zones bleues : matrice de SiO 2 .

Figure 8 -

 8 Figure 8 -« Fit » de la puissance de sortie par la courbe théorique. Croix noires : résultats numériques avec une source de vitesse croissante ; cercles noirs : résultats numériques avec une source de vitesse décroissante ; ligne bleue : courbe en « S » calculée à partir des paramètres de simulation ; ligne rouge : transmission complète (P out = P in ).

Figure 9 -

 9 Figure 9 -Un système masse-ressort comportant des ressorts de constantes K 1 et K 2 .

Figure 10 -

 10 Figure 10 -Inclusions non-linéaires dans une matrice linéaire. (a) En 3D, nous définissons les consantes non-linéaires quadratiques A, B, C et cubiques E, F, G, H. (b) En 2D, celles-ci se simplifient en e et f pour les quadratiques et g, h et l pour les cubiques.

  Les vitesses longitudinales et transverses et le coefficient de compressibilité pour ces matériaux sont données par le tableau 3. Les paramètres non-linéaires du polymère poreux sont : f = 10 × K 2D = 1.28 GPa et h = 100 × K 2D = 12.8 GPa. Les constantes effectives linéaires (µ eff et K eff ) et non-linéaires (f eff et h eff ) du matériau effectif obtenu par la procédure d'homogénéisation sont calculées pour un ensemble de porosités du polymère. La valeur de porosité 0.0087 est choisie car elle optimise les amplifications du paramètre cubique. Pour cette porosité, la concentration optimale en polymère dans le PDMS est 0.0067. Pour cette concentration, l'amplification cubique atteint 1042, 6 et l'amplification quadratique 4, 69, comme le montre la figure 11.

Figure 11 -

 11 Figure 11 -Variation des amplifications quadratique et cubique avec la concentration.

Figure 12 -

 12 Figure 12 -Validation linéaire et non-linéaire de l'opérateur élastodynamique. (a) Propagation d'une onde élastique dans un demi-espace linéaire, avec une perturbation produite par une source directionnelle ponctuelle dans Hedge comparée avec le résultat analytique. (b) Évolution de la valeur de la troisième harmonique avec la distance x parcourue depuis la source dans Hedge comparée avec le résultat analytique. Enfin, nous avons simulé la propagation d'une onde élastique à travers une structure hétérogène composée d'un matrice linéaire de PDMS et d'inclusions non-linéaire en polymère poreux, dont les constantes sont données dans la section précédente.Dans cette structure, représentée sur la figure13, la concentration de polymère poreux dans le PDMS est de 0,67 %, concentration censée encourager l'amplification des non-linéarités cubiques au détriment des non-linéarités quadratiques.En enregistrant l'évolution des deuxième et troisième harmoniques en fonction de la distance parcourue par l'onde depuis le point source, nous avons pu estimer la valeur des constantes non-linéaires d'ordre deux et trois grâce aux formules de la littérature précédemment utilisées pour la validation du code.Nous avons ainsi obtenu une amplification du paramètre quadratique par un facteur 4,7 alors que

Figure 13 -

 13 Figure 13 -La structure comporte 200 inclusions, avec un diamètre r = 0.04618 m et un espacement a = 1 m entre chaque centre. Une NPML est présente à chaque extrémité de la structure. Les positions des source et récepteur sont indiquées.

Figure 14 -

 14 Figure 14 -Example of a phononic crystal. This example shows the lattice parameters for the Bragg resonance on the left and the resulting band diagram on the right, as well as the directions of propagation in the Brillouin zone.

Figure 15 -

 15 Figure 15 -Basic diagram of a precompressed granular phononic crystal.

Figure 1 .

 1 Figure 1.1 -A 1D superlattice with alternate layers with impedances Z 1 and Z 2 .

Figure 1 . 2 -

 12 Figure 1.2 -Dispersion curves calculated with the dispersion law. (a) Dispersion curves of hydroxy apatite inclusions in collagen d 1 = d 2 = 33.5 nm, Z 1 ≈ 13.8 × 10 6 kg m -2 s -1 and Z 2 ≈ 0.353 × 10 6 kg m -2 s -1 . (b) Dispersion curves of steel inclusions in silica d 1 = d 2 = 3.0665 µm, Z 1 ≈ 46.8 × 10 6 kg m -2 s -1 and Z 2 ≈ 13.1 × 10 6 kg m -2 s -1 . Red dots: analytical solution for a = 1.

Figure 1 . 3 -

 13 Figure 1.3 -Graphical solutions of equation (1.47). The points corresponding to the solutions are shown in red. The bold blue lines between those points correspond to the gaps.

Figure 1 . 4 -

 14 Figure 1.4 -Graphical solutions of equation (1.68). The solutions correspond to the case where a = 2. The points corresponding to the solutions are shown in red.

Figure 1

 1 Figure 1.6.

Figure 1 . 7 -

 17 Figure 1.7 -Transmission curves for 5 and 7 bilayers. Transmission curves for (a) five / (b) seven steel inclusions in silica d 1 = d 2 = 3.0665 µm, Z 1 ≈ 46.8 × 10 6 kg m -2 s -1 and Z 2 ≈ 13.1 × 10 6 kg m -2 s -1 . The input and output media are both silica. A band gap appears between ~300 MHz and ~700 MHz.

Figure 1 . 8 -

 18 Figure 1.8 -A 1D system composed of N bilayers and a central defect with width d C and impedance Z C .

Figure 1 . 9 -

 19 Figure 1.9 -Transmission curves for 5 and 7 bilayers with a central defect. Transmission curves for (a) five / (b) seven steel inclusions in silica d 1 = d 2 = 3.0665 µm, Z 1 ≈ 46.8 × 10 6 kg m -2 s -1 and Z 2 ≈ 13.1 × 10 6 kg m -2 s -1 . The input and output media are both silica. The central defect is made of steel. Its width is d C = 2 × d 1 .

Figure 1 .

 1 Figure 1.10 -A 1D system composed of N bilayers with external layers with impedances Z C and Z L .

Material ρ (kg m - 3 )

 3 C 11 (GPa) Z = ρC 11 (kg m -2 s -1

Figure 1 . 11 -

 111 Figure 1.11 -Quality factor versus the number of bilayers for different materials. The solid lines are the result of equation (1.120). The black asterisks are the results of numerical simulations.

Figure 1 . 12 -

 112 Figure 1.12 -Staggered grid as used in the Pseudospectral method. The stress and particle velocity are defined on different grid points and different timesteps.

Figure 2 . 1 -

 21 Figure 2.1 -Young's modulus Y , stress σ and elastic energy E versus strain for collagen in the presence of water.

  Material Density ρ (kg m -3 ) Young's modulus Y (GPa) c = Y (0)/ρ (m s -1 )

  21 wave vectors ranging from 0 to π/a by step of 0.05 π/a. The values used for the initial displacement of the masses are ∆x/45 and ∆x/150. The wave number interval [0, π/a] constitutes the first Brillouin zone of the phononic periodic structure.

Figure 2 . 4 -

 24 Figure 2.4 -SED spectrum of the system at low and high initial displacements. Solid black lines: spectral energy density of the phononic crystal structure with an initial displacement of (a) ∆x/150 and (b) ∆x/45. Blue dashes: linear band diagram obtained with a matrix transfer method. The SED modes match almost perfectly the linear band diagram for low initial displacements. At high initial displacements, a large frequency shift can be observed between the SED modes and the linear modes.

Figure 2 . 5 -

 25 Figure2.5 -Taylor series of the elastic energy versus strain for the collagen/water open system. The solid line is the same as that in Fig.1, namely the integral of the Taylor series decomposition of the stress as a function of strain. The dotted and dashed lines are the terms of orders 2, 4, 6 and 8 of this Taylor series. The term of order 6 is particularly interesting as it is gives a concave down curve characteristic of an unstable system (its second derivative is negative).

Figure 3 .

 3 Figure 3.1 represents the resonance curve plotted for different values of χ and a. We chose ω res = 100 Hz.

  linear case. When a increases (a = 20 m, Figure 3.1c), the curve shows a critical point where the slope is infinite. For even larger values of a (a = 50 m, Figure 3.1d), the resonance curve shows an bistable zone and folds over itself. Some of the points in this bistable zone (curve between points C and D) have no physical meaning. When ξ increases from the left side of the curve, the output amplitude b will "jump" directly from point C to point E and when it decreases from the right side of the curve, it will jump directly from point D to point G, as the red arrows show.

Figure 3 . 1 -

 31 Figure 3.1 -Comparison of resonance curves. Resonance curve of a system with (a) χ = 0 and a = 50 m (linear case), (b) χ = 0.2 s 2 m -1 and a = 5 m (quasi-linear case), (c) χ = 0.2 s 2 m -1 and a = 20 m and (d) χ = 0.2 s 2 m -1 and a = 50 m.

Figure 3 . 2 -

 32 Figure 3.2 -Phononic crystal composed of N bilayers with a central defect.

Figure 3 . 3 -

 33 Figure 3.3 -Two views of the "S" curve. (a) Curve of P in vs. P out . (b) Curve of P out vs. P in with remarkable points A (full transmission) and B and C (jumps).

1 .Table 3 . 1 -

 131 Material ρ (kg m -3 ) C 11 (GPa) c L (m/s) Density, linear elastic constant and longitudinal elastic wave velocity for steel and silica.In this structure, represented on Figure3.4, the steel layers have quarter-wavelength thicknesses. With a basis frequency f 0 = 500 MHz, the wavelength in steel is c L /f 0 = 12.266 µm. Thus, the thicknesses of the steel layers are d 1 = 3.0665 µm. For the sake of simplicity, we chose the SiO 2 layers' thicknesses d 2 such as d 2 = d 1 .

Figure 3 . 4 -

 34 Figure 3.4 -Representation of the considered 1D PnC without (a) and with (b) a defect. Red: steel; blue: SiO 2 . In the defected structure, represented on Figure 3.4b, the central inclusion thickness d C is chosen twice of d 2 i.e. d C = 6.133 µm and constitutes a defect inserted inside the perfect PnC of Figure 3.4a. The space resolution is dz = 0.306 65 µm, which means that the thickness of each inclusion is 10 dz (resp. 20 dz for the defect).

Figure 3 . 5 -

 35 Figure 3.5 -Transmission spectrum of the structure. Blue: transmission spectrum of the perfect structure obtained with the transfer matrix method. Red: transmission spectrum of the defected structure obtained with the pseudo-spectral method. Inset: zoom on the defect peak, whose central frequency and quality factor are f res = 491.95 MHz and Q = 608.36, respectively.

Figure 3 . 8 -

 38 Figure 3.8 -Fit of the output power by the theoretical curve. Black crosses: numerical result with increasing source levels; black circles: numerical result with decreasing source levels; blue line: "S" curve with simulation parameters; red line: full transmission (P out = P in ). The parameters of the theoretical curve are Q = 608, ∆ = 4.57 and β = -4.7 × 10 31 m -2 s -2 .

Figure 3 . 9 -

 39 Figure 3.9 -Variation of the shape of the S curve according to the variation of β, ∆ and Q. (a) Variation according β with ∆ = 4.51 and Q = 600. (b) Variation according ∆ with β = -4.5 × 10 31 m -2 s -2 and Q = 600. (c) Variation according Q with ∆ = 4.51 and β = -4.5 × 10 31 m -2 s -2 . Straight line: full transmission (P in = P out ).

Figure 3 . 10 -

 310 Figure 3.10 -Variation of the threshold power and width of the "S" curve. (a) Variation according to β with ∆ = 4.51 and Q = 600. (b) Variation according to ∆ with β = -4.5 × 10 31 m -2 s -2 and Q = 600. (c) Variation according to Q with ∆ = 4.51 and β = -4.5 × 10 31 m -2 s -2 .

Figure 4 . 1 -

 41 Figure 4.1 -Mass-spring system model containing nonlinear springs with constants K 1 and K 2 .

  and δ m c opt (n = 3) The global nonlinearity of the system for each concentration of air is represented on Figure 4.2. The values observed in the curve correspond to the expected results.

Figure 4 . 2 -

 42 Figure 4.2 -Global nonlinear parameters β m and δ m of a medium composed of air bubbles in water.

(

  1951), written l, m and n. Alternatively, Taupin and Bernstein (1961) use the three constants c 1 , c 2 and c 3 . Bland (1969) uses another notation with constants α, β and γ. Finally, Eringen and Suhubi (1974) use a notation with constants l E , m E and n E . Tables 4.2 and 4.3 provide the relations between the components of the C IJK tensor, the Landau constants and the constants used in notations mentioned above.

Figure 4 .

 4 Figure 4.3 defines the linear and nonlinear constants for those materials.

Figure 4 . 3 -

 43 Figure 4.3 -Quadratic nonlinear inclusions in a linear matrix. (a) In 3D, we define the quadratic nonlinear parameters A, B and C. (b) In 2D, we define the quadratic nonlinear parameters e and f .

. 64 )

 64 Since Tr( Î ) = 2, the trace of the expression in equation (4.64) is:2 det( ε) -Tr 2 ( ε) + Tr( ε2 ) = 0. (4.65) So, by replacing det( ε) from equation (4.65) into equation (4.64), we get:

Figure 4 . 4 defines

 44 the linear and nonlinear constants for those materials.

Figure 4 . 4 -

 44 Figure 4.4 -Cubic nonlinear inclusions in a linear matrix. (a) In 3D, we define the quadratic A, B, C and cubic E, F, G, H nonlinear parameters. (b) In 2D, we define the quadratic e, f and cubic g, h, l nonlinear parameters.

(4. 109 ). 110 ) 2

 1091102 So, finally, the tensor T is written as:T = 2µ ε + λ Tr( ε) Î + A ε2 + BTr( ε2 ) Î + 2B ε Tr( ε) + C Tr 2 ( ε) Î + E Tr( ε3 ) Î + 3E ε2 Tr( ε) + 2F Tr( ε) Tr( ε2 ) Î + 2F ε Tr 2 ( ε) + 4G ε Tr( ε2 ) + 4H Tr 3 ( ε) Î . (4Restriction to 2D We replace expressions (4.66) and (4.67) of ε2 and Tr( ε3 ) respectively in equation (4.110) and get: T = 2µ ε + λ Tr( ε) Î + A ε Tr( ε) -1 Tr 2 ( ε) -Tr( ε2 ) Î + B Tr(ε 2 ) Î + 2B ε Tr( ε) + C Tr 2 ε) -Tr( ε2 ) Î Tr( ε) + 2F Tr( ε) Tr( ε2 ) Î + 2F ε Tr 2 ( ε) + 4G ε Tr( ε2 ) + 4H Tr 3 ( ε) Î .

  .113) or written as a tensor:T = 2µ ε + λ Tr( ε) Î + e Tr( ε2 ) Î + 2e ε Tr( ε) + 3f Tr 2 ( ε) Î + 4g ε Tr( ε2 ) + h Tr 3 ( ε) Î + l ε Tr 2 ( ε) + l Tr( ε) Tr( ε2 ) Î .

  ε∞ = L εd + M Tr( εd ) Î + N εd Tr( εd ) + O Tr( ε2 d ) Î + P Tr 2 ( εd ) Î + Q εd Tr( ε2 d ) + R Tr( εd ) Tr( ε2 d ) + S εd Tr 2 ( εd ) + T Tr 3 ( εd ) Î , (4.116)where L, M, N , O, P, Q, R, S, T depend on the first-, second-and third-order elastic constants of the inclusions and the linear elastic constants of the matrix. Their expressions are given in the appendix, in equation (A.10).

ε = c εd + ( 1 -

 1 c)ε ∞ . This gives us the following expression for the average strain ε :ε = (c + (1 -c)L) εd + (1 -c)M Tr( εd ) Î + (1 -c)N εd Tr( εd ) + (1 -c)O Tr( ε2 d ) Î + (1 -c)P Tr 2 ( εd ) Î + (1 -c)Q εd Tr( ε2 d ) + (1 -c)RTr( εd ) Tr( ε2 d ) + (1 -c)S εd Tr 2 ( εd ) + (1 -c)T Tr 3 ( εd ) Î .

Figure 4

 4

Figure 4 . 5 -

 45 Figure 4.5 -Variation of the quadratic and cubic amplifications with concentration. Variation of the quadratic and cubic amplifications vs. concentration of porous polymer in the composite system.

101

 101

5. 2 A

 2 tool for 2D Nonlinear Elastodynamics: Hedge

  Klöckner et al. has implemented a DG method in a software framework called Hedge. This enables nodal DG to run on multi-CPUs and multi-GPUs. The authors applied their algorithm to the resolution of Maxwell's equation with polynomial orders from 1 to 9. The speedup between CPUs and GPUs varies from 14 to 65, the maximum being reached for the order 4[START_REF] Klöckner | Nodal discontinuous Galerkin methods on graphics processors[END_REF].

Figure 5 .

 5 2 and 5.3, the curves represent the propagation of an elastic wave in a solid with density ρ 0 = 2200 kg m -3 , longitudinal velocity c l = 3200 m s -1 , transverse velocity c t = 1847.5 m s -1 with a Ricker wavelet source using parameters f c = 7.25 Hz, t D = 0.16 s and a 1 = -(πf c ) 2 in the equation:

Figure 5 . 1 -

 51 Figure 5.1 -Rectangular 2D structure used for the validation of the Elastodynamics Operator

Figure 5 . 2 -

 52 Figure 5.2 -Horizontal (a) and vertical (b) particle displacement recorded at the first receiver.

Figure 5 . 3 -

 53 Figure 5.3 -Horizontal (a) and vertical (b) particle displacement recorded at the second receiver.

Figure 5 . 4 -

 54 Figure 5.4 -Numerical error when simulating wave propagation in our system between orders 1 and 9.

  kg/m 3 959.58 m/s 109.65 m/s 41.39 GPa 91920 GPa Table

5. 3 . 1

 31 Validation of the Nonlinear ParametersIn order to validate the principle of quadratic and cubic nonlinear parameters measurements, we study the propagation of a longitudinal plane wave in a homogeneous medium. We consider the structure represented on Figure5.5 using the cubic Elastodynamics Operator in Hedge. The total length of the structure is 300 m and its width is 20 m. The structure is periodical respectively to the y-axis. Two NPML with a thickness of 20 m are positioned at each end of the structure. The positions of the line source and the 20 receivers are also depicted on Figure5.5.

Figure 5 . 5 -

 55 Figure 5.5 -Schematic of the model used for the validation of nonlinearities measurement in Hedge. The simulation domain consists in a 300 m × 20 m elongated rectangle made of the effective isotropic medium described in section 4.6.3. The properties of the effective medium are given in table 5.1. Two NPML, with a thickness of 20 m, are positioned at each end of the structure. The positions of the source and the 20 receivers are depicted.

Figure 5 . 7 -

 57 Figure 5.7 -3rd harmonic vs. distance in a homogeneous medium with cubic nonlinearities only. Evolution of the third harmonic component when the quadratic nonlinear parameter β = 0 as a function of distance x calculated with equation (5.24) (blue line) and with Hedge (stars).

Figure 5 . 8 -

 58 Figure 5.8 -3rd harmonic vs. distance in a homogeneous medium with quadratic and cubic nonlinearities. Evolution of the third harmonic component as a function of distance x calculated with equation (5.23) (blue line) and with Hedge (stars).

5. 3 . 2

 32 Validation of the Nonlinear Mixing Law in a Propagative SystemWe will now study the validity of the mixing law developed in chapter 4 when used in the prediction of the propagation of an elastic wave in an heterogeneous nonlinear medium. Thus, we study the structure represented on Figure5.9 using the cubic Elastodynamics Operator in Hedge. The total length of the structure is 200 m and its width is 1 m. The structure is composed of 200 unit cells with size 1 m × 1 m. The structure is periodical respectively to the y-axis. Two 50-m-wide NPMLs are located on each x side of the structure.

Figure 5 . 9 -

 59 Figure 5.9 -Heterogeneous medium constituted of nonlinear inclusions in a linear matrix. The structure contains 200 inclusions, with size r = 0.04618 m and a = 1 m spacing between each center. A NPML is present at each end of the structure. The positions of the source and the receiver are indicated.

Figure 5 . 10 -

 510 Figure 5.10 -Band structure of a PnC made of cylindrical inclusions of porous polymer in a PDMS matrix. Band structure of a PnC made of 0.04618 m radius inclusions of porous polymer in a PDMS matrix, and positioned in a square lattice with lattice parameter 1 m.

Figure 5 . 11 -

 511 Figure 5.11 -Particle velocity comparison in PnC and effective homogeneous medium. The particle velocity of a 20 Hz modulated sine wave is recorded for a propagation through a PnC (red) made of 0.04618 m radius inclusions of porous polymer in a PDMS matrix, and positioned on a square lattice with lattice parameter 1 m, and through an effective homogenous media (black), 90 m away from the source. Inset: zoom on the center part of the modulated sine wave.

3 .

 3 The evolutions of the second (a, b) and third harmonic (c, d) components as a function of distance x calculated with equations (5.16) or (5.23) (blue line) and with Hedge (stars) are compared for the heterogeneous structure of Figure5.9 and the effective homogeneous medium on Figure5.13. The evolutions calculated with Hedge are for both second and third harmonics in excellent agreement with the theoretical predictions.We use expressions (5.22) and (5.25) to evaluate the quadratic and cubic nonlinear parameters and amplifications and we compare them to the theoretical values f eff , h eff , A f and A h in table 5.2.All these results validate the possibility to use the mixing law developed in chapter 4 in the simulation of a nonlinear elastic medium up to third order.

Figure 5 . 12 -

 512 Figure 5.12 -Spectrum of the particle velocity of a 20 Hz modulated sine wave propagating through a PnC. Spectrum of the particle velocity of a wave propagating through a PnC made of 0.04618 m radius inclusions of porous polymer in a PDMS matrix, and positioned on a square lattice with lattice parameter 1 m, 90 m away from the source. Structure f eff (GPa) h eff (GPa) A f A h Heterogeneous 41.49 -94390.7 4.70 1070 Homogeneous effective 41.39 -91920.7 4.69 1042

Figure 5 . 13 -

 513 Figure 5.13 -2nd and 3rd harmonics vs. distance in homogeneous and heterogeneous NL media. Evolution of the second (a, b) and third harmonic (c, d) components as a function of the distance x calculated with equations (5.16) or (5.23) (blue line) and with Hedge (stars). (a) Second harmonic evolution in an heterogeneous system. (b) Second harmonic evolution in an effective homogeneous system. (c) Third harmonic evolution in an heterogeneous system. (b) Third harmonic evolution in an effective homogeneous system.

  We calculate the following expressions:Tr ε = (L ′ + 2M ′ ) Tr( εd ) + (2P ′ + N ′ ) Tr 2 ( εd ) + 2O ′ Tr( ε2 d ) + (Q ′ + 2R ′ ) Tr( εd ) Tr( ε2 d ) + (S ′ + 2T ′ ) Tr 3 ( εd ) ε Tr ε = L ′ (L ′ + 2M ′ ) εd Tr( εd ) + 2(N ′ L ′ + L ′ P ′ + N ′ M ′ ) εd Tr 2 ( εd ) + 2O ′ L ′ εd Tr( ε2 d ) + M ′ (L ′ + 2M ′ ) Tr 2 ( εd ) Î + (4P ′ M ′ + P ′ L ′ + M ′ N ′ ) Tr 3 ( εd ) Î + (4O ′ M ′ + O ′ L ′ ) Tr( εd ) Tr( ε2 d ) Î ε 2 = L ′2 ε2 d + 2L ′ M ′ εd Tr( εd ) + 2L ′ N ′ ε2 d Tr( εd ) + 2L ′ O ′ εd Tr( ε2 d ) + 2(L ′ P ′ + M ′ N ′ ) εd Tr 2 ( εd ) + M ′2 Tr 2 ( εd ) Î + 2M ′ O ′ Tr( εd ) Tr( ε2 d ) Î + 2M ′ P ′ Tr 3 ( εd ) Î Tr ε 2 = L ′2 Tr( ε2 d ) + 2M ′ (L ′ + M ′ ) Tr 2 ( εd ) + 2(L ′ N ′ + L ′ O ′ + 2M ′ O ′ ) Tr( εd ) Tr( ε2 d ) + 2(L ′ P ′ + M ′ N ′ + 2M ′ P ′ ) Tr 3 ( εd ) Tr 2 ε = (L ′ + 2M ′ ) 2 Tr 2 ( εd ) + 2(L ′ + 2M ′ )(2P ′ + N ′ ) Tr 3 ( εd ) + 4O ′ (L ′ + 2M ′ ) Tr( εd ) Tr( ε2 d ) ε Tr ε 2 = L ′3 εd Tr( ε2 d ) + 2L ′ (L ′ M ′ + M ′2 ) εd Tr 2 ( εd ) + M ′ L ′2 Tr( εd ) Tr( ε2 d ) Î + 2M ′2 (L ′ + M ′ ) Tr 3 ( εd ) Î ε Tr 2 ε = L ′ (L ′ + 2M ′ ) 2 εd Tr 2 ( εd ) + M ′ (L ′ + 2M ′ ) Tr 3 ( εd ) Î Tr ε Tr ε 2 = L ′2 (L ′ + 2M ′ )Tr( εd ) Tr( ε2 d ) + 2M ′ (L ′ + M ′ )(L ′ + 2M ′ ) Tr 3 ( εd ) Tr 3 ε = (L ′ + 2M ′ ) 3 Tr 3 ( εd ).

Td = 2µ 2 εd + λ 2

 2 Tr( εd )I + C Tr 2 ( εd )I + 4H Tr 3 ( εd )I , (A.23) where we have considered A = B = E = F = G = 0. This leads to: Tr( Td ) = 2µ 2 εd + Tr( εd )I + C Tr 2 ( εd )I + 4H Tr 3 ( εd )I . (A.24) and: T = c Td + C 1 [ εc εd ] . (A.36) We calculate the following expressions: ε = [c + (1 -c)L] εd + (1 -c)M Tr( εd )I + (1 -c)Q Tr 2 ( εd )I + (1 -c)W Tr 3 ( εd )I = L ′ εd + M ′ Tr( εd )I + Q ′ Tr 2 ( εd )I + W ′ Tr 3 ( εd )I Tr ε = (L ′ + 3M ′ ) Tr( εd ) + 3Q ′ Tr 2 ( εd ) + 3W ′ Tr 3 ( εd ) Tr 2 ε = (L ′ + 3M ′ ) 2 Tr 2 ( εd ) + 6Q ′ (L ′ + 3M ′ ) Tr 3 ( εd ) Tr 3 ε = (L ′ + 3M ′ ) 3 Tr 3 ( εd ), (A.37)

2µ 2 3

 3 Tr( εd )I + C Tr 2 ( εd )I + 4H Tr 3 ( εd )I+ 2µ 1 (L ′ -C) εd + M ′ Tr( εd )I + Q ′ Tr( εd )I + Q ′ Tr 2 ( εd )I + W ′ Tr 3 ( εd )I + K 1 -2µ 1 3 (L ′ -C + 3M ′ ) Tr( εd ) + 3Q ′ Tr 2 ( εd ) + 3W ′ Tr 3 ( εd ) I .(A.39)

4 1. 1 A 1 . 2 31 2. 1

 4112311 xxii 10 Inclusions non-linéaires dans une matrice linéaire. . . . . . . . . . . . . . . . . . . . . . . . xxiii 11 Variation des amplifications quadratique et cubique avec la concentration. . . . . . . . . . xxiv 12 Validation linéaire et non-linéaire de l'opérateur élastodynamique. . . . . . . . . . . . . . . xxv 13 Milieu hétérogène constitué d'inclusions non-linéaires dans une matrice linéaire. . . . . . xxvi 14 Example of a phononic crystal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 15 Basic diagram of a precompressed granular phononic crystal. . . . . . . . . . . . . . . . . . 1D superlattice with alternate layers with impedances Z 1 and Z 2 . . . . . . . . . . . . . . 9 Dispersion curves calculated with the dispersion law. . . . . . . . . . . . . . . . . . . . . . . 12 1.3 Graphical solutions of equation (1.47). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.4 Graphical solutions of equation (1.68). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 1.5 A 1D bilayer with impedances Z 1 and Z 2 with external layers with impedances Z 0 and Z L . 21 1.6 A 1D system composed of N bilayers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.7 Transmission curves for 5 and 7 bilayers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 1.8 A 1D system composed of N bilayers and a central defect with width d C and impedance Z C . 22 1.9 Transmission curves for 5 and 7 bilayers with a central defect. . . . . . . . . . . . . . . . . . 23 1.10 A 1D system composed of N bilayers with external layers with impedances Z C and Z L . . . 23 1.11 Quality factor versus the number of bilayers for different materials. . . . . . . . . . . . . . 27 1.12 Staggered grid as used in the Pseudospectral method. . . . . . . . . . . . . . . . . . . . . . 29 1.13 Mass-spring system as it is used by the SED method. . . . . . . . . . . . . . . . . . . . . . . Elastic curves for collagen in the presence of water. . . . . . . . . . . . . . . . . . . . . . . . 40 2.2 Schematic representation of the system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.3 Band structure of the AB superlattice. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 2.4 SED spectrum of the system at low and high initial displacements. . . . . . . . . . . . . . . 45 2.5 Taylor series of the elastic energy versus strain for the collagen/water open system. . . . . 46 2.6 SED spectrum for a specific wave vector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 139

  Tableau 3 -Vitesses longitudinale et transverse et coefficient de compressibilité du polymère et du PDMS.

					GPa)
	PDMS	970	1100	110	1,1737
	Polymère plein	1000	2000	1000	3,0
	Polymère poreux	991,3	131,2	65,6	0,128

  (Pa) 9.565 × 10 7 1.543 × 10 11 -1.571 × 10 13 6.958 × 10 14 -1.650 × 10 16

			a 2 4	ε 4 +	a 4 6	ε 6 +	a 6 8	ε 8 .	(2.25)
	The coefficients used in the Taylor's expansion series of the Young's modulus and the elastic energy are
	given in table 2.1. They were obtained by fitting experimental data from [13].
	i	0	2	4				6	8
	a i

Table 2 .

 2

1 -Coefficients used in Taylor's expansion series of Young's modulus and elastic energy.

Table 3 . 2

 32

	Material C 111 (GPa) C 1111 (GPa)	Γ	δ
	Steel	-3000	-30000	(3.73) 27.37

3.2

. For this study, we consider a purely cubic nonlinear steel, i.e. we take Γ = 0. -Nonlinear elastic constants of steel.

Table 4 .

 4 2 -Comparison between the different notations for the nonlinear quadratic constants.

	C IJK	Versus A, B, C Versus l, m, n	Versus C IJK
	C 111 , C 222 , C 333 C 112 , C 113 , C 122 , C 133 , C 223 , C 233 C 123 C 144 , C 255 , C 366 C 155 , C 166 , C 244 , C 266 , C 344 , C 355 C 456	2A + 6B + 2C 2B + 2C 2C B A 2 + B A 4	2l + 4m 2l 2l -2m + n m -n 2 m n 4	1 2 (C 112 -C 123 ) 1 4 (C 111 -C 112 ) 8 (C 111 -3C 112 + 2C 123 ) 1

Table 4 .

 4 3 -Summary of the relations between C IJK tensor, the Landau constants and l, m, n notation.

	which gives:

  .125) It is easy to observe that we obtain the larger values of A f for small values of x and large values of y. It can be seen by developing the function for x → 0 and y → ∞:

	A f =	y 4	-	xy 2 2	+	xy 2	+	x 2	+	1 4	+	x 2y	+	1 4y	+ ...,	(4.126)

table 4 .

 4

	8.

Table 4 .

 4 7 -Index of references for the elastic constants of some materials.

	Material	ρ 0 (kg m -3 ) (m s -1 ) (m s -1 ) (GPa) c l c t K 2D
	PDMS	970	1100	110	1.1737
	Porous polymer	991.3	131.2	65.6	0.128

Table 4 .

 4 8 -Density, speed and compressibility constants for PDMS and an arbitrary polymer.

  The difficult part in implementing NPMLs in Hedge is due to the need for a new partial differential equation, since NPML use the stretched fluxes F ′′ .

	The state equation, initially:
				∂Q ∂t	-	∂F x ∂x	-	∂F y ∂y	-	∂F z ∂z	= S(t),	(5.10)
	is transformed into:	∂Q ∂t	-	∂ ∂x

Table 5 . 2 -

 52 Comparison of observed and expected nonlinear parameters and amplifications. f eff , h eff , A f and A h are obtained with the analytical model.

1

  Exemple de cristal phononique 2D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 2 Un super-réseau 1D avec des couches d'impédances Z 1 et Z 2 . . . . . . . . . . . . . . . . . . xv 3 Courbes de dispersion calculées avec la loi de dispersion. . . . . . . . . . . . . . . . . . . . xvi 4 Spectre en transmission pour une structure avec et sans défaut. . . . . . . . . . . . . . . . . xvi 5 Représentation schématique du modèle de la structure osseuse simulée. . . . . . . . . . . . xvii 6 Graphique de densité d'énergie spectrale (DES) pour le vecteur d'ondes π/a. . . . . . . . . xix 7 Représentation du cristal phononique 1D considéré. . . . . . . . . . . . . . . . . . . . . . . xx 8 « Fit » de la puissance de sortie par la courbe théorique. . . . . . . . . . . . . . . . . . . . . xxi 9 Un système masse-ressort comportant des ressorts de constantes K 1 et K 2 . . . . . . . . . . .

ou 8, le comportement du système ne change pas significativement, ce qui prouve que c'est le terme d'ordre quatre de l'énergie qui joue le plus dans ce comportement. Nous pouvons en déduire que le comportement non-linéaire de la structure osseuse est associée à une fonction d'ordre quatre de l'énergie élastique et donc aux interactions quatre-ondes (quatre-phonons). En effet, si nous considérons que la déformation est une superposition d'ondes planes avec différents vecteurs d'ondes et fréquences, alors l'expression de ε 4 conduit à des interactions quatre-phonons. Ces interactions de diffusion conservent le moment et la fréquence et incluent différents processus comme la séparation d'un phonon en trois autres, la diffusion de deux phonons en en formant deux autres, etc.[START_REF] Manktelow | Multiple scales analysis of wave-wave interactions in a cubically nonlinear monoatomic chain[END_REF][START_REF] Swinteck | Phonon Scattering in One-Dimensional Anharmonic Crystals and Superlattices: Analytical and Numerical Study[END_REF]. Ces nouvelles interactions de phonons offrent davantage de canaux pour la dispersion de l'énergie mécanique.Nous pouvons comparer le comportement de cet super-réseau non-linéaire au comportement d'un système masse-ressorts à puits d'énergie multiples. Un système à puits d'énergie multiples est un système pour lequel la représentation de l'énergie élastique en fonction de la déformation n'est pas parabolique mais est une superposition d'une infinité de paraboles. À chaque parabole correspond une seule valeur de la constante de ressorts. La représentation d'une fonction d'énergie non quadratique par une multiplicité de puits paraboliques est équivalente à l'introduction d'un continuum de valeurs de constantes de ressorts. Aux hautes amplitudes, les ondes vont explorer des zones plus grandes de valeurs de déformation, visitant ainsi une multiplicité de puits d'énergie. Le diagramme de bandes montrera une continuité de modes au-dessus des fréquences primaires du système linéaire. Ces bandes passantes s'élargissent et remplissent finalement les bandes interdites. Les bandes passantes de ce système seront donc augmentées.Dans le cas présent, la raison pour laquelle les constantes de ressorts augmentent en fonction de l'amplitude de la déformation est que le modèle choisi pour le module d'Young est symétrique et son développement en série de Taylor ne contient que des puissances paires de ε.

L'utilisation de l'effet magnéto-élastique est ensuite envisagée pour permettre l'adaptation dynamique du système : cet effet permet de changer les constantes élastiques non-linéaires des matériaux et pourrait donc conduire à la réalisation d'un commutateur élastique contrôlé à distance par une onde magnétique. Enfin, l'instabilité de la transition de réorientation de spin pourrait encore améliorer les performances du système en permettant une variation importante des paramètres à partir d'une faible variation du champ magnétique.

Acknowledgements

thank my father Philippe, my mother Isabelle, my sisters Sylvie and Maylis, my wife Pauline and all my family-in-law for their patience and their support during all my life and particularly during the three years of my PhD. I acknowledge partial financial support from Centre National de la Recherche Scientifique I acknowledge partial financial support from the Bourse Mobilité Recherche from Conseil

governed by a similar expression as the one derived for fluids [START_REF] Adler | Determination of the Nonlinearity Parameter B/A for Water and m-Xylene[END_REF]:

xP 2 1 (0), (5.16) where P 1 (0) is the pressure at the source, and β = 1 + B l / (2A l ) is the quadratic nonlinear coefficient of fluid.

It appears that the amplitude of the second harmonic component increases linearly with the distance x, as it was obtained using Hedge. To be more quantitative in this comparison, we need to find the link between the nonlinear parameters defined in fluids, i.e. β and η, and the one introduced in the numerical simulation, i.e.

f and h. Following the procedure developed by Kostek et al. [START_REF] Kostek | Third-order elastic constants for an inviscid fluid[END_REF], and extended up to third order by Hamilton et al. [START_REF] Hamilton | Separation of compressibility and shear deformation in the elastic energy density (L)[END_REF], it is possible to link the nonlinear Landau's parameters to the linear and nonlinear parameters A l , B l , and C l of fluids: µ = 0,

(5.17)

Appendix

A.1 Derivation of the Landau Coefficients to the Third Order Derivation to the third order in 2D for all terms

If we write Ŝ the Eshelby's tensor and εd the strain inside the nonlinear inclusions 2, in a linear matrix 1, then, according to the Eshelby's theory of elasticity, we can write the two following results [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF]83]:

and:

where

The stress-strain relation inside the linear matrix is written:

We transform the previous relation to express ε1 as a function of T1 :

We replace Ŝ εd and Ŝ Ĉ-1 1 Td in equation (A.1) using equation (A.2):

119

We now replace Ĉ-1 1 Td and Tr Ĉ-1 1 Td using expression (A.4): with:

(A.10)

A.1. Derivation of the Landau Coefficients to the Third Order 121

We separate the 2D space V into two domains: the matrix V m and the inclusions V e . We define the ratio between the surfaces of the inclusions and the total surface:

We write the exact expression of the stress in the inhomogeneous medium:

This is an exact expression for which we had to make no approximation. Now, because c ≪ 1, we make the approximation that the stress and strain are homogeneous inside the inclusions, i.e. Td ≈ Td and εd ≈ εd . So:

and we use the following approximation:

Using equation (A.9) and expression (A.14), we obtain:

with:

The compression modulus is:

with the condition:

The Eshelby tensor is:

and the Eshelby equation is:

εd -S εd + SC -1 1 Td = ε∞ where Td = C 2 ( εd ) εd , (A. [START_REF] Cabaret | Amplitude-dependent phononic processes in a diatomic granular chain in the weakly nonlinear regime[END_REF] where:

and:

By writing the Eshelby tensor (A.27) in the Eshelby equation (A.28), we obtain: εd -6 5

which, when we replace C -1 1 by its expression, reduces to:

Tr( εd )I + 3 5 with:

5(3K 1 + 4µ 1 )

List of Tables

Theoretical and numerical study of nonlinear phononic crystals

Abstract

This work is dedicated to the theoretical and numerical study of nonlinear phononic crystals. The studied nonlinearities are those due to the second (quadratic) and third (cubic) order elastic constants of the materials that constitute the crystals. Nonlinear effects are studied by the means of finite element methods, used to simulate the propagation of an elastic wave through the crystals.

A first research project concerns the study of a bone structure, namely the dispersion of elastic waves in a structure composed of collagen and hydroxy apatite alternate constituent layers. Simulations showed that it exists a strong link between bones hydration and their ability to dissipate the energy.

The second study relates to an elastic resonator. A structure composed of steel inclusions in a silica matrix shows a switch behavior when the cubic nonlinearities of steel are taken into account. This strong nonlinear effect appears when the amplitude of the incident wave reaches a threshold. A full analytical model is provided.

The last study demonstrates the design of composite materials with both strong cubic nonlinearities and weak quadratic nonlinearities. The derivation of the mixing laws of the elastic parameters of a nonlinear material inside a linear one is performed up to order three. Equations show a strong amplification of the nonlinear parameters of the material for some concentrations. Numerical simulations allow to conclude that the above mentioned resonator can be produced.

For this thesis, an innovative tool based on the Discontinuous Galerkin (DG) finite element method is developed for the simulation of elastic wave propagation, in linear and nonlinear systems and in finite and semi-infinite media. The implementation of this DG code for 2D and 3D simulations benefits from the efficient exploitation of modern computer infrastructure (GPU units, clusters) using the property of massive parallelization of DG algorithms. This thesis is part of a joint agreement for an international PhD degree between École Centrale de Lille and the Materials Science and Engineering department of the University of Arizona at Tucson.

Keywords