Pascal Poullet
résumé
L’essentiel de mes activités de recherche peut s’articuler autour de 3 thèmes distincts :
– Modélisation numérique d’écoulements de fluides incompressibles et applications
– Méthodes des inconnues incrémentales et applications
– Résolution de lois de conservation en formulation non conservative
Dans le premier thème, des travaux sur l’étude d’une des méthodes les plus utilisées pour simuler des écoulements de fluides incompressibles ont été réalisés. Il s’agit d’un domaine de recherche très concurrentiel, dans la mesure où de nombreuses applications industrielles, environnementales et militaires requièrent le calcul de la solution de problèmes de ce type. Dans des configurations très générales, les méthodes de projection sont connues pour permettre de calculer de manière très efficace les solutions des équations instationnaires de Navier–Stokes (NSE), sous forme vitesse–pression, pour un fluide incompressible visqueux. La difficulté de traitement de ces équations est importante du fait d’être des équations aux dérivées partielles non linéaires dépendant du temps et de l’espace et surtout couplant les inconnues vitesse et pression au moyen de la contrainte d’incompressibilité. Les méthodes de projection sont basées sur la décomposition de Hodge de l’espace des fonctions de carré intégrable, exploitant le fait que cette décomposition peut s’écrire comme somme directe orthogonale de l’espace de Sobolev H 1 0 et de l’espace fonctionnel constitué de gradient de fonctions de carré intégrable (projecteur de Leray). Ces méthodes utilisent donc une technique à pas fractionnaires, consistant à calculer lors d’une première étape, une vitesse approchée à un temps intermédiaire, puis, lors d’une seconde étape, de la corriger au moyen d’une projection sur ’espace à divergence nulle. Lors de cette seconde étape, la pression est souvent directement calculée et la nouvelle vitesse est actualisée ensuite (méthode de projection incrémentale).
Visant à améliorer ce type de méthodes, une première étude a consisté à faire une étude comparative d’une des principales variantes de ces méthodes (méthodes de correction de la pression). Parmi cette variante, la méthode de projection-pénalisation a été introduite sous une forme différente de celle qui était connue jusqu’alors. Une étude théorique a permis de calculer l’ordre d’approximation du schéma en temps proposé, et cette étude a été complétée par des simulations numériques. En fait, cette méthode de projection-pénalisation n’avait jamais été implémentée auparavant, du fait que le paramètre de pénalisation, r 1 , n’était originellement censé valoir que l’inverse du carré du pas de temps .t (si .t est souvent petit, r 1 est donc une valeur relativement grande). Ce terme de pénalisation étant un opérateur mal conditionné, il nous a fallut utiliser les techniques les plus sophistiquées de l’algèbre linéaire (factorisation incomplète de type IKJ) pour rendre l’implémentation de la méthode efficace. Des travaux supplémentaires ont été poursuivis et nous ont menés à proposer un nouvel algorithme visant, au cours des simulations, à omettre l’étape de projection qui s’avérait être la plus coûteuse de l’algorithme d’origine. De ce fait, on obtient une nouvelle méthode un peu moins précise que la méthode d’origine mais beaucoup moins coûteuse en temps de calcul. Les estimations théoriques de l’erreur de fractionnement ont attestés de la perte de précision et les simulations pour des problèmes tests prouvent que les résultats obtenus demeurent acceptables.
Le deuxième thème porte sur des méthodes multiniveaux pour l’approximation d’équations aux dérivées partielles (EDP) provenant de la Physique. Plus particulièrement, la méthode des Inconnues Incrémentales (II), dont l’efficacité est connue pour les EDP elliptiques comme le problème de Poisson avec conditions aux limites de Dirichlet, a été souvent comparée aux méthodes multigrilles classiques. La méthode des II a été introduite au début des années 90 dans le but de permettre une hiérarchisation des inconnues des équations de la Mécanique des
fluides, après discrétisation par différences finies. C’est donc au sein de l’étude des systèmes dynamiques provenant de la Mécanique des fluides, qu’a émergé l’idée d’un traitement différent des structures de taille plus ou moins grande d’un écoulement. Tout un concept a été développé à partir de cette idée, et a donné naissance aux méthodes de type Galerkin non linéaires. Une littérature abonde dans le sujet et le challenge d’une méthode efficace de simulation directe d’écoulements turbulents est toujours très actuel.
Un de mes travaux, réalisé sur les méthodes utilisant les II, a consisté à utiliser les II classiques pour la résolution d’un problème non linéaire stationnaire dans un domaine borné (introduit par R. Temam sous forme évolutive). Ce problème, présentant un terme visqueux et une non linéarité quasiment identiques à ceux des équations NSE, l’approche développée a consisté à présenter des stratégies de méthodes de type Newton inexact et Newton tronquées utilisant les II. Plus précisément, le problème considéré ayant une jacobienne non symétrique,
la méthode de Newton a été implémentée par le biais d’algorithmes adéquats de projection sur un sous-espace de Krylov. Faisant varier le nombre de grilles d’II (ou faisant varier le nombre de points sur la grille la plus grossière), l’efficacité du préconditionneur est obtenue et croît avec l’ellipticité de l’équation. Ces résultats sont confirmés par la théorie et pour une viscosité en-dessous d’un certain seuil, on obtient une stratégie optimale avec un nombre de niveaux de grille d’II pas le plus élevé (correspondant à la grille la plus grossière non réduite à un point).
Une autre étude a consisté à améliorer, en utilisant les II, les algorithmes classiques de résolution des problèmes de Stokes et Stokes généralisé. Ces problèmes sont obtenus après simplifications des équations de NSE, le premier étant le modèle pour un fluide incompressible à très faible nombre de Reynolds, tandis que le second s’obtient après discrétisation semi-implicite des équations évolutives. La difficulté principale du problème provenant du couplage des inconnues vitesse-pression (V. thème 1), une des approches classiques est d’utiliser un algorithme d’optimisation sous contraintes comme l’algorithme d’Uzawa. Appliquer cet algorithme dans notre cas, revient à utiliser une méthode de Krylov pour résoudre un système linéaire qui a pour matrice le complément de Schur de la matrice de rigidité. L’apport significatif de mon travail a été de proposer une nouvelle définition des II (appelées SIU) basée sur une décomposition triadique des inconnues discrètisées par di.érences finies. L’intérêt de cette nouvelle définition est qu’elle octroie naturellement (sans aucune interpolation) la possibilité de hiérarchiser un maillage décalé de type MAC, tout en conservant la propriété inhérente aux cellules de type MAC sur les différents niveaux de grille. Cette propriété, garantissant la stabilité de la discrétisation (liée à la condition LBB en éléments finis mixtes), est très importante si l’on souhaite calculer les solutions à grand nombre de Reynolds. Par des résultats numériques, nous prouvons l’efficacité de deux nouvelles approches basées sur l’algorithme d’Uzawa : l’un utilisant les SIU en vitesse comme préconditionneur du problème en vitesse pour le problème de Stokes, l’autre utilisant les SIU en pression pour le problème en pression (le complément de Schur) pour le problème de Stokes généralisé.
Mon intérêt s’est porté aussi sur le calcul des solutions stationnaires de l’équation de propagation des ondes et principalement quand cette propagation était perturbé par corps étranger de forme non connexe, un diffuseur. C’est un problème modèle provenant de l’ingénierie posé avec un corps solide métalique en forme de fer à cheval, pour lequel le challenge réside dans l’efficacité de la méthode pour de basses fréquences. Ce problème, nécessitant la résolution de l’équation d’Helmholtz dans un domaine non borné, peut être traité en utilisant les espaces de Sobolev à poids pour lesquels l’inégalité de Poincaré n’est pas valable. Mais nous avons choisi l’approche souvent considéré numériquement, de compléter le problème par l’ajout de conditions aux limites absorbantes (ABC pour l’acronyme anglais), destiné à négliger les réfractions non physiquement admissibles de la solution, au bord d’un domaine fictif englobant le diffuseur. Lors d’une première étude, les II classiques ont été adaptées à la résolution du problème elliptique indéfini complété de conditions aux limites de type Robin (ABC au second ordre), posé dans le plan complexe. Un algorithme a été proposé, basé sur une méthode de gradient pour l’équation normale préconditionnée par les II, du fait du caractère non hermitien de l’opérateur obtenu. Les résultats obtenus attestent de la robustesse et de l’efficacité de l’algorithme même pour de relativement faibles fréquences. Toutefois, on retrouve le comportement similaire des II pour le problème non linéaire de type Navier-Stokes, à savoir que la grille la plus grossière d’II doit comporter un nombre suffisant de points de discrétisation pour obtenir la meilleure efficacité de l’algorithme (cinq noeuds par longueur d’onde). C’est ce verrou que nous avons cherché à analyser et à surmonter par la suite. Lors d’une étude théorique, nous avons appliqué une stratégie introduite par H. Yserentant pour résoudre le problème modèle simplifié 1D suivant :
–uxx–k2u = f, in ]0,1[ (1)
u(0) = 0, (2)
ux(1) = iku(1) (3)
Il s’agit bien du modèle simplifié provenant de l’équation des ondes en 1D, quand on considère que sa solution se propage de la gauche vers la droite et qu’après changement d’échelle sa solution vérifie la condition de Sommerfeld à l’infini (. représentant le complexe tel que . 2 = - 1).
Après avoir introduit, l’espace fonctionnel V des fonctions de carré intégrable, ayant leur dérivée première de carré intégrable et vérifiant la condition (1), à partir de la forme faible du problème (1-2-3), on prouve aisément que la forme bilinéaire sous-jacente est continue sur V ×V . Bien qu’elle ne soit pas coercive, on prouve qu’elle vérifie une inégalité de Gärding et que pour f fonction de carré intégrable et k > 0, le problème possède une unique solution dans V . Au moyen du formalisme inhérent aux bases hiérarchiques d’éléments finis et à un argument introduit par Schatz (utilisé par Yserentant), nous prouvons qu’il est possible d’obtenir un gain significatif de conditionnement du problème, dès que grille la plus grossière est choisie suffisamment représentative. Bien que les résultats théoriques sou.rent d’optimalité (on trouve que la grille grossière doit être de l’ordre de k 4 inconnues), on comprend bien que c’est le manque de V -ellipticité qui explique la limitation sur la grille grossière.
Lors d’une nouvelle étude, une méthode numérique est proposée pour surmonter cette limitation causée par les phénomènes survenant à basse fréquence. Une redéfinition des II a été proposée à partir de la solution fondamentale en 1D de l’équation de Helmholtz discrétisée, permettant d’avoir une interpolation exacte pour cette équation. En 1D, on trouve des poids d’interpolation égaux à sin(kh)/ sin(2kh), avec h le pas discrétisation au lieu de 1/2 dans la définition classique. Et à l’aide de cette nouvelle définition, l’amélioration du conditionnement de la matrice de rigidité dans la base des II est obtenue même avec la grille grossière la plus petite. L’extension aux problèmes de dimensions supérieures ne peut pas se faire directement dans la mesure où la solution fondamentale (en l’absence de source) de l’équation de Helmholtz s’écrit comme une intégrale suivant toutes les directions du problème. A la manière des méthodes intégrales, la stratégie proposée est de calculer les coefficients d’interpolation en certains points de la cellule, à partir des valeurs aux bords de cette cellule (que constituent la classe de référence). En dimension 2, par exemple, on choisit donc la construction d’une cellule à un niveau l, que l’on divise en 4 au niveau précédent, et les points intérieurs au niveau l sont obtenus comme l’union des points du bord des 4 cellules du niveau l - 1, privée du bord de la cellule du niveau l. Par conséquent, l’incrément à un certain niveau est associé aux points de la croix intérieure de la cellule courante en fonction des valeurs aux points du bord de celle-ci. Du fait de la densité des classes de référence et incrémentale, les poids sont donc calculés à partir de décompositions singulières des matrices d’interpolation sous–jacentes. L’algorithme proposé combine alors les II classiques lorsqu’elles sont efficaces et les nouvelles II basées sur l’équation fondamentale pour les grilles les plus grossières. Les résultats ont été prouvés comme significatifs, car ils attestent de l’efficacité de la nouvelle stratégie. L’extension à un problème tridimmensionnel ne pose aucun problème particulier sinon que technique et est un travail en cours. Le troisième thème porte sur l’étude des systèmes d’EDP hyperboliques mises sous forme non conservative. Nous considérons plus précisément des systèmes non conservatifs provenant de l’élastodynamique et de la dynamique des gaz. Tout d’abord, la forme non conservative souvent intéressante pour le calcul de solutions approchées pose des problèmes d’ordre théorique. En fait, à cause de solutions physiquement admissibles présentant souvent des discontinuités, la forme faible des équations fait intervenir des produits de distributions qui n’ont pas de sens dans l’espace des distributions. Plusieurs théories ont été développées pour pallier cette difficulté, et c’est à l’aide de l’espace des fonctions généralisées de Colombeau que nous avons choisi de travailler. Selon cette théorie, il est possible de définir une relation plus faible que l’égalité, appelée l’association, entre fonctions appartenant à une classe de fonctions indéfiniment dérivables. Pour les besoins de notre étude c’est l’algèbre simplifiée de Colombeau (quotient algébrique de fonctions indéfiniment différentiables par un idéal de fonctions à croissance lente) qui a été retenue. Lors d’une première étude le système modélisant l’élasticité non linéaire en coordonnées eulériennes en 1D a été considéré. Pour celui-ci, le problème de Riemann a été résolu en utilisant une approche par pas fractionnaires où les termes de propagation ont été regroupés au sein d’un premier système hyperbolique, tandis que les termes convectifs ont été regroupés au sein d’un système non linéaire non strictement hyperbolique. Des approches de type Godunov en deux phases ont été proposées, en substituant la solution discontinue dans une équation de transport à coefficient discontinu ou encore en considérant le problème convectif dégénéré comme limite de problèmes bien posés paramétrés. Alors que le problème de propagation, présentant deux discontinuités de contact, est résolu par projection classique.
Une autre étude propose la résolution d’un des problèmes de la dynamique des gaz, parmi ceux qui ont le plus motivé la communauté, en formulation non conservative. Il s’agit du système d’équations d’Euler que nous avons considéré en formulation vitesse–pression–volume spécifique. Toujours par l’approche de pas fractionnaires, on regroupe les termes convectifs au sein d’un système non linéaire non strictement hyperbolique. Et le second système regroupant les termes de propagation est encore non linéaire mais hyperbolique. Les approches développées pour l’équation précédente peuvent encore s’appliquer. Et les résultats numériques du tube à choc de Sod attestent d’un bon comportement du schéma numérique dans plusieurs configurations.