Driss Yakoubi
Post Doctorant à REO Team
INRIA, Centre de Recherche de Paris Rocquencourt.

La modélisation de la respiration chez l'homme nécessite une bonne connaissance de la géométrie du poumon humain et de sa mécanique. Il est globalement constitué d'un arbre bronchique dyadique à 23 générations. Cet arbre est enveloppé du parenchyme qui se comporte comme un milieu visco-élastique troué par des alvéoles. Lors de l'inspiration, l'arbre bronchique est soumis à une force musculaire générant une augmentation du volume pulmonaire tandis qu'à l'expiration, il retourne à sa position d'équilibre de part sa nature élastique.
Le but de mon travail de recherche au cours de mon Post Doc au sein de l'équipe REO de l'INRIA, centre de recherche de Paris Rocquencourt est la modélisation des écoulements de l'air dans l'appareil respiratoire et de déterminer des algorithmes efficaces relatifs au modèle.
Le mod&aecute;le utilis&aecute consiste en une description selon trois parties de l'arbre respiratoire.
L'arbre bronchique a une géométrie très complexe et la partie distale (à partir de la 9ème génération) n'est pas visible par les technologies de l'imagerie médicale. Par conséquence, il est nécessaire d'élaborer un modèle mathématique simple et réalistique afin de fournir une meilleure compréhension des différentes pathologies du poumon humain.

Pour modéliser le couplage océan/atmosphère, nous avons utilisé un modèle de turbulence de type Reynolds Averaged Navier-Stokes (RANS), où les deux fluides (air, mer) sont couplés par une condition non linéaire à l'interface, et les inconnues vitesse-pression-énérgie cinétique turbulente (ECT) sont couplées par les viscosités turbulentes et par un terme source.
Nous avons proposé un schéma itératif pour approcher la solution de ce modèle. Ainsi, on montre qu'en présence des viscosités turbulentes suffisament grandes et pour des vitesses et ECT assez régulières, l'algorithme considéré est contractant. Nous montrons aussi que la limite de la solution du problème discret est une solution du problème continu.
Sous les mêmes conditions de convergence, nous avons établi que la solution du modèle est unique. Enfin, nous avons réalisé des simulations en méthodes spectrales et des éléments finis, montrant que la convergence de notre schéma ne dépend pas de la méthode de discrétisation choisie.

Babuska et Brezzi ont montré (séparément) que l'existence de la constante de la condition inf-sup est nécessaire pour obtenir l'existence de la pression dans les équations de Stokes et de Navier-Stokes. Plusieurs travaux ont été réalisés afin d'établir l'existence de cette constante dans différents domaines.
D'autre part, lorsqu'on considère une discrétisation du problème de Stokes ou de Navier-Stokes, la constante de la condition inf-sup dans les espaces discrets, dépend de la constante de la condition inf-sup "cas continu". Rappelons par exemple, que si on résout le problème discret en utilisant l'algorithme d'Uzawa, le conditionnement de la matrice et le nombre des itérations du gradient conjugué dépendent fortement de cette constante. Pour toutes ces raisons, nous nous sommes intéressé à sa caractérisation.
Nous avons trouvé un encadrement de cette constante, qui ne dépends que du relèvement harmonique sur le domaine du travail.

Les méthodes spectrales telles introduites par D. Gottlieb et S. Orszag, reposent sur l'approximation de solutions d'équations aux dérivées partielles initialement par des séries de Fourier tronquées puis par des polynômes de haut degré et sur l'utilisation de bases tensorisées des espaces d'approximation. Pour ces raisons, le domaine de base de ces méthodes est construit par tensorisation~(rectangle, parallélépipéde rectangle). Toutefois l'extension à un certain nombre d'autres géométries par transformation et par décomposition de domaines est maintenant bien connue, aussi bien en dimension 2 qu'en dimension 3.
Ces méthodes sont de précision infinie, au sens où l'ordre de l'erreur entre les solutions des problèmes continus et discrets n'est limité que par la régularité de la solution continue.
On propose une nouvelle approche de cette méthode dans des domaines en géométrie plus complexe (non-tensorielle). Cette méthode s'appuie sur deux idées~: un traitement des conditions aux limites de Dirichlet par pénalisation frontière (méthode de Nitsche), et une approximation de la géométrie du domaine par un pavage en parallélépipèdes.
Nous avons prouvé que cette nouvelle méthode conserve les propriétés spectrales, à savoir, les bons ordres d'erreurs d'approximation polynomiale et les ordres des estimations a priori, en norme L2 et H1. En effet, si on désigne par N, le degré des polynômes de Legendre dans les trois direction, et par epsilon le paramètre de pénalisation, alors l'erreur en norme L2 entre la solution exacte et la solution discréte est la somme de deux quantités : une esy dûe à la pénalisation et la seconde provient de des estimations a priori dans le cas tensoriel (donc en fonction du degré N).
Cette méthode à été validée par des tests numériques en utilisant le logiciel FreeFEM3D (voir rubrique simulations numériques ).

Nous avons implémenté la méthode spectrale décrite ci-dessus, en C++ dans le logiciel FreeFEM3D. Nous avions un cahier des charges à respecter. En effet, il fallait intégrer cette méthode de sorte que :

• le passage des éléments finis à spectrale, soit transparent,
• le mélange des éléments finis et spectrale pour le même calcul, soit permis.
• Résolution des EDP dans des domaines tensoriels et non-tensoriels

Une stratégie a été mise en place :

• écriture d'un code 1D pour commencer, (la méthode spectrale est tensorielle!)
• Diviser pour régner :
  • implémentation d'opérateurs de bases, c-à-d, programmer les différentes intégrales provennant des formes linéaires et bilinéaires:( Int(f*v), Int( mu grad u .grad v), etc... )
  • Nous avons mis en place des tests de non-régression.

Conception C++

• Factorisation du code tensoriel et non tensoriel
  nous avons tiré partie de la structure objet du langage C++, qui nous a permis par exemple d'intégrer de manière transparente de nouveaux types de fonctions et de maillages dans FreeFEM3D. Ainsi, nous avons utilisé un modèle de programmation objet : la méthode spectrale dans les domaines tensoriels et celle dans les domaines complexes utilisent la même brique de base. Il est donc souhaitable d'utiliser les mêmes lignes de code, afin d'éviter la réécriture du code, donc d'éventuelles erreurs, et de permettre qu'une optimisation unique profite aux deux méthodes.
• Interopérabilité Eléments finis et Spectrale
  Le logiciel FreeFEM3D est un code permettant la résolution des EDP par la méthode des éléments finis, donc les seules fonctions discrètes étaient de type éléments finis (femfunction dans le langage utilisateur). L'intégration de la méthode spectrale dans FreeFEM3D nous a conduit à créer un nouveau type de fonctions : les fonctions spectrales notées sfunction dans les programmes FreeFEM3D. Ainsi, la classe SpectralFunction hérite de la classe de base ScalarFunction, qui permet de manipuler de manière abstraite tout type de fonction scalaire.

Tests de validation
Nous avons validé cette méthode en réalisant plusieurs simulations académiques, (voir rubrique simulations numériques ).