<figure><img src="image_1.png" /></figure>

HAL
open science

# Méthodes hybrides grille-particules pour la simulation 3D d'écoulements complexes non newtoniens

Robin Chatelin, Philippe Poncet

► To cite this version:

Robin Chatelin, Philippe Poncet. Méthodes hybrides grille-particules pour la simulation 3D d'écoulements complexes non newtoniens. 12e colloque national en calcul des structures, CSMA, May 2015, Giens, France. hal-01400458

HAL Id: hal-01400458

https://hal.science/hal-01400458v1

Submitted on 22 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Méthodes hybrides grille–particules pour la simulation 3D d’écoulements complexes non newtoniens

Robin Chatelin<fnref n="1" />, Philippe Poncet<fnref n="2" />

**Résumé** — Quelques méthodes numériques récemment développées dans le cadre du projet ANR Bio-FiReaDy sont regroupées dans cette communication. L'objectif est de simuler en 3D des écoulements de fluides complexes interagissant avec des obstacles immergés. L'algorithme repose sur une discrétisation hybride grille-particules, des algorithmes à pas fractionnaires et l'utilisation de solveurs rapide afin de réduire les coûts de calcul. Cela ouvre la porte à des simulations 3D de grande dimension, tout en permettant d'étudier l'influence de nombreux jeux de paramètres sur l'écoulement.

**Mots clés** — Fluides complexes, Méthodes numériques pour les écoulements 3D, Méthodes hybrides grille–particules, Interaction fluide–structure, Calcul scientifique

Cet exposé se concentre sur la présentation de stratégies algorithmiques adaptées au calcul 3D d'écoulements complexes, très visqueux et incompressibles. Le fluide considéré n'est pas homogène et ses propriétés dépendent de la fraction massique $\alpha$ d'un certain agent. Il peut en particulier avoir des propriétés non newtoniennes de type loi de puissance. De plus, le fluide interagit avec des obstacles mobiles immergés, ce qui conduit à la résolution du problème non linéaire (1).

$$ 
\left\{ 
\begin{array}{ll} 
-\mathrm{div}(\tilde{\sigma}(u)) + \nabla p = f & \text{dans } \Omega \setminus \mathcal{B}(t) \\ 
\mathrm{div}u = 0 & \text{dans } \Omega \setminus \mathcal{B}(t) \\ 
u(t) = \bar{u}(t) & \text{dans } \mathcal{B}(t) \\ 
\displaystyle\frac{\partial \alpha}{\partial t} + u \cdot \nabla \alpha - \eta \Delta \alpha = 0 & \text{dans } \Omega 
\end{array} 
\right. 
\quad (1) 
$$

$$ \text{où } \tilde{\sigma} = K(\alpha)\dot{\gamma}^{n(\alpha)-1}D(u), \quad D(u) = \frac{\nabla u + \nabla u^T}{2} \quad \text{et} \quad \dot{\gamma} = 2\|D\| = 2\sqrt{D:D} $$

$\Omega$ est le domaine de calcul complet, $\mathcal{B}(t)$ correspond au domaine occupé par les obstacles immergés à l'instant $t$, $\Omega \setminus \mathcal{B}(t)$ correspond donc au domaine occupé par le fluide à cet instant. Dans ce système $u$ est le champs de vitesse, $p$ la pression, $\sigma = \tilde{\sigma} + pId$ le tenseur des contraintes, $D$ le tenseur des déformations et $\bar{u}(t)$ la vitesse des obstacles à l'instant $t$. Enfin $K$ et $n$ sont les paramètres de la loi rhéologique qui peuvent dépendre de la fraction massique $\alpha$, transportée par l'écoulement et qui peut également se diffuser (en pratique $\eta$ reste "petit").

La méthode de pénalisation [1] permet de gérer l'interaction entre le fluide et les obstacles. Elle fait partie de la grande famille des méthodes de "domaine fictif" pour lesquelles il n'est pas nécessaire de connaître précisément la position de l'interface entre le fluide et les obstacles (dont la position peut être très délicate à calculer). Le premier avantage de ce type de méthodes est que le maillage n'a pas besoin d'être conforme aux obstacles, il n'est donc pas nécessaire de le déplacer avec des méthodes de type ALE qui sont très coûteuses en 3D. Le second est que la première équation de (1) est étendue dans tout le domaine de calcul $\Omega$ en forçant l'égalité $u(t) = \bar{u}(t)$ dans $\mathcal{B}(t)$ :

$$ 
\left\{ 
\begin{array}{ll} 
-\mathrm{div}(\tilde{\sigma}(u)) + \nabla p + \frac{\chi}{\epsilon}(u - \bar{u}) = f & \text{dans } \Omega \\ 
\mathrm{div}u = 0 & \text{dans } \Omega \setminus \mathcal{B}(t) \\ 
\displaystyle\frac{\partial \alpha}{\partial t} + u \cdot \nabla \alpha - \eta \Delta \alpha = 0 & \text{dans } \Omega 
\end{array} 
\right. 
\quad (2) 
$$

<footnote n="1">¹ Université de Lyon, ENI Saint-Étienne, LTDS, UMR 5513 CNRS, 58 rue Jean Parot, 42023 Saint-Étienne, France, robin.chatelin@enise.fr</footnote>

<footnote n="2">² Université de Pau et des Pays et l'Adour, LMAP, UMR 5142 CNRS, IPRA, avenue de l'Université, 64013 Pau, France, philippe.poncet@univ-pau.fr</footnote>
avec $\chi$ la fonction caractéristique de $\mathcal{B}(t)$ et $\varepsilon \ll 1$ le paramètre de pénalisation. La vitesse des obstacles $\bar{u}(t)$ est supposée connue a priori. L'interaction fluide-structure complète peut être étendue dans le cas de solides rigides [2] ou élastiques [3].

Il a été montré [4] que le problème (2) admet une unique solution et elle converge vers la solution de (1) lorsque $\varepsilon \to 0$.

Remarque : bien qu'il s'agisse d'un problème de Stokes quasi-statique, la première équation de (2) dépend implicitement du temps à travers $\alpha$ (qui intervient dans $\sigma$), $\bar{u}$ et les conditions aux limites (non précisées ici pour alléger les notations).

Pour résoudre l'équation de convection-diffusion sur $\alpha$, une méthode de splitting est utilisée [5, 6]. Cela sépare la résolution de la convection et de la diffusion, permettant ainsi d'utiliser une discrétisation bien adaptée pour chacun de ces phénomènes. Cela conduit à la résolution d'une équation de convection, puis d'une équation de la chaleur.

La fraction massique $\alpha$ est d'une part discrétisée sur un nuage de $N$ particules [6, 7], chacune étant localisée au point $\xi_p$, portant une fraction massique $\alpha_p$ et ayant un volume $v_p$ constant (car l'écoulement est incompressible) : $\alpha = \sum_{p=1}^{N} \alpha_p \delta_{\xi_p} v_p$

L'équation de convection, discrétisée avec cette méthode particulière, est alors transformée en un système d'équations différentielles. L'avantage de ce système est qu'il peut être résolu de manière explicite sans être contraint par une condition de stabilité de type CFL, de grands pas de temps peuvent ainsi être utilisés :

$$ \frac{d\alpha_p}{dt} = 0, \quad \frac{d\xi_p}{dt} = u(\xi_p) $$

Une discrétisation sur une grille est par ailleurs utilisée pour la résolution de la diffusion (équation de la chaleur). Comme le coefficient de diffusion reste petit (la convection est le phénomène dominant dans les applications visées) une résolution explicite ne pose pas de problème particulier.

Le coût de calcul de ces méthodes est linéaire par rapport au nombre de particules, ce qui les rendent très attractives pour des calculs 3D de grande dimension. Le point clé pour ne pas ruiner cette efficacité est de faire des interpolations efficaces entre la grille et le nuage de particule. L'utilisation de noyaux à support compact [8] permet de conserver ce coût de calcul qui fait la force des méthodes hybrides grille-particules.

Afin d'exploiter pleinement ces performances de calcul il reste alors à calculer efficacement le champs de vitesse $u$, solution du problème de Stokes. Pour alléger les notations on note $\mu_{NN}(\alpha, u) = K(\alpha)\dot{\gamma}^{n(\alpha)-1}$ qui peut être vu comme une “viscosité non newtonienne équivalente”.

Les méthodes de projection sont traditionnellement utilisées pour la résolution des équations de Navier-Stokes instationnaires [9] : ce sont des algorithmes à pas fractionnaire, qui permettent de séparer le calcul du champs de vitesse et de la pression dans un algorithme de type prédiction/correction. Une méthode similaire a été développée pour ce problème quasi-statique à viscosité variable [10]. Il s'agit d'une méthode itérative qui permet de traiter de manière explicite, dans un point fixe multi-critères, les termes issus de la non homogénité du fluide et la non linéarité des termes non newtoniens. La méthode garanti également une bonne précision du champs de vitesse proche des bords (des obstacles et du domaine de calcul), problème inhérent aux méthodes de projection traditionnelles [9].

À chaque itération $i$, l'algorithme s'écrit en trois étapes (voir [10, 11] pour plus de détails), il est présenté ici pour des conditions aux limites de Dirichlet $u = g$ sur $\partial\Omega$ :

1. Prédiction, résoudre le problème :

$$ \begin{cases} -\mu_{NN}(\alpha, u_i)\Delta u_{i+1}^* + \frac{\chi}{\varepsilon}(u_{i+1}^* - \nabla\zeta_i - \bar{u}) = f + [2D(u_i) + (\text{div}u_i^*)\text{Id}]\nabla\mu_{NN}(\alpha, u_i) & \text{dans } \Omega \\ u_{i+1}^* = g - \nabla\zeta_i & \text{sur } \partial\Omega \end{cases} $$

2. Calcul du projecteur $\zeta$ en résolvant le problème de Poisson : $$ \begin{cases} -\Delta\zeta_{i+1} = -\text{div}u_{i+1}^* & \text{dans } \Omega \\ \frac{\partial\zeta_{i+1}}{\partial n} = 0 & \text{sur } \partial\Omega \end{cases} $$
3. Correction : $u_{i+1} = u_{i+1}^* - \nabla \zeta_{i+1}$

Une dicrétisation avec une grille cartésienne régulière est utilisée, ce qui permet d'utiliser des sol-
veurs rapides pour la résolution du problème de Poisson [12] à l'étape 2. Le coût de calcul est quasi
linéaire par rapport au nombre de points de grille, i.e. $O(N \log N)$.

Pour la résolution du problème de l'étape 1, un algorithme basé sur la formule de Shermann-Morrison-Woodbury a été développé [13]. Il consiste à voir le problème de l'étape 1 comme un problème de Poisson perturbé du type $-\mu\Delta u + cu = g$ (où $c$ est à support dans $\mathcal{B}$). Cet algorithme permet de résoudre un autre problème itératif posé sur les points pénalisés (i.e. situés dans les obstacles) et qui est donc de dimension bien inférieure au problème initial. Il se décompose en une suite de problèmes de Poisson, ce qui permet de profiter de l'efficacité du solveur précédent et d'avoir un algorithme complet au coût de calcul quasi linéaire, particulièrement bien adapté aux calculs 3D de grande dimension.

Les algorithmes de résolution présentés ci-dessus ont été implémentés dans un code de calcul qui a principalement été utilisé pour simuler l’écoulement du mucus pulmonaire. L’objectif est d’étudier le phénomène de clairance mucociliaire qui permet d’évacuer les impuretés inhalés (agents pathogènes,

<figure><img src="image_3.png" /><figcaption>FIGURE 1 – (En haut) Simulation de l’écoulement du film de mucus autour d’un tapis de 128 cils avec battement asynchrone (Résolution : 1024 × 512 × 256 points de grille). En vert est affichée une isosurface de fraction massique de mucines. L’échelle de couleur présente la norme de vitesse sur les bords du domaine (conditions aux limites périodiques). (En bas) Exemples d’écoulements d’un scalaire passif à travers deux milieux poreux (Résolution : 1024³ points de grille). En gris est représentée une isosurface de la fonction caractéristique de la géométrie du milieu. L’échelle de couleurs, utilisée pour représenter le scalaire, quantifie la norme de la vitesse de l’écoulement. (En bas à gauche) Échantillon d’un empilement de microsphères de silice. (En bas à droite) Vue multi-échelles de l’écoulement dans un échantillon de grès de Castlegate en Utah (seule la partie supérieur est affichée)</figcaption></figure>
poussières, pollution, ...). Ces simulations ont été réalisées dans le cadre du projet ANR BioFiReaDy<fnref n="1" />,
en se focalisant sur l'écoulement à l'échelle des cellules épithéliales, équipées de cils vibratiles, dont
le battement entraîne le film de mucus. Dans ce fluide biologique, α représente la fraction massique
de mucines dans le mucus, ce sont des protéines qui influencent directement la viscosité et la rhéolo-
gie [14]. L'objectif de ces simulations est d'étudier l'influence des différents paramètres du modèle sur
l'efficacité de la clairance mucociliaire, une stagnation du mucus entraînant de nombreuses infections
secondaires, notamment dans le cas de pathologies comme la mucovicidose. Grâce aux bonnes perfor-
mances algorithmiques présentées ci-dessus de nombreux jeux de paramètres peuvent être simulés de
manière systématique.

Ce code de calcul a également été utilisé pour simuler des écoulements en milieux poreux à l'échelle des pores (figure 1). Cet écoulement est modélisé par le problème de Stokes-Brinkman qui est exactement un problème de Stokes pénalisé (le paramètre de pénalisation quantifie ici la porosité de la phase solide). Les bonnes performances de calcul de ce solveur de Stokes permettent en outre de calculer les tenseurs de perméabilité équivalente du milieu à très haute résolution (testé jusqu'à un milliard de points de grille). Ils peuvent ensuite être utilisés dans le problème de Darcy, beaucoup plus simple à résoudre.

Remerciements

Ce travail a été partiellement financé par le projet ANR BioFiReaDy (ANR-2010-JCJC-0113-01).
Ce travail a bénéficié d'un accès aux ressources de calcul hautes performances CALMIP (P0314).

Références

[1] Philippe Angot, Charles-Henri Bruneau, and Pierre Fabrie. A penalization method to take into account obstacles in incompressible viscous flows. *Numerische Mathematik*, 81(4):497–520, 1999.

[2] M. Coquerelle and G.-H. Cottet. A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies. *Journal of Computational Physics*, 227(21):9121–9137, November 2008.

[3] Mattia Gazzola, Philippe Chatelain, Wim M. van Rees, and Petros Koumoutsakos. Simulations of single and multiple swimmers with non-divergence free deforming geometries. *Journal of Computational Physics*, 230(19):7093–7114, August 2011.

[4] Robin Chatelin, Davis Sanchez, and Philippe Poncet. Analysis of penalized variable viscosity 3d stokes equations coupled with diffusion and transport. *Submitted*, 2014.

[5] J. Thomas Beale and Andrew Majda. Rates of convergence for viscous splitting of the navier-stokes equations. *Mathematics of Computation*, 37(156):243–259, October 1981.

[6] Georges-Henri Cottet and Petros D. Koumoutsakos. *Vortex Methods: Theory and Practice*. Cambridge University Press, March 2000.

[7] G.-H. Cottet and P. Poncet. Advances in direct numerical simulations of 3D wall-bounded flows by vortex-in-cell methods. *Journal of Computational Physics*, 193(1):136–158, January 2004.

[8] Michael Bergdorf and Petros Koumoutsakos. A lagrangian Particle–Wavelet method. *Multiscale Modeling & Simulation*, 5(3):980–995, January 2006.

[9] J.L. Guermond, P. Minev, and Jie Shen. An overview of projection methods for incompressible flows. *Computer Methods in Applied Mechanics and Engineering*, 195:6011–6045, September 2006.

[10] Robin Chatelin and Philippe Poncet. A hybrid grid-particle method for moving bodies in 3D stokes flow with variable viscosity. *SIAM Journal on Scientific Computing*, 35(4):B925–B949, August 2013.

[11] Robin Chatelin. Méthodes numériques pour l'écoulement de Stokes 3D : fluides à viscosité variable en géométrie complexe mobile ; application aux fluides biologiques. PhD thesis, Université Toulouse 3 Paul Sabatier, 2013.

[12] Paul Swarztrauber and Roland Sweet. Efficient FORTRAN subprograms for the solution of elliptic partial differential equations (abstract). *SIGNUM Newsl.*, 10(4), December 1975.

[13] Robin Chatelin and Philippe Poncet. Hybrid grid-particle methods and penalization : A Sherman–Morrison–Woodbury approach to compute 3D viscous flows using FFT. *Journal of Computational Physics*, 269:314–328, July 2014.

[14] D.J. Smith, E.A. Gaffney, and J.R. Blake. Modelling mucociliary clearance. *Respiratory Physiology & Neurobiology*, 163:178–188, November 2008.

<footnote n="1">1. https://sites.google.com/site/biofiready/</footnote>