Code de calcul :
Simulation d'écoulements visco-élastiques 2D ou 3D avec le modèle d'Oldroyd par une méthode de splitting
Contexte :
Les équations de Navier-Stokes / Oldroyd-B permettent de modéliser des écoulements visco-élastiques. Les non-linéarités introduites dans le modèle pour décrire les propriétés rhéologiques de ces fluides non-newtoniens se traduisent mathématiquement et numériquement par des difficultés qui nécessitent un traitement spécifique. Dans ce contexte, nous nous sommes intéressés à plusieurs aspects de ce modèle :
Modèle diffusif et limite de viscosité évanescente. Le modèle d'Oldroyd possède une variante qui consiste à ajouter un terme de difusion dans l'équation constitutive de la loi élastique. Le lien entre le modèle diffusif et le modèle usuel est simple d'un point de vue formel mais l'analyse mathématique de l'asymptotique associée à une viscosité évanescente ne l'est pas du tout ; le calcul direct d'écoulements 2D ou 3D permet de conjecturer le comportement des solutions associées à cette étude asymptotique.
Interaction fluide-structure dans un fluide visco-élastique. La présence d'inclusions rigides dans un fluide se traduit par un couplage entre le fluide et les objets rigides immergés dans le fluide. La simulation numérique d'un tel système nécessite de mettre en œuvre des approches adaptées à la détermination de la dynamique de ce couplage. Les difficultés reposent sur la résolution des équations de la mécanique des fluides (équations de Navier-Stokes / Oldroyd-B dans le cas présent) au sein d'un domaine mobile en temps et à la détermination des interactions entre le fluide et l'inclusion rigide.
Algorithme :
La méthode de calculs est basée sur une discrétisation par éléments finis (P1b ou P2 pour le champ de vitesse, P1 pour le champ de pression, P0 pour le champ de contrainte élastique dans la version usuelle, P1 pour le champ de contrainte élastique dans la version diffusive) et un algorithme de splitting en temps :
L'interaction fluide-structure entre des objets rigides et fluide visco-élastique peut être facilement étudiée en intégrant une approche de type domaine fictif, par pénalisation du tenseur de déformation associé au champ de vitesse et du tenseur de contrainte élastique.
l'équation bilan de la quantité de mouvement et l'équation de continuité sont résolues par traitement explicite du terme de contrainte élastique ; il s'agit donc de résoudre un problème de type Navier-Stokes. La non-linéarité associée au terme de dérivée totale (décrivant les effets inertiels) est résolue par la méthode des caractéristiques.
l'équation constitutive est un système de transport(-diffusion). Dès lors, le champ de vitesse étant connu par l'étape précédente, le système de transport non-linéaire est résolu par une approximation associée à la méthode des caractéristiques.
Programmation :
Le programme, implanté en Freefem++, est constitué de plusieurs versions commentées, permettant de prendre en compte différentes géométries : 2D (cavité entraînée, écoulement de Poiseuille dans un tube contenant un obstacle), 3D (cavité entraînée), interaction fluide-structure en 2D (sédimentation de sphères rigides en 2D, cisaillement d'une collection de sphères ou ellipses rigides, avec algorithme de gestion de contacts). Les fichiers de résultats sont édités en format .vtk.
Auteurs :
Laurent Chupin (Clermont-Ferrand), Astrid Decoene (Orsay), Sébastien Martin (Paris) & Bertrand Maury (Orsay) | version originale : juillet 2014 | version augmentée : juillet 2017
Illustrations numériques :
Fig. 1 : Solution stationnaire associée à un écoulement visco-élastique dans tube avec obstacle (Re=4, We=1, r=0.5). Représentation du module du champ de vitesse, de la contrainte de cisaillement élastique et de la contrainte d'élongation élastique autour de l'obstacle.
 |  |  |
| a) Champ de vitesse | b) Contrainte de cisaillement élastique | c) Contrainte d'élongation élastique |
Fig. 2 : Solution stationnaire du probléme de la cavité entraînée 2D (Re=10, We=1, r=0.5). Représentation du module du champ de vitesse, de la contrainte de cisaillement élastique et de la contrainte d'élongation élastique.
 |  |  |
| a) Champ de vitesse | b) Contrainte de cisaillement élastique | c) Contrainte d'élongation élastique |
Fig. 3 : Solution stationnaire du probléme de la cavité entraînée 3D (Re=4, We=1, r=0.5).
 |  |  |
| Champ de vitesse | Cisaillement élastique σ13 | Cisaillement élastique σ23 |
 |  |  |
| Champ de vitesse | Élongation élastique σ11-σ33 | Élongation élastique σ22-σ33 |
Fig. 4 : Sédimentation d'une bille dans un fluide visco-élastique (We=1, r=0.5) en régime non inertiel (Re=0). L'interaction fluide-structure est gérée par pénalisation du tenseur de déformation et du tenseur de contrainte élastique. Les films (au format .gif animé) associés aux images ci-dessous représentent l'évolution en temps des inconnues du système dans le cas newtonien (Re=0, We=0, r=0) et dans le cas visco-élastique (Re=0, We=1, r=0.5).
 |  |  |
| a) Champ de vitesse | b) Cisaillement élastique | c) Élongation élastique |
Fig. 5 : Sédimentation de deux billes dans un fluide newtonien ou visco-élastique : phénomènes de drafting, kissing and tumbling dans le cas newtonien et de chaining dans le cas visco-élastique.
 |
 |
|
| a) Cas newtonien | b) Cas visco-élastique |
Publications :
L. Chupin, S. Martin : Stationary Oldroyd model with diffusive stress: mathematical analysis and vanishing diffusion process, Journal of Non-Newtonian Fluid Mechanics, 218, 27–39 (2015).
A. Decoene, S. Martin, B. Maury : Direct simulation of rigid particles in a viscoelastic fluid, Journal of Non-Newtonian Fluid Mechanics, 260, 1–25 (2018).