Contributions à la simulation numérique en élastodynamique

C’est au cœur des années 1970 qu’est véritablement né le développement de ces méthodes pour répondre aux besoins de l’industrie pétrolière en matière de simulation d’expériences de détection de réservoirs de pétrole par ondes sismiques [Boore 72], [Alford 74], [Kelly 76]. Des méthodes de différences finies ont d’abord été développées mais, très vite, les méthodes variationnelles pour la discrétisation en espace se sont imposées, en particulier parce qu’elles permettaient de traiter des milieux hétérogènes et des conditions aux limites (de surface libre notamment) sans générer de problèmes d’instabilité numérique [Ciarlet 78], [Bamberger 80], [Virieux 86]. C’est ainsi que, progressivement, des codes de calcul ont été développés à partir de nombreuses méthodes comme les éléments finis lagrangiens en déplacement, les éléments finis mixtes en contrainte-déplacement et plus récemment les méthodes de Galerkin discontinues (liste non exhaustive), chaque méthode ayant ses propres avantages et inconvénients, de sorte qu’aucune ne s’est encore imposée aux autres [Komatitsch 97], [Fauqueux 03].

En revanche, il est aujourd’hui communément admis qu’il est utile, tant pour des raisons de coût de calcul que de précision, d’utiliser les versions d’ordre élevé de ces méthodes, que l’on pourrait regrouper sous le nom de méthodes éléments finis (pour évoquer le choix d’inconnues polynomiales par morceaux sur un maillage spatial) spectraux (terminologie plus ou moins consacrée pour parler d’ordre élevé, au moins sur les maillages de type héxaédriques ou quadrangulaires [Cohen 01]) possiblement discontinus (pour suggérer la diversité des stratégies de discrétisation en espace). Il est également admis que, pour pouvoir traiter des problèmes de très grande taille sans faire exploser le temps de calcul, il faut se tourner vers des stratégies de discrétisation explicites et, le cas échéant, vers des techniques de condensation de masse appropriées [Fauqueux 03], [Joly 08], qui sont alors soumises à des conditions de stabilité de type CFL.

Les champs d’application de ces méthodes n’ont cessé de se diversifier avec le temps, pour s’étendre notamment à la géophysique terrestre ou au contrôle non destructif. Les applications potentielles au contrôle non destructif par ultrasons ont motivé nos recherches. Elles sont à l’origine du développement de cette thèse qui s’intéresse à des milieux constitués d’assemblages de domaines homogènes isotropes par morceaux, séparés par des interfaces minces, permettant d’aborder un large spectre d’applications réelles. Ainsi, le lecteur pourra se figurer un assemblage de pièces métalliques, constitué de matériaux éventuellement différents et séparés par une fine couche de collage ou de soudure ayant des propriétés mécaniques distinctes des deux matériaux qu’elle sépare.

De cette façon, on peut représenter une cuve de réacteur nucléaire ou une aile d’avion formée de matériaux composites dans lesquelles on souhaiterait, par exemple, détecter la présence d’une fissure. Un tel cas de figure peut également trouver des applications dans le domaine médical : on considère alors que le milieu de propagation est une partie du corps humain que l’on peut souvent assimiler à un assemblage de zones homogènes séparées par de fines membranes. Dans ce cas, on cherchera à détecter par ultrasons la présence d’une tumeur, d’un calcul, etc. [Audière 11].

Les méthodes numériques évoquées plus haut et, par conséquent, les nombreux codes de calcul qui y sont associés, sont bien entendu susceptibles de résoudre ce type de problèmes. Elles seront toutefois fortement pénalisées dans deux situations particulières :

♦ Premier cas : on considère ici des interfaces infiniment minces et on se place dans le cadre de milieux homogènes où les ondes S, ou ondes de cisaillement, se propagent à une vitesse VS très inférieure à la vitesse VP des ondes P, ou ondes de pression. C’est ce qui se passe dans les milieux mous à l’instar des tissus vivants, comme le montre l’exemple du foie traité sous l’angle physique dans [Audière 11]. Dans ce cas – nous l’expliquerons de façon détaillée plus tard – les méthodes classiques sont pénalisées par la petite valeur du rapport VS/VP , ne serait-ce que par son influence sur la condition de stabilité. En outre, les deux types d’ondes ont des dynamiques et des longueurs d’onde très différentes. Nous aimerions donc pouvoir les traiter séparément, avec des maillages ou des discrétisations adaptés aux propriétés respectives de chaque type d’onde, ce que les méthodes standard, qui travaillent sur un seul maillage de calcul, ne permettent pas de faire a priori.

C’est pourquoi nous avons cherché à exploiter numériquement une technique employée dans d’autres contextes comme la mécanique des fluides et l’électromagnétisme : la décomposition d’un champ de vecteurs en potentiels (décomposition de Helmholtz) qui sont deux quantités scalaires quand on se place en dimension 2. Cette méthode permet de découpler les ondes P et les ondes S et, par conséquent, de les traiter a priori de façon adaptée en utilisant des approximations numériques différentes pour chacun des deux potentiels qui vérifient dans chaque sous-domaine une équation des ondes scalaire. C’est essentiellement au niveau des conditions aux limites sur les bords physiques et des conditions de transmission entre sous-domaines que se situe, comme nous le verrons, toute la difficulté du problème. Ce travail fait l’objet de la première partie de la thèse.

♦ Second cas : l’épaisseur η des interfaces minces est très petite devant les longueurs d’onde caractéristiques du problème, suffisamment toutefois pour influencer de façon non négligeable les phénomènes de réflexion-transmission des ondes. Dans ce cas, l’utilisation d’un maillage – ne serait-ce que localement – de calcul à l’échelle η va fortement pénaliser les méthodes classiques en raison de l’augmentation du nombre de degrés de liberté et de l’incidence sur la condition CFL. Cette problématique « multi-échelle » est commune à de nombreux problèmes de propagation d’ondes. Une parade purement numérique consiste à se tourner vers des techniques de raffinement de maillage non conforme couplées à des méthodes de pas de temps local, au sujet desquels nous renvoyons le lecteur à [Bécache 05], [Rodriguez 04]. Dans cette thèse, nous allons nous tourner vers une solution plus analytique afin d’éviter le recours à tout maillage à l’échelle η : il s’agit de concevoir des conditions de transmission dites effectives permettant de calculer une approximation de la solution exacte à un certain ordre par rapport au petit paramètre η. La construction de ces conditions repose sur une analyse asymptotique de cette solution par rapport à η. Ces travaux font l’objet de la seconde partie de la thèse.

Introduction à la décomposition en potentiels

L’idée générale de l’utilisation du découplage en potentiels du champ de déplacement découle de la difficulté du traitement de l’opérateur différentiel apparaissant dans les équations volumiques, tant en élastodynamique que dans d’autres domaines comme l’électromagnétisme ou l’élastostatique. Décomposer le champ recherché en potentiels permet de résoudre un problème avec moins d’inconnues, ou alors de se ramener à des équations aux dérivées partielles ayant une forme plus familière, comme des équations de Poisson ou de d’Alembert par exemple.

De façon générale, on parlera de décomposition en potentiels d’un champ de vecteur u dès que l’on sera capable de décomposer ce champ de vecteur comme la somme d’un gradient et d’un rotationnel.

En pratique, ceci ne se traduira que par une hypothèse sur les sources du problème. Une seconde question est l’unicité de cette décomposition. Nous verrons que, de façon générale, cette décomposition n’est pas unique et que, par exemple dans un domaine borné, le choix d’une décomposition est lié à un choix de conditions aux limites à imposer aux potentiels. Là encore, la question de ce choix ne se posera pas pour nous et le choix des conditions à imposer sur les potentiels nous sera imposé par la physique du problème que nous voulons traiter, en l’occurrence les conditions aux limites physiques de Dirichlet ou de Neumann.

La suite de ce chapitre s’organise comme suit : nous commençons par rappeler les résultats bien connus sur la (les) décomposition(s) de Helmholtz-Hodge d’un champ de vecteur et son (leurs) application(s) dans le contexte de l’électromagnétisme. Même si ce ne sont pas ces décompositions que nous utiliserons ensuite, ces résultats nous permettront d’évoquer les résultats de nonunicité de la décomposition de Helmholtz. Nous nous intéresserons dans un second temps à la recherche des solutions de l’équation de l’élastodynamique dans le plan homogène R² entier sous la forme d’une décomposition en potentiels  . Sous réserve d’hypothèses adéquates sur les données du problème, nous établirons en particulier l’équivalence entre le problème initial en déplacement  et sa formulation en potentiels, et montrerons en particulier que cette équivalence permet d’illustrer de façon élégante la décomposition de l’onde élastique comme la somme d’une onde de pression et une onde de cisaillement se propageant indépendamment l’une de l’autre.

Influence d’un bord ou d’une interface

Influence d’un bord

Dans le volume, les ondes de pression et de cisaillement se propagent de façon découplée. Cependant, cette propriété est perdue dès que l’onde touche un bord ou une interface. Pour s’en rendre compte, il suffit par exemple de regarder comment se réécrivent les conditions de type Dirichlet (bord encastré) ou Neumann (bord libre) dans la formulation en potentiels.

On a vu que les conditions de bord sur le champ de déplacement sont données par .

Conditions de transmission effectives en élastodynamique

Ce problème est loin d’être nouveau en mathématiques, notamment dans le contexte de la propagation des ondes. Dès le début des années 1990, plusieurs auteurs se sont posé la question de proposer des conditions aux limites équivalentes, à des fins numériques, pour résoudre des problèmes de diffraction d’ondes (acoustiques ou électromagnétiques) par des obstacles recouverts par une couche mince. Cela s’est d’abord fait en régime harmonique dans le cas de couches minces homogènes [Engquist 93], [Bendali 96], puis dans le cas de couches minces périodiques lorsque la période est proportionnelle à η : on parle alors de surfaces rugueuses [Ammari 98], [Delourme 12]. Les problèmes transitoires, qui posent des difficultés supplémentaires liées à des questions de stabilité, ont notamment été abordés dans [Haddar 02], y compris dans le cas de matériaux non linéaires. Dans le cas de couches homogènes, si l’on excepte [Engquist 93], la méthode d’analyse qui semble s’être imposée comme référence repose sur une technique de scaling dans la couche mince puis sur un ansatz de représentation de la solution sous la forme d’un développement en séries de puissances de η. Le cas de couches périodiques est un peu plus délicat et nécessite de combiner des techniques de type couches minces avec des techniques d’homogénéisation, ce qui se fait à l’aide de développements asymptotiques plus sophistiqués comme les développements multi-échelles [Dauge 10] ou les développements asymptotiques raccordés [Maz’ya 91]. Dans tous les cas, on aboutit à des conditions d’impédance généralisées (Generalized Impedance Boundary Conditions, GIBC en anglais) largement répandues dans la communauté des ingénieurs [Senior 91]. De telles conditions apparaissent aussi dans le contexte de la diffraction d’ondes par des obstacles fortement conducteurs [Haddar 07], [Antoine 05], où l’effet couche mince est alors remplacé par l’effet de peau.

Si le cas des conditions de transmission a été moins abordé, nous pouvons toutefois citer les travaux de [Chun 10] en électromagnétisme transitoire dans le cas d’interfaces homogènes – travaux qui ont largement inspiré les nôtres – et ceux de [Delourme 12] pour la transmission d’ondes acoustiques et électromagnétiques en régime harmonique dans le cas d’interfaces hétérogènes rugueuses. Il s’avère que le traitement des problèmes de transmission est plus délicat d’un point de vue purement technique, pour l’analyse d’erreur notamment. Il mène rapidement à des conditions plus complexes et peut soulever des questions de stabilité qui sont nouvelles. Enfin, la question des conditions de transmission en élasticité linéaire a déjà été abordée dans le cas statique [Geymonat 14] mais pas dans le cas dynamique (fréquentiel) et encore moins en instationnaire.

Table des matières

Introduction
1 Rappels sur les équations de l’élastodynamique et motivations
1.1 Le modèle élastodynamique classique
1.1.1 Notations préliminaires
1.1.2 Équations du problème
1.1.3 Formulation variationnelle
1.1.4 Conservation d’énergie
1.2 Propagation des ondes élastiques en milieu isotrope
1.2.1 Loi de Hooke et coefficients de Lamé
1.2.2 Introduction au découplage des ondes élastodynamiques
1.2.3 Ondes planes élastiques et relation de dispersion
1.2.4 Comparaison des vitesses des ondes P et des ondes S
1.3 Une méthode numérique classique et ses inconvénients
1.3.1 Semi-discrétisation en espace
1.3.2 Discrétisation en temps
1.3.3 Conservation d’énergie discrète et limite du pas de temps
1.3.4 Inconvénients des méthodes classiques en déplacement
2 Principe de la décomposition en potentiels
2.1 Introduction à la décomposition en potentiels
2.1.1 La décomposition de Helmholtz-Hodge
2.1.2 À propos de la non-unicité de la décomposition
2.2 Application aux équations de l’élastodynamique homogène dans l’espace libre
2.3 Influence d’un bord ou d’une interface
2.3.1 Influence d’un bord
2.3.2 Influence d’une interface
2.4 La décomposition en potentiels en physique
2.4.1 Les potentiels en électromagnétisme
2.4.2 Les potentiels en élasticité statique
2.4.3 Les potentiels de Lamé
2.4.4 Le potentiel de Somigliana
3 Formulation en potentiels du problème de Dirichlet
3.1 Conditions de bord adaptées
3.1.1 Équivalence des formulations
3.1.2 Une première formulation variationnelle du problème couplé
3.2 Énergie associée à la formulation en potentiels avec bords encastrés
3.2.1 Une intuition qui se révèle fausse
3.2.2 L’énergie effectivement conservée
3.2.3 Comparaison avec l’énergie élastodynamique
3.2.4 Retour sur le bon cadre fonctionnel
3.3 Discrétisation du problème de Dirichlet
3.3.1 Semi-discrétisation en espace
3.3.2 Discrétisation totale en espace-temps par différences finies
3.4 Résultats numériques
3.4.1 Présentation des simulations
3.4.2 Comparaison avec la solution de la méthode en déplacement
3.4.3 À propos de la source
3.4.4 Maillages et paramètres
3.4.5 Simulations dans un carré
3.4.6 Simulations dans un anneau
3.4.7 Simulations dans une géométrie avec un coin rentrant
3.5 Étude du problème fréquentiel
3.5.1 Formulation du problème
3.5.2 Les difficultés théoriques
3.6 Conclusion
4 Le problème de Neumann et ses défis
4.1 Formulation du problème de Neumann avec inconnue auxiliaire en vitesse
4.1.1 Construction de la formulation
4.1.2 Équivalence des formulations
4.1.3 Formulation variationnelle
4.1.4 Conservation d’énergie
4.1.5 Semi-discrétisation en espace
4.1.6 Discrétisation en temps
4.1.7 Un phénomène d’instabilité numérique
4.1.8 Étude du problème fréquentiel
4.2 Une formulation avec élimination de l’inconnue auxiliaire
4.2.1 Construction de la formulation
4.2.2 Formulation variationnelle et identité d’énergie
4.3 Autres formulations
4.3.1 Formulation avec deux inconnues auxiliaires
4.3.2 Comparaison des valeurs propres pour les différentes formulations
4.3.3 Un mot sur la décomposition du domaine en une partie déplacement et une partie potentiels
4.4 Conclusion
Conclusion

Télécharger le rapport complet