Méthodes numériques pour l’écoulement de Stokes 3D

Modèle à viscosité variable

     Les mucus ne sont pas des fluides homogènes : ils sont essentiellement composés d’eau, comme tous les fluides  biologiques, ainsi que de plusieurs agents : protéines, particules inhalées, diverses cellules… Parmi ces agents, une protéine joue un rôle clé : la mucine. Elle est sécrétée par les cellules gobelets (elles aussi présentes sur la paroi bronchique) et est responsable de l’augmentation de la viscosité du fluide – voir figure 1.2 pour une illustration de l’écoulement de mucus autour des cils des cellules épithéliales (cette description sera reprise plus en détails au début du chapitre 6). À la surface, le film a tendance a être déshydraté par l’air qui circule dans les bronches, cela augmente localement la concentration de mucines et donc la viscosité. Certains auteurs considèrent d’ailleurs deux couches de fluide : la couche inférieure est appelée fluide périciliaire et la couche supérieure plus visqueuse garde l’appellation de mucus. Ce choix de modélisation n’a pas été fait dans ce travail pour deux raisons. Compte tenu de la composition de ces deux fluides (mélange d’eau et de protéines), l’existence d’une interface (et donc de deux couches distinctes) n’est pas certaine. Dans l’hypothèse où elle existerait, elle serait extrêmement difficile à identifier et il n’y a pas de mesures expérimentales qui permettent d’ajuster convenablement ce modèle. Deuxièmement, même au sein de chacune de ces couches la composition du fluide reste non homogène et il semble que la transition entre un fluide périciliaire très aqueux et le mucus très visqueux soit continue. C’est sur cette non homogénéité que l’étude se concentrera. Suite à des études expérimentales (détaillées plus loin dans ce manuscrit), d’autres auteurs ont cherché à modéliser les effets non newtoniens de cet écoulement. Bien que plusieurs équipes indépendantes aient fait plusieurs campagnes de mesure, l’identification de la loi de comportement du mucus reste un problème ouvert. Cela est dû notamment à la difficulté de faire des mesures expérimentales sur des matériaux vivants qui se dégradent dès qu’ils ne sont plus in-vivo. Dans ce travail on s’intéressera notamment à l’influence de la variation de viscosité des fluides biologiques sur les écoulements. En particulier, on ne considérera qu’une seule phase de fluide. La viscosité µ dépendra directement de la fraction massique α d’un certain agent responsable de ces variations de viscosité, dans le cas des mucus il s’agira de la fraction massique de mucines. Ainsi on introduira une loi du type : µ(t, x) = Φ (α(t, x)) (1.1) où α évoluera avec l’écoulement. Cette évolution est modélisée par une équation de convection-diffusion dans laquelle la convection domine car les mucines sont peu solubles dans l’eau. Ce modèle est un bon compromis entre un modèle à viscosité constante et l’utilisation de lois de rhéologie non newtoniennes dont on ne connaît pas encore les domaines de validité pour ce type de fluide. Les développements effectués dans ce sens ne ferment pas la porte à un ajout ultérieur d’une loi de rhéologie non newtoniennes, comme par exemple µ(t, x) = Φ (α(t, x),D(u)) (avec u la vitesse du fluide et D le tenseur des déformations.

Méthodes de résolution pour les équations de convection

     La résolution des équations de convection est un problème récurent en mécanique des fluides, il peut s’agir de quantités scalaires (conservation de la masse par exemple) ou vectorielles (terme d’advection dans Navier-Stokes). Contrairement aux équations paraboliques ou elliptiques, il n’y a pas de termes régularisants et des méthodes robustes, stables et non oscillantes doivent être utilisées pour faire des calculs précis. Deux approches peuvent être choisies pour résoudre ce type d’équations : une approche eulérienne utilisant une grille de discrétisation ou une approche lagrangienne utilisant une discrétisation particulaire. L’utilisation de particules n’est pas une idée récente puisque dans [Rosenhead, 1931] des particules sont introduites et les calculs sont faits analytiquement pour étudier la déformation d’une feuille de vorticité. Son utilisation en vue d’applications numériques a été reprise par [Harlow, 1962] qui a introduit la technique Particle-in-Cell. Cette technique a reçu de nombreuses améliorations depuis ce travail précurseur, et a notamment été utilisée intensément dans la communauté des méthodes vortex par les équipes de Anderson, Cottet, Leonard, Koumoutsakos, Winckelmans… Deux ouvrages rassemblent une bonne partie de ces résultats [Cottet et Koumoutsakos, 2000; Koumoutsakos et al., 2009]. Le cadre particulaire permet de transformer la résolution d’une équation de convection en un système d’équations différentielles posé sur des particules. Cela supprime les conditions de stabilité de type CFL. Les résultats de convergence datant des années 1980 [Anderson et Greengard, 1985; Beale et Majda, 1981, 1982; Cottet et Mas-Gallic, 1990] ont montré que les constantes devant les termes d’erreurs sont faibles. Cela permet de faire de grands pas de temps en restant stable et précis. De plus le coût de cette méthode est linéaire par rapport au nombre de particules ce qui la rend très attractive pour des systèmes de grande dimension. D’autres travaux se sont intéressés à la résolution des opérateurs de diffusion directement sur les particules [Degond et Mas-Gallic, 1989a,b] ou bien en développant un stratégie hybride grille particules [Couet et al., 1981] performante en utilisant des opérateurs d’interpolation d’ordre élevé et à support compact entre les deux discrétisations. Le problème des méthodes particulaires réside dans la déformation du nuage de particules avec le champ de vitesse, ce qui peut induire des pertes de recouvrement (c’est à dire des zones d’accumulation ou de raréfaction de particules). Dans ce cas, de l’information est perdue et la solution est faussée. Il faut donc périodiquement réinitialiser les particules pour garder le recouvrement. Cela est très naturel dans le cas hybride grille-particules grâce aux interpolations. Ces méthodes ont connus des extensions pour traiter des problèmes avec un maillage adaptatif [Bergdorf et al., 2005], de la décomposition de domaine[Cottet, 1990], des simulations multi-échelles [El Ossmani et Poncet,2010], des problèmes à plusieurs milliards de degrés de liberté [Chatelain et al., 2008], une implémentation sur GPU [Rossinelli et Koumoutsakos, 2008], parmi d’autres exemples…

Méthodes lagrangiennes basées sur les solutions fondamentales

    Enfin un quatrième type de méthode d’interaction fluide structure repose sur le développement des solutions en fonction des solutions fondamentales. Cette méthode adaptée aux écoulements de Stokes exploite la linéarité du problème lorsque la viscosité est constante. Cette méthode appelée Stokeslets a été introduite dans [Cortez, 2001] en 2D et dans [Cortez et al., 2005] en 3D. La méthode repose sur le calcul de point de forces répartis dans les obstacles ou sur leurs interfaces à partir des solutions fondamentales du problème de Stokes à viscosité constante, basées sur les fonctions de Green. La méthode a été adaptée à des conditions aux limites périodiques [Leiderman et al., 2013]. Deux arguments s’opposent à l’utilisation de cette méthode dans ce travail. Tout d’abord la méthode n’est valable qu’à viscosité constante ce qui n’est pas compatible avec le couplage non linéaire décrit précédemment. Deuxièmement, le coût de calcul dépend du nombre de points de force créés (en résolvant un système linéaire plein) et la méthode est limitée lorsqu’on s’intéresse à des géométries trop lourdes.

Résolution de la diffusion de la fraction massique

    Grâce à la méthode de splitting introduite dans le paragraphe 2.1.1, il a été possible de séparer la résolution de la convection et de la diffusion de la fraction massique. La convection est discrétisée de manière lagrangienne avec une méthode Particle-in-Cell, comme cela a été présenté dans la section précédente. La résolution de la diffusion peut être faite de manière explicite : le phénomène dominant dans cette équation est la convection et l’explicitation du terme diffusif ne contraint pas la la stabilité de la méthode. Pour cette discrétisation un schéma aux différences finies standard à sept points est utilisé. Le principal argument qui justifie ce choix est bien entendu sa facilité d’implémentation sur la grille cartésienne régulière. Cela se combine très bien avec la méthode d’intégration temporelle décrite dans le paragraphe 2.1.2 puisque l’effet régularisant de ce terme de diffusion est recherché avant chaque calcul du champ de vitesse u, la solution du problème de Stokes. Ce calcul est également effectué sur la grille cartésienne et il n’y a pas besoin d’ajouter un transfert entre la grille et les particules. Un autre choix possible aurait été d’utiliser une méthode dite “Particle Strength Exchange” (PSE) introduite par Degond et Mas-Gallic à la fin des années 1980 [Degond et Mas-Gallic, 1989a,b]. Cette méthode permet de résoudre les opérateurs de diffusion directement sur le nuage particulaire, ce qui permet de s’affranchir des interpolations récurrentes entre la grille et les particules (mais pas des problèmes liés au remaillage particulaire, voir section 2.1.5). Une extension de cette méthode pour des opérateurs de dérivées spatiales quelconque a été introduite dans [Eldredge et al., 2002b], qui ont été utilisées dans la communauté des méthodes vortex ([Eldredge et al., 2002a] par exemple). Ces méthodes approchent l’opérateur de diffusion par une formule de quadrature faisant intervenir un noyau discret de diffusion. En pratique, il s’agit de faire interagir (à travers ce noyau de diffusion discrétisé) la quantité portée par chaque particule avec les particules voisines. La taille du voisinage est directement liée à l’importance du coefficient de diffusion. Bien que la méthode PSE soit un peu plus précise qu’une méthode aux différences finies standard [Poncet, 2006], il a été décidé de ne pas utiliser cette méthode pour le calcul de la diffusion de la fraction massique. Le premier argument étant que comme ce phénomène n’est pas dominant il n’est pas nécessaire de le calculer avec la méthode la plus robuste. Le second argument a déjà été avancé : compte tenu du splitting et de la méthode d’intégration temporelle que nous avons développé, une méthode explicite de diffusion sur grille structurée s’adapte très bien et est peu coûteuse ; la méthode PSE est surtout très efficace sur des grilles non structurées où aucune hypothèse d’agencement des particules peut être faite, ce qui interdit l’utilisation de méthodes de différences finies et de solveurs rapides. Enfin, la méthode PSE est plus difficile et plus coûteuse à implémenter puisqu’elle demande des calculs de voisinages qui peuvent s’avérer très onéreux en 3D s’ils ne sont pas méticuleusement développés, en se basant sur des algorithmes de tri rapide notamment.

Modélisation de la rhéologie du mucus

   Les études rhéologiques qui ont été réalisées sur le mucus pulmonaire tendent à s’accorder sur le fait que ce n’est pas un fluide newtonien. Cependant de fortes variations existent entre ces résultats. Certains auteurs pensent qu’il s’agirait plutôt d’un fluide viscoélastique [Puchelle et al., 1987; Mitran, 2007] à mémoire de forme alors que d’autres études ont conclu à un comportement viscoplastique [Edwards et Yeates, 1992; Craster et Matar, 2000]. Dans une revue récente, les auteurs ne tranchent pas sur ce point [Lai et al., 2009]. De nouvelles campagnes de mesure de la rhéologie du mucus, qui viendront par la suite enrichir le modèle, sont prévues dans le cadre du projet ANR BioFiReaDy . Ce point sera abordé dans les perspectives de ce travail. Il a été observé que le mucus est plus hydraté au voisinage de la paroi bronchique qu’à la surface du film où l’air a tendance à le déshydrater. La fraction de mucines est donc plus importante à cet endroit. La partie inférieure du film, la plus hydratée, est appelée fluide périciliaire, voir figure 6.1. Certains travaux numériques considèrent un modèle à deux fluides : une couche de fluide périciliaire en dessous d’une couche plus visqueuse riche en mucines [Smith et al., 2007b]. La principale difficulté dans ce modèle consiste à régler la hauteur (et les lois de comportement) de l’interface entre les deux phases ; cette interface est difficile à identifier compte tenu du fait que les deux fluides sont essentiellement composés d’eau. Le modèle à viscosité variable est une bonne alternative à ces problèmes puisqu’il n’impose pas de suivit de cette interface, tout est piloté par la fraction massique. Dans ce travail, nous avons décidé de ne pas étudier les effets non newtoniens en raison de cette incertitude dans le choix du modèle à utiliser. Ce point sera repris dans les perspectives, dans le dernier chapitre de ce manuscrit. Notre étude s’est focalisée sur le problème 3D en étudiant l’influence de la viscosité variable, avec une description précise de la géométrie ciliaire. Ce point n’ayant pas encore été abordé dans la littérature.

Table des matières

Chapitre 1 Introduction 
1.1 Position du problème
1.1.1 Contexte
1.1.2 Problématique
1.1.3 Choix de modélisation
Modèle à viscosité variable
Généralités sur le problème de Stokes
Interaction entre le fluide et la géométrie en mouvement
Interaction entre le fluide et le cycle respiratoire
1.2 État de l’art, choix des méthodes de discrétisation et de résolution
1.2.1 Méthodes de résolution pour les équations de convection
1.2.2 Méthodes de résolution pour l’interaction fluide-structure
Maillage conforme
Méthodes de frontières immergées
Méthodes de domaine fictif volumique
Méthodes lagrangiennes basées sur les solutions fondamentales
1.2.3 Méthodes de résolution pour le problème de Stokes
1.2.4 Méthodes de discrétisation
1.2.5 Applications aux fluides biologiques
1.3 Contributions
1.4 Organisation du manuscrit
Chapitre 2 Algorithmes numériques 
2.1 Méthodes de résolution pour l’équation de convectiondiffusion
2.1.1 Méthode de splitting pour le problème de convection-diffusion
Cas linéaire
Splitting du problème de convection-diffusion
2.1.2 Intégration temporelle de la méthode de splitting
2.1.3 Discrétisation lagrangienne de la convection
2.1.4 Discussion sur les noyaux d’interpolation
Choix des noyaux
Implémentation
Calcul parallèle
2.1.5 Remaillage particulaire
2.1.6 Résolution de la diffusion de la fraction massique
2.2 Algorithme rapide pour le problème de Stokes pénalisé
2.2.1 Méthode de pénalisation pour l’interaction fluide-structure
2.2.2 Une méthode de projection sur les champs à divergence nulle
2.2.3 Méthode de point fixe multi-critère
Présentation de la méthode
Adaptation à d’autres conditions aux limites
2.2.4 Traitement explicite des termes de viscosité non homogène
2.2.5 Implémentation numérique
Solveurs utilisés
Projection et grilles décalées
Calcul parallèle
2.3 Résumé de l’aglgorithme complet de résolution
Chapitre 3 Validations 
3.1 Tourbillon de Green-Taylor
3.1.1 Convergence de l’algorithme de point fixe
3.1.2 Ordre de la méthode
3.1.3 Temps de calcul
3.2 Cas test de la sphère dans un écoulement incident
3.3 Convergence en temps et discussion sur les intégrateurs temporels
3.4 Conservation de la masse
3.5 Temps de calcul pour le problème complet
Chapitre 4 Améliorations des algorithmes de résolution 
4.1 Approche Shermann-Morrison-Woodbury pour la pénalisation
4.1.1 La formule de Sherman-Morrison-Woodbury
4.1.2 Présentation pour un problème de Poisson pénalisé
Réécriture du problème pour faire apparaître la formule SMW
Approche itérative de type GMRES et couplage avec les solveurs FFT
Temps de calcul et occupation mémoire, comparaison avec Mudpack
4.1.3 Application au problème de Stokes pénalisé
Reformulation de l’étape d’estimation
Étude des performances en fonction de la géométrie
4.1.4 Extension aux équations de Navier-Stokes
4.2 Accélération de la convergence du point fixe
4.2.1 Préconditionnement multigrille
4.2.2 Extrapolation de Richardson des conditions aux limites
4.2.3 Améliorations des temps de calcul
4.3 Régularisation de la fonction caractéristique
4.4 Raffinement doublé des variables transportées
Chapitre 5 Résultats d’existence et d’unicité pour le problème pénalisé 
5.1 Notations et définitions préliminaires
5.2 Énoncé du résultat principal
5.3 Arguments techniques
5.4 Analyse du problème de Stokes pénalisé
5.4.1 Existence et unicité
5.4.2 Prolongement “harmonique” de u
5.4.3 Estimations sur (u, p)
5.4.4 Régularité de la solution
5.5 Problème de Stokes pénalisé dans un domaine dépendant du temps
5.6 Résultats pour le problème de convection-diffusion
5.7 Étude du couplage entre le problème de Stokes pénalisé et le problème de convection-diffusion
5.8 Preuve du théorème 2
5.9 Preuve du théorème 3
5.9.1 Étude de α
5.9.2 Étude de µ
5.9.3 Étude de (u, p)
5.9.4 Unicité
Chapitre 6 Applications aux fluides biologiques
6.1 Généralités sur les fluides biologiques
6.2 Transport du mucus pulmonaire par les cils des cellules épithéliales
6.2.1 Position du problème
6.2.2 Modélisation de la rhéologie du mucus
6.2.3 Dynamique d’oscillation de la cellule épithéliale
6.2.4 Conditions aux limites pour l’écoulement du fluide
6.2.5 Simulation autour d’un cil
6.2.6 Simulations autour de différentes cellules
Cellule à seize cils
Cellule de forme elliptique et battement circulaire
Tapis de cils et synchronisation métachronale
6.2.7 Influence des paramètres biologiques
Profil de viscosité linéaire à un paramètre
Profil de viscosité logistique à trois paramètres
Fréquence de battement
Épaisseur du film de mucus
Diffusion de la fraction massique
Synchronisation métachronale
Conclusion sur l’influence des différents paramètres
6.3 Exemple d’un micro-nageur : le spermatozoïde
6.3.1 Contexte et position du problème
6.3.2 Modélisation du mouvement de battement
6.3.3 Simulations numériques
6.4 Calculs d’écoulements dans les milieux poreux
6.4.1 Position du problème
6.4.2 Transport dans les milieux poreux
Chapitre 7 Conclusions et perspectives 
7.1 Conclusions
7.1.1 Modélisation
7.1.2 Développement d’algorithmes rapides de résolution
Méthodes utilisées
Performances
7.1.3 Modélisation de l’écoulement du mucus pulmonaire
7.1.4 Autres applications
7.2 Perspectives
7.2.1 Enrichissement du modèle
Modèle two-ways pour l’interaction fluide-structure
Prise en compte de l’interface air-liquide
Modélisation des effets non newtoniens
Modélisation des sécrétions de mucines
7.2.2 Plate-forme intégrative
7.2.3 Analyse mathématique
7.2.4 Extension à d’autres écoulements
Bibliographie

Télécharger le rapport complet