# Ondes locales dans les milieux hétérogènes

## Aspects numériques

Destuynder Philippe
Wilk Olivier

Département Mathématiques
CNAM
Paris
France
destuynd@cnam.fr
wilk@cnam.fr

**RÉSUMÉ.** Lorsqu'un milieu présente une vitesse d'onde inférieure à celles de ceux qui l'entourent, un mécanisme de localisation d'énergie vibratoire peut apparaître. Cela provient du fait qu'une partie plus importante de l'énergie est réfléchie du côté du milieu le plus souple. Or la transition conduit à des contraintes d'interface très fortes qui sont susceptibles de créer un endommagement local. Notre ambition est de proposer une stratégie d'évaluation de l'énergie de ces surcontraintes sous forme d'un taux de restitution dynamique de l'énergie totale du système qui ne nécessite pas un calcul très précis au voisinage de l'interface.

**ABSTRACT.** At the interface between two media with different wave velocities, local stationary waves can appear. The mechanical explanation is that the reflection of the waves is more important inwards the softest media. The overstressing which appears at the interface between the two media can be at the origin of a damage mechanism. Our goal is to suggest an energetical method in order to be able to compute accurately the energy release rate due to this overstressing without being obliged to use a very refined mesh in the neighbourhood of the interface.

**MOTS-CLÉS** : Ondes locales, Taux de restitution de l'énergie, Endommagement, Méthodes numériques

**KEYWORDS** : Local waves, Energy release rate, Damage mechanisms, Numerical methods
# 1. Introduction

Considérons l'équation des ondes posée sur un ouvert $\Omega$ dont la frontière est noté $\Gamma$. On désigne par $\Omega_+$ et $\Omega_-$ deux sous-ouverts complémentaires dans $\Omega$ où les vitesses d'onde sont respectivement $c_+$ et $c_-$ (voir figure 1). Supposons pour simplifier que les

<figure><img src="image_1.png" /><figcaption>**Figure 1.** L'ouvert $\Omega$ sur lequel est posé l'équation des ondes</figcaption></figure>

conditions aux limites soient du type Dirichlet homogènes. Les conditions initiales étant notées $u_0$ et $u_1$. Le modèle consiste alors à trouver $y$ solution de :

$$ egin{cases} 
\ddot{u} - \operatorname{div}(c^2 \nabla u) = 0 \quad \text{dans} \ \Omega \times ]0, T[, \\ 
u(x, 0) = u_0(x), \quad \dot{u}(x, 0) = u_1(x) \quad \text{dans} \ \Omega, \\ 
u = 0 \quad \text{sur} \ \Gamma \times ]0, T[. 
\end{cases} 
\quad (1) $$

L'existence et l'unicité d'une solution sous des hypothèses raisonnables sont classiques et nous avons par exemple le résultat suivant.

**Théorème 1** Soient $u_0 \in H_0^1(\Omega)$ et $u_1 \in L^2(\Omega)$. Le système (1) admet une solution unique dans l'espace $C^1([0, T]; L^2(\Omega)) \cap C^0([0, T]; H_0^1(\Omega))$. En outre si on note $\{w_n\}_{n \ge 1}$ les modes propres solution de

$$ egin{cases} 
-\operatorname{div}(c^2 \nabla w) = \lambda w \quad \text{dans} \ \Omega, \\ 
w = 0 \quad \text{sur} \ \Gamma \quad \text{et} \quad \int_{\Omega} |w(x)|^2 dx = 1, 
\end{cases} 
\quad (2) $$

alors cette solution s'écrit :

$$ u(x, t) \sum_{n \ge 1} \alpha_n(t) w_n(x), $$
avec :

$$
\alpha_n(t) = \left[ \int_{\Omega} u_0(x)w_n(x)dx \right] \cos(\sqrt{\lambda_n}t) + \left[ \frac{\int_{\Omega} u_1(x)w_n(x)dx}{\sqrt{\lambda_n}} \right] \sin(\sqrt{\lambda_n}t).
$$

En fait, nous allons voir que l’ensemble des modes propres peut être séparé en deux sous-espaces de $H_0^1(\Omega)$. L’un constitué de vecteurs propres dits locaux que nous noterons $V^L$ et l’autre constitué de vecteurs dits intérieurs et que nous notons $V^I$. Nous verrons comment définir cette décomposition. Remarquons cependant que grâce à l’orthogonalité des modes propres ces deux sous-espaces sont bien orthogonaux à la fois au sens du produit scalaire de l’espace $L^2(\Omega)$ ainsi que celui de $H_0^1(\Omega)$. La décomposition de la solution sur ces deux sous-espaces nous amène à poser :

$$
u(x,t) = u^I(x,t) + u^L(x,t), \text{ où } \forall t \ge 0: \ u^I \in V^I \text{ et } u^L \in V^L. \quad (3)
$$

Si nous notons *E* l'énergie totale de la solution *u* (qui est constante vis-à-vis de *t*), et en
utilisant l'orthogonalité précédente, nous obtenons :

$$
\forall t \ge 0, \quad E(u) = \frac{1}{2} \int_{\Omega} |\dot{u}(x,t)|^2 dx + \frac{1}{2} \int_{\Omega} c^2 |\nabla u(x,t)|^2 dx = E(u^I) + E(u^L). \quad (4)
$$

Bien entendu, on pourrait envisager de calculer $E(u^L)$ en utilisant la décomposition de $u$ sur les modes propres locaux. Mais cela conduit à des calculs prohibitifs et dont la précision est discutable. C'est pour cela que nous proposons une méthode énergétique du type $\theta$-méthode [DD 81], pour évaluer l'énergie des ondes locales qui peuvent être responsables d'un endommagement de l'interface entre les deux milieux. Nous validons ensuite cette méthode sur quelques exemples numériques.

2. Un problème modèle

Considérons un ouvert rectangulaire comme celui de la figure 1 où le matériaux mou
occupe la portion pour laquelle $-a < x_2 < 0$. Les notations sont précisées sur cette
figure.

Compte tenu de la géométrie particulière de l'ouvert on peut rechercher les modes
propres solution de (2), sous forme de fonctions à variables séparées en x₁ et x₂. Pour
cela nous posons a priori :

$$
w(x_1, x_2) = \sqrt{\frac{2}{L}} \sin\left(\frac{n\pi x_1}{L}\right) q^n(x_2),
$$

et en reportant dans (2) nous obtenons que $q^n$ doit être solution de :

$$
\left\{
\begin{array}{l}
[\lambda - c_+^2(\frac{n\pi}{L})^2]q^n - c_+^2\frac{d^2q^n}{dx_2^2} = 0 \quad \text{pour } 0 < x_2 < H, \\
\\
[\lambda - c_-^2(\frac{n\pi}{L})^2]q^n - c_-^2\frac{d^2q^n}{dx_2^2} = 0 \quad \text{pour } -a < x_2 < 0, \\
\\
q^n(0^+) = q^n(0^-), \quad c_+^2\frac{dq^n}{dx_2}(0^+) = c_-^2\frac{dq^n}{dx_2}(0^-).
\end{array}
\right.
\qquad (5)
$$
<figure><img src="image_2.png" /><figcaption>**Figure 2.** L'ouvert du problème modèle</figcaption></figure>

On trouve deux types de solutions suivant la valeur du nombre $\xi = \frac{\sqrt{\lambda L}}{n\pi}$. Si $c_- < \xi < c_+$ alors les solutions sont des ondes locales qui présentent une décroissance exponentielle dans l'ouvert $\Omega_+$ depuis l'interface $\Gamma_i$. Les valeurs de $\lambda$ (pulsation propres) et donc de $\xi$ sont obtenues en résolvant l'équation transcendante :

$$ \frac{tg(J_{-}^{n}a)}{c_{-}^{2}J_{-}^{n}a} = -\frac{H}{a}\frac{th(K_{+}^{n}H)}{c_{+}^{2}K_{+}^{n}H} \quad (6) $$

où nous avons posé :

$$ [K_{+}^{n}]^{2} = c_{+}^{2} \frac{n\pi^{2}}{L} - \lambda \quad \text{et} \quad [J_{-}^{n}]^{2} = \lambda - c_{-}^{2} \frac{n\pi^{2}}{L}. $$

Les solutions sont de la forme suivante :

$$ \begin{cases} q^n(x_2) = A_+^n \mathrm{sh}(K_+^n(x_2 - H)) & \text{pour } x_2 \ge 0, \\ q^n(x_2) = A_-^n \mathrm{sin}(J_-^n(x_2 + a)) & \text{pour } x_2 \le 0. \end{cases} \quad (7) $$

Si par contre $\xi > c_+$ les solutions sont du type ondes intérieures et on les obtient en résolvant l'équation transcendante :

$$ \frac{tg(aJ_{-}^{n})}{aJ_{-}^{n}c_{-}^{2}} = -\frac{H}{a}\frac{tg(HJ_{+}^{n})}{HJ_{+}^{n}c_{+}^{2}} \quad (8) $$

où cette fois nous avons posé :

$$ [J_{+}^{n}]^{2} = \lambda - c_{+}^{2} \frac{n\pi^{2}}{L} \quad \text{et} \quad [J_{-}^{n}]^{2} = \lambda - c_{-}^{2} \frac{n\pi^{2}}{L}. $$

Dans ce cas les solutions sont :

$$ \begin{cases} q^n(x_2) = B_+^n \mathrm{sin}(J_+^n(x_2 - H)) & \text{pour } x_2 \ge 0, \\ q^n(x_2) = B_-^n \mathrm{sin}(J_-^n(x_2 + a)) & \text{pour } x_2 \le 0. \end{cases} \quad (9) $$
La discussion des solutions peut se faire aisément en utilisant une représentation graphique [DES 07]. On peut aisément compter les solutions locales en fonction de n. Nous avons par exemple le résultat suivant.

**Théorème 2** Il y a p solutions correspondant à des ondes locales pour une valeur donnée de n si

$$n \geq \frac{(2p + 1)L}{2a} \frac{1}{\sqrt{\left(\frac{c_+}{c_-}\right)^2 - 1}}.$$

Nous avons représenté quelques solutions locales et intérieures sur la figure 3.

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Ondes intérieures (à gauche) et locales (à droite) dans le cas de l'ouvert modèle</figcaption></figure>

### 3. Généralisation

Plaçons-nous maintenant dans le cas où l'ouvert $\Omega$ est celui de la figure 4. On utilise alors une carte locale pour construire une base de fonctions (locales) permettant de faire la synthèse de la solution $y$ de l'équation des ondes. En utilisant un système de coordonnées locales $(s, \zeta)$ (voir figure 4), nous pouvons réécrire le problème de valeurs propres (2) mais posé sur un voisinage tubulaire $\mathcal{T}$ de l'interface $\Gamma_i$ en notant $R(s)$ le rayon de courbure :

$$\left\{ \begin{array}{l} -\Delta_{s\xi} w - \frac{\xi}{R^2(1 + \frac{\xi}{R})} \frac{\partial R}{\partial s} \frac{\partial w}{\partial s} + \frac{1}{R(1 + \frac{\xi}{R})} \frac{\partial w}{\partial \xi} = \lambda w \text{ dans } \mathcal{T}, \\ w = 0 \text{ sur } \partial \mathcal{T}. \end{array} \right. \qquad (10)$$

Dans la mesure où la courbure $C = 1/R$ de $\Gamma_i$ ainsi que sa variation sont faibles, ce problème de valeurs propres est une perturbation du modèle modèle que nous avons discuté dans la section 1. En appliquant les résultats classiques des méthodes de perturbations compactes des modèles de valeurs propres (voir K.O. Friedrichs [FRI 61] et T. Kato
<figure><img src="image_4.png" /><figcaption>**Figure 4.** *Cas d'un ouvert arbitraire et passage en coordonnées locales*</figcaption></figure>

[KAT 66]), on peut conclure que les modes locaux et intérieurs trouvés à la section 1 sont génériques.

**4. Le taux de restitution dynamique**

Une fois connue, la base des modes propres, une possibilité théorique pour évaluer l'énergie qui leur est associée est de calculer les produits scalaires entre les conditions initiales et les modes en question. En ne retenant que ceux qui sont locaux, on obtient l'énergie qui sera *a priori* consacrée à l'endommagement de l'interface. Mais cette démarche n'est pas souhaitable car les modes locaux sont distribués tout le long de l'interface $\Gamma_i$ et ne permettent pas de retrouver une notion de force thermodynamique spécifique (attachée à un point matériel) et surtout, le nombre de modes locaux étant infini, leur calcul est illusoire et même si l'on dispose d'un critère de troncature, il semble peu réaliste d'envisager une telle approche. Pour ces raisons, nous proposons une approche énergétique globale du type *taux de restitution dynamique* dont le calcul est très stable et précis et qui ne nécessite pas une investigation locale et coûteuse de la solution de l'équation des ondes au voisinage de l'interface $\Gamma_i$. L'origine de cette méthode est la formulation thermodynamique des milieux continus. Elle est devenue courante dans les problèmes d'endommagement. On introduit une énergie par rapport à des paramètres observables, comme la position de la ligne de discontinuité $\Gamma_i$ et on calcule la dérivée de cette énergie par rapport à un mouvement virtuel de la variable. Le résultat répond au critère de l'état local si cette dérivée ne dépend que des différentes variables d'état locales. Mais dans le cas d'un système dynamique il est souvent préférable d'utiliser les équations de Lagrange pour construire les forces thermodynamiques. Imaginons un mouvement virtuel de la ligne $\Gamma_i$ qui est représenté par un champ de vecteur $\theta$. On prolonge ce champ à l'intérieur de l'ouvert $\Omega$ sur lequel est posé l'équation des ondes et on le note de la même
façon dans la mesure où aucune confusion n’est possible. Considérons ensuite l’action S du système mécanique qui est définie par :

$$S(\Gamma_i) = \frac{1}{2T} \int_0^T \int_{\Omega} |\dot{u}(x, t)|^2 dxdt - c^2 |\nabla u(x, t)|^2]dxdt, \quad (11)$$

et le dérivée de S par rapport à $\Gamma_i$ dans la direction $\theta$ est définie par :

$$G^{\theta} \equiv \lim_{\eta \to 0} \frac{S(\Gamma_i + \eta\theta) - S(\Gamma_i)}{\eta}. \quad (12)$$

Dans cette écriture, on note $\Gamma_i + \eta\theta$ la frontière $\Gamma_i$ déplacée par le champ de vecteurs $\eta\theta$. Un simple calcul mais un peu long pour être détaillé ici, nous permet d’écrire (des calculs similaires mais en régime permanent, sont effectués dans [DES 07] et [DES 06]) :

$$G^{\theta} = \frac{1}{2T} \int_{0}^{T} \int_{\Omega} [|\dot{u}|^{2} + c^{2} |\nabla u|^{2}] \mathrm{div}(\theta) - 2(D\theta \nabla u \cdot \nabla u) + \frac{1}{T} \left[ \int_{\Omega} \dot{u} \nabla u \cdot \theta \right]_{0}^{T} \quad (13)$$

où par exemple $\nabla u \cdot \theta$ est ici le produit scalaire entre les vecteurs $\nabla u$ et $\theta$ dans $\mathbb{R}^2$ et $D\theta$ est la matrice Jacobienne du champ $\theta$. En appliquant la formule de Stokes [RAV 83] on obtient une autre écriture de $G^\theta$ que nous avons résumé dans l’énoncé suivant.

**Théorème 3** Soit u la solution de l’équation des ondes (1) et $\theta$ un champ de vecteurs arbitraire sur $\Omega$ mais dont les composantes sont au moins éléments de l’espace $W^{1,\infty}(\Omega)$. Notons $\nu$ la normale unitaire à la frontière $\Gamma_0$ de l’ouvert $\Omega$ et ainsi que la normale unitaire à la ligne $\Gamma_i$ dirigée de $\Omega_-$ vers $\Omega_+$. On a alors :

$$G^{\theta} = \Lambda^{\theta} = \frac{c^2}{2T} \int_{0}^{T} \int_{\Gamma_{0}} \left| \frac{\partial u}{\partial \nu} \right|^2 \theta.\nu + \frac{[c_+^2 - c_-^2]}{2T} \int_{0}^{T} \int_{\Gamma_i} \left[ \frac{c^4}{c_+^2 c_-^2} \left| \frac{\partial u}{\partial \nu} \right|^2 + \left| \frac{\partial u}{\partial s} \right|^2 \right] \theta.\nu$$

s désignant l'abscisse curviligne le long de $\Gamma_i$ et on notera que $c^4|\frac{\partial u}{\partial \nu}|^2$ est continue à la traversée de $\Gamma_i$.

Le résultat précédent nous apporte trois informations très précieuses et sur lesquelles re-
pose la méthode que nous proposons.

1. La première est que $\Lambda^\theta$ d’après son expression même ne dépend que de $\theta.\nu$ sur $\Gamma_0$ et sur $\Gamma_i$. Il en donc de de même de $G^\theta$. En d’autre termes, peu importe le relèvement de $\theta.\nu$ que nous choisissions à l’intérieur de l’ouvert $\Omega$, le taux de restitution de l’énergie dynamique $G^\theta$ ne dépend que de $\theta.\nu$ sur $\Gamma_0 \cup \Gamma_i$.

2. La seconde information concerne l'aspect énergétique du taux de restitution de l'énergie dynamique lorsque $T \to \infty$. Pour cela nous faisons un choix particulier de $\theta$. C'est celui suggéré par C. Morawetz [MOR 77] et qui permet d'obtenir une estimation de l'énergie totale du système. Posons : $\theta = x - x_0$ où $x_0$ est un point arbitraire du plan contenant l'ouvert $\Omega$. On a alors $\mathrm{div}(\theta) = 2$ et $D\theta = I_d$ (matrice identité de $\mathbb{R}^2$). L'expression de $G^\theta$ devient :

$$G^{\theta} = \frac{1}{2} \int_{0}^{T} \int_{\Omega} |\dot{u}|^{2} + \frac{1}{T} \left[ \int_{\Omega} \dot{u} \nabla u \cdot \theta \right]_{0}^{T}, \quad (14)$$

et en notant qu’en multipliant l’équation des ondes (1) par $u$ et en intégrant sur $]0, T[ \times \Omega$ on obtient l’invariant de Maupertuis [GER 86] pour obtenir le même résultat que J. L. Lions dans [JLL 88] :

$$\left[ \frac{1}{2} \int_{\Omega} \dot{u} u \right]_{0}^{T} - \frac{1}{2} \int_{0}^{T} \int_{\Omega} |\dot{u}|^{2} + \frac{c^{2}}{2} \int_{\Omega} |\nabla u|^{2} = 0, \quad (15)$$
ce qui joint à la relation (14) nous donne, puisque l'énergie E est constante par rapport au temps (invariant classique obtenu en multipliant l'équation des ondes par $\dot{u}$) :

$$E + \frac{1}{T} \left[ \int_{\Omega} \dot{u} (\nabla u \cdot \theta + \frac{u}{2}) \right]_{0}^{T} = \Lambda^{\theta} = G^{\theta}. \quad (16)$$

Puisque $\dot{u}$ et $\nabla u$ sont bornés dans l'espace $C^0([0, T]; L^2(\Omega))$ on en déduit que pour ce choix de $\theta (= x - x_0)$ :

$$\lim_{T \to \infty} \Lambda^{\theta} = \lim_{T \to \infty} G^{\theta} = E. \quad (17)$$

Mais aussi pour tout champ $\theta$ égal à $\theta \cdot \nu$ sur les frontières $\Gamma_0$ et $\Gamma_i$ (puisque $\Lambda^\theta$ ne dépend que de cette quantité) :

$$E = \lim_{T \to \infty} \Lambda^{\theta} = \lim_{T \to \infty} \frac{1}{2T} \int_{0}^{T} \int_{\Omega} [|\dot{u}|^2 + c^2 |\nabla u|^2] \mathrm{div}(\theta) - 2(D\theta \nabla u \cdot \nabla u) \quad (18)$$

Nous noterons de la façon suivante les deux quantités ci-dessous qui sont égales :

$$\left\{
\begin{array}{l}
G^{\theta\infty} = \lim_{T \to \infty} \frac{1}{2T} \int_0^T \int_{\Omega} [|\dot{u}|^2 + c^2 |\nabla u|^2] \mathrm{div}(\theta) - 2(D\theta \nabla u \cdot \nabla u), \\
\\
\Lambda^{\theta\infty} = \lim_{T \to \infty} \Lambda^{\theta}.
\end{array}
\right.
\quad (19)$$

3. La troisième information que nous retirons de l'analyse précédente est que $G^{\theta\infty}$ (tout comme $\Lambda^{\theta\infty}$) sont des formes bilinéaires équivalentes à l'énergie et que les modes propres sont deux à deux orthogonaux vis-à-vis de chacune d'elles. Ce dernier point joint à la propriété de décroissance des modes locaux nous permet d'établir le résultat suivant que nous énonçons sous forme d'un théorème par commodité, mais qui en fait est surtout une recette de calcul. Une justification partielle est dans [DES 07]. Une autre plus complète est dans [WIL 07].

**Théorème 4** Désignons par $\theta_{\Gamma_i}$ un champ de vecteurs sur l'ouvert $\Omega$ égal $(x - x_0) \cdot \nu$ sur la frontière $\Gamma_i$ entre les deux milieux qui est supposée strictement intérieure à l'ouvert $\Omega$. Le milieux mou est celui qui est à l'intérieur de $\Omega$. Par ailleurs le support de $\theta_{\Gamma_i}$ est un voisinage de $\Gamma_i$. L'énergie des ondes stationnaires locales peut alors être approchée par $G^{\theta_{\Gamma_i}}$ et ce d'autant plus que les fréquences sont élevées et que la frontière $\Gamma_0$ de l'ouvert $\Omega$ est éloignée de $\Gamma_i$.

# 5. Quelques remarques sur l'approximation numérique

D'après les résultats de la section précédente, nous disposons de deux expressions identiques en continue pour évaluer les taux de restitution des énergies locale et intérieure. Mais sur le problème discret, ces deux quantités n'ont pas les mêmes propriétés et ne sont plus égales car la formule de Stokes n'est plus valable. La première est celle utilisant l'intégrale de bord et donc nécessite l'évaluation de la dérivée normale ce qui se fait à l'aide d'une formulation faible sur le schéma numérique. La seconde utilise une intégrale distribuée et ne nécessite pas l'évaluation de la dérivée normale. En outre, en utilisant un champ $\theta$ quasi constant au voisinage des zones où la solution est plus singulière (ce qui
est le cas des ondes locales car le gradient tend vers l'infini au voisinage de l'interface $\Gamma_i$ lorsque l'on monte en fréquences), on peut éliminer la perte de précision due à ces singularités. Les deux formules sont utilisables, mais la seconde est plus stable pour les hautes fréquences. Pour justifier ces propos, on peut donner l'explication suivante. Pour les comparer simplement prenons un exemple, celui d'une corde homogène de longueur $L$. Mais il est possible de compléter cette étude dans le cas général [WIL 07]. Le modèle est le suivant :

$$ egin{cases} 
\ddot{u} - c^2 \frac{\partial^2 u}{\partial x^2} = 0, \text{ pour } 0 < x < L \text{ et } t > 0, \quad u(0,t) = u(L,t) = 0, \\ 
\text{et la condition initiale : } u(x,0) = \sqrt{\frac{2}{L}} \sin\left(\frac{n\pi x}{L}\right), \quad \dot{u}(x,0) = 0, \quad \forall x \in ]0,L[. 
\end{cases} 
\quad (20) $$

Imaginons un maillage uniforme de taille $h$ du segment $]0, L[$ et notons $U_i, i = 1, N$ les valeurs aux nœuds de discrétisation de la solution approchée $u^h$ (en espace) de $u$. Par contre, nous conservons l'opérateur d'ordre deux en temps (semi-discrétisation). Notons alors que la solution exacte de l'équation des ondes est :

$$ u(x,t) = \sqrt{\frac{2}{L}} \sin\left(\frac{n\pi x}{L}\right) \cos\left(\frac{cn\pi t}{L}\right). \quad (21) $$

Les deux formules possibles pour le taux de restitution de l'énergie et après discrétisation en espace, sont les suivantes :

$$ egin{cases} 
\Lambda^{\theta h} = \frac{c^2 \theta(L)}{2T} \int_0^T |\frac{\partial u^h}{\partial x}(L,t)|^2 dt, \\ 
G^{\theta h} = \frac{1}{2T} \int_0^T \int_0^L [\dot{u}^h(x,t) + c^2 |\frac{\partial u^h}{\partial x}|^2(x,t)] dx dt. 
\end{cases} 
\quad (22) $$

Mais dans le cas d'une approximation numérique utilisant des polynômes du premier degré sur chaque élément et globalement continue vérifiant les conditions de Dirichlet homogènes aux deux extrémités, nous avons l'approximation classique de la dérivée normale en $x = L$ par (simple calcul analytique) :

$$ \left[ \frac{\partial u}{\partial x}(L,t) \right]^h = - \frac{U_N(t)}{h} + h \frac{\ddot{U}_N(t)}{6c^2} \quad (23) $$

Ceci nous conduit à la première estimation du taux de restitution de l'énergie :

$$ \Lambda^{\theta h} = \frac{c^2 L}{2T} \int_0^T \left[ -\frac{U_N(t)}{h} + h \frac{\ddot{U}_N(t)}{6c^2} \right]^2 dt. \quad (24) $$

Si nous appliquons cette formule approchée à la solution exacte (21) et que nous faisons tendre $T$ vers l'infini (pour évaluer l'ordre de la méthode), nous obtenons $(U_N)$ est la valeur en $x = L - h)$ :

$$ \Lambda^{\theta h} = \frac{c^2 L}{4} \left(\frac{U_N}{h}\right)^2 \left[1 + \frac{n^2 \pi^2 h^2}{6L^2}\right]^2. \quad (25) $$

Certains auteurs utilisent maladroitement une autre approximation de la dérivée normale qui n'est pas à conseiller, et qui consiste à poser :

$$ \left[ \frac{\partial u}{\partial x}(L,t) \right]^{hh} = - \frac{U_N(t)}{h}. \quad (26) $$
Cela conduit à une approximation peu stable et dont l'erreur ne peut pas être contrôlée en faisant tendre *h* vers zéro lorsque par exemple l'indice *n* augmente (les hautes fréquences ne sont pas capturées par un maillage plus fin). L'autre méthode est celle qui utilise l'expression $G^θ$. Toujours dans le cas d'une méthode d'éléments finis d'ordre un, nous obtenons cette fois :

$$G^{\theta h} = \frac{1}{2T} \int_{0}^{T} \int_{0}^{L} [|\dot{u}^h|^2(x, t) + c^2|\frac{\partial u^h}{\partial x}(x, t)|^2] \frac{\partial \theta}{\partial x}(x) dx dt. \quad (27)$$

Or, nous avons vu que pour le modèle continu, seule importait la valeur de $\theta$ en $x = L$. Choisissons par exemple

$$\theta(x) = L(x - L + h) \text{ pour } L - h < x < L \text{ et } \theta(x) = 0 \text{ sinon.} \quad (28)$$

Un simple calcul nous donne cette fois :

$$G^{\theta h} = \frac{c^2 L}{4} \left(\frac{U_N}{h}\right)^2 \left[1 + \frac{n^2 \pi^2 h^2}{3L^2}\right]. \quad (29)$$

Pour comparer les deux approches, calculons la valeur exacte de l'énergie pour la solution continue. A partir de (21), nous obtenons en faisant intervenir la valeur de $u$ en $x = L-h$ :

$$E = \frac{c^2 n^2 \pi^2}{2L^2} = \frac{c^2 L}{4} \left(\frac{U_N}{h}\right)^2 \left[\frac{\frac{nh\pi}{L}}{\sin\left(\frac{nh\pi}{L}\right)}\right]^2 \approx \frac{c^2 L}{4} \left(\frac{U_N}{h}\right)^2 \left[1 + \frac{n^2 \pi^2 h^2}{3L^2} + \frac{n^4 \pi^4 h^4}{15L^4}\right] \quad (30)$$

Ceci nous conduit aux conclusion suivantes. En basses fréquences les deux méthodes d'évaluation de l'énergie sont du même ordre (à condition de ne pas utiliser l'approximation inconsistante de la dérivée normale). Par contre, ce résultat suppose que la variation en temps est connue exactement (semi-discrétisation en espace). Lorsqu'un schéma numériserique est utilisé pour l'intégration en temps la méthode utilisant $\Lambda^\theta$ fait intervenir la dérivée seconde en temps qui est n'est pas estimée de façon fiable dans la norme de l'espace $C^0([0, T]; L^2(\Omega))$, alors que celle basée sur $G^\theta$ n'utilise que la dérivée première dont l'approximation est stable et consistante avec tous les schémas d'intégration classiques (voir P.A. Raviart- J.M. Thomas [RAV 83]). Dans les calculs qui suivent, nous avons utilisé un schéma aux différences centrées. Cela nous permet d'affirmer que la méthode $G^\theta$ est la plus recommendable.

# 6. Quelques exemples numériques

Commençons par un grand carré à l'intérieur duquel il y a un autre carré plus petit et dans un matériau plus mou. Les géométries sont visualisées sur la figure 5 ainsi que le maillage utilisé pour chacun des deux cas (zone molle centrée et décentrée).

Cette géométrie est plus pénalisante que le cas d'un disque lorsque la zone molle est centrée car on perd la symétrie de révolution et les angles sont à l'origine de réfraction complexes. Nous avons utilisé différents types de conditions initiales suivant que ces dernières sont localisées dans la partie molle ou celle qui est la plus dure. Dans chaque cas nous avons tracé en fonction du temps l'évolution de $G^{\theta I}$ (énergie des ondes intérieures lorsque $T \to \infty$) et $G^{\theta L}$ (énergie des ondes locales) ainsi que la somme des deux qui est
<figure><img src="image_5.png" /><figcaption>Figure 5. Maillages utilisés dans les calculs et allures des conditions initiales</figcaption></figure>

$G^θ$. Les résultats sont reportés sur la figure 6 (excitation intérieure) et 7 (excitation extérieure). On constate que la stabilisation (obtention de la valeur asymptotique en temps) est très rapide et correspond sensiblement au temps nécessaire à l'onde pour faire un aller retour dans la direction la plus grande de l'ouvert en tenant compte d'un rebond sur l'interface.

Nous avons reporté sur les figures 8 les résultats obtenus pour une condition initiale située dans la zone dure. On voit que les conclusions sont similaires et que la zone molle fonctionne bien comme un piège à onde. Cet aspect souligne évidemment la fragilisation de la pièce mécanique soumise à des sollicitations dynamiques. C'est le coté négatif. Mais dans certaines applications, ce piège à ondes peut aussi catalyser des procédés de contrôle anti-bruit en couplant les zones molles avec des actionneurs permettant de supprimer les mouvements. Dans ce cas ils devraient d'ailleurs être localisés dans la zone molle. Mais c'est une autre étude que nous discuterons dans [WIL 07].

Lorsque la zone molle est décentrée les symétries du problème sont perdues et les résultats sont sur les figures 8 avec des conditions initiales dont le support est dans la zone molle (en haut) et des conditions initiales dans la zone dure (en bas). Dans les deux cas la condition initiale est ici en haute fréquence (longueur d'onde de l'ordre du dixième de la dimension de la zone molle). L'évaluation de l'énergie des ondes locales est donc très convenablement approchée par la méthode du taux de restitution de l'énergie dynamique. Mais on peut aller plus loin dans l'interprétation des résultats, compte tenu du caractère local de ces ondes. En effet l'énergie de ces ondes se formule sous forme d'une intégrale distribuée le long de l'interface entre les deux milieux. L'intégrande est notée *g* et l'on a :

$$ \Lambda^{\theta} = \frac{1}{2T} \int_{0}^{T} \int_{\Gamma_i} g(s, t)(\theta, \nu)dsdt. \qquad (31) $$

Il est donc naturel d'associer *g* à une force thermodynamique liée à la décohésion de l'interface. Le calcul peut être fait en utilisant un champ $\theta$ local et en utilisant la méthode
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Evolution des taux de restitution dynamiques de l'énergie pour une condition initiale dans la zone molle (en haut basses fréquences et en bas hautes fréquences)</figcaption></figure>

$G^θ$ de façon à assurer une meilleure stabilité, principalement pour les hautes fréquences. Nous avons reporté cette grandeur sur les figures 9 pour une excitation hautes fréquences dans le cas de conditions initiales localisées dans la zone molle (graphique en haut) ou dans la zone dure (graphique en bas).

## 7. Conclusion

Dans cet article, avons mis en évidence la localisation d'ondes stationnaires (modes propres) à l'interface entre deux matériaux de caractéristiques distinctes. Le point im-
<figure><img src="image_8.png" /></figure>

Ec+Ep (trait epais) – Githetas (traits fins avec symboles)

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Evolution des taux de restitution dynamiques de l'énergie pour une condition initiale dans la zone dure (en haut basses fréquences et en bas hautes fréquences)</figcaption></figure>

portant est la possibilité d'extraire l'énergie de ces ondes locales qui ici sont du type de celles découvertes par A.E. Love et H. Lamb [LAM 19] [FUN 65] pour des ondes progressives, à l'aide d'un taux de restitution de l'énergie ne nécessitant pas le calcul des modes propres locaux. Grâce au fait que les ondes locales soient à décroissance rapide (exponentielle) depuis l'interface entre les deux matériaux et dans la direction du milieu le plus dur. Cela permet de découpler les contributions des ondes locales et intérieures dans le taux de restitution dynamique. On peut aussi en extraire une grandeur locale le long de l'interface entre les deux milieux qui s'apparente à une force thermodynamique et qui par conséquent est un candidat pour une étude de l'endommagement de cette interface. C'est pourquoi nous avons utilisé la terminologie de taux de restitution dynamique
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Evolution des taux de restitution dynamiques de l'énergie pour une condition initiale hautes fréquences en haut dans la zone molle et en bas dans la zone dure (cas d'un carré décalé)</figcaption></figure>

de l'énergie. Cependant il faut retenir que la grandeur qui est mécaniquement significative est la moyenne quadratique vis-à-vis du temps d'une grandeur définie sur la ligne d'interface. La motivation physique de notre étude était de trouver une grandeur énergétique, donc fiable dans des calculs numériques, permettant d'envisager une analyse de l'endommagement des interfaces. Mais il y a aussi les aspects liés à la réduction des vibrations qui apparaîssent comme des bénéficiaires de cette localisation d'énergie vibratoire en créant des pièges à onde dans des zones molles où seraient disposés des systèmes de contrôle dont nous proposons l'analyse dans un travail à venir [WIL 07].
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Représentation des taux de restitution dynamiques de l'énergie locaux à différents instants pour une condition initiale hautes fréquences dans la zone molle</figcaption></figure>

## 8. Bibliographie

[DES 07] PH. DESTUYNDER, « Vibration des systèmes couplés », *Hermès, Londres*, 2007.

[DES 06] PH. DESTUYNDER, « Analyse et traitement des images numériques », *Hermès, Londres*, 2006.

[DD 81] PH. DESTUYNDER, M. DLAOUA, « Sur une interprétation mathématique de l'intégrale de Rice en mécanique de la rupture fragile », *Math. Meth. in the Applied Scie.*, vol. 3, 1981.

[WIL 07] PH. DESTUYNDER, O. WILK, « Analysis of local waves at the interface between two media », *à paraître, Rapport CNAM*, 2006.

[FRI 61] K.O. FRIEDRICHS, « Perturbation of spectra in Hilbert spaces », *Courant institute, New York*, 1961.

[FUN 65] Y. C. FUNG, « Foundations of solid méchanics », *Prentice Hall, New-York*, 1965.

[GER 86] P. GERMAIN, « Mécanique », *Ellipses, Paris*, 1986.

[KAT 66] T. KATO, « Perturbation theory for linear operators », *Springer-Verlag, Berlin*, 1966.

[JLL 88] J.L. LIONS, « Contrôlabilité exacte perturbations et stabilisation de systèmes distribués. Tome 1. », *Collection RMA, Masson, Paris*, 1988.

[LAM 19] H. LAMB, « On the propagation of tremors over the surface of an elastic solid », *Phil. Trans. Roy. soc. Londres, ser. A*, vol. 203, p. 1-42, 1916.

[MOR 77] C. MORAWETZ, J. RALSTON, JW. STRAUSS, « Decay of the solution of wave equation outside a non trapping obstacle », vol. 0-30, p. 447-508, 1977.

[RAV 83] P.A. RAVIART, J.M. THOMAS, « Introduction à l'analyse numérique des équations aux dérivées partielles », *Masson, Paris*, 1983.

[RIE 72] F. RIESZ, JB.S.Z NAGY, « Leçons d'analyse fonctionnelle », *Gauthier-Villars, Paris*, 2006.