Simulation de fluides sur des surfaces

De Ensiwiki
Aller à : navigation, rechercher


Simulation de fluides sur des surfaces
Projet Projet de spécialité 2A
Sous-projet Image
Étudiants Mathieu BAILET (MMIS-MCS)

Bastien FRISULLI (MMIS-MCS)

Joan MARTÍN HERNÁNDEZ (MMIS-IRVM)

Victor MICHEL (MMIS-IRVM)

Promo 2011
Tuteur Stefanie HAHMANN stefanie.hahmann@grenoble-inp.fr

Introduction

La problématique initiale était de pouvoir "simuler des fluides sur des surfaces à topologie quelconque". Ce sujet nous semblait intéressant puisqu'il allié à la fois des problèmes de résolution d'équation différentielle et d'imagerie 3D.

Objectifs attendus

Les objectifs n'étaient pas explicites puisqu'une certaine liberté nous été proposée et nous pouvions ainsi explorer plusieurs aspects du problème. Cependant il était nécessaire d'arriver à un résultat concret de visualisation des fluides sur des formes simples telles que des boules ou des tores.

Résolution 2D

La première étape de notre projet était la simulation de fluides sur une surface plane. On considère ainsi un carré de taille 1 de côté comme surface de base de simulation.

Pour cela il a été nécessaire de résoudre les équations qui modélisent l'évolution des fluides et ensuite de simuler de manière intuitive. Nous allons donc séparer la première étape en deux sous-étapes: résolution des équations et visualisation des solutions.

L'équation de fluides

L'équation qui décrit l'évolution du système dynamique où les fluides se déplacent est donné par les équations de Navier-Stokes :

\left\{
    \begin{array}{ll}
  \frac{\partial u}{\partial t} = -(u \cdot \nabla )u + \frac{1}{\rho}\nabla p + \nu \nabla^2 u +f \\ \nabla \cdot u = 0 
    \end{array}
\right.


Ces équations peuvent être séparer en quatre parties représentant les caractéristiques d'un mouvement de fluide :

Force externe : Dans l'espace du fluide, plusieurs forces externes s'appliquent, ce qui influe sur le déplacement du fluide.

  f  : la variation du valeur du fluide par rapport au temps dépend directement de la force externe.

Advection : Lorsqu'un fluide se déplace, il suit une certaine trajectoire dictée par la partie "advection" de l'équation.

 -(u \cdot \nabla )u 

Diffusion : Un fluide tend toujours à se propager dans le fluide externe afin de s'homogénéiser. La diffusion est dictée par cette partie de l'équation.

 \nu \nabla^2 u 


Projection : Il est nécessaire d'assurer la conservation de masse de notre fluide. L'équation suivante doit donc être conservé au fur et à mesure.

 \nabla \cdot u = 0  : La variation de notre fluide ne doit pas dépendre de sa position, mais seulement de l'évolution par rapport au temps.

Résolution de l'équation

L'équation est résolue depuis en partant de l'état initial u_0(x) = u(x,0) en parcourant le temps par pas de \Delta t. Lorsqu'on a la soluion au temps t, on calcule la solution au temps t+\Delta t en quatre étapes en partant de w_0(x) = u(x,t).

Addition de force

A chaque itération de notre système dynamique nous devons ajouter la force externe (poussée d'Archimède, gravité ...)

  w_1(x) = w_0(x) + \Delta t \times f(x,t) 

Advection

Afin de résoudre le système  \frac{\partial u}{\partial t} = - u\cdot \nabla u \quad nous utilisons la méthode des caractéristiques (voir ci-après) avant une approximation numérique.

Méthode des caractéristiques

Soit  v(x,t) un champ vectoriel,  a(x,t) un champ scalaire et  a_0(x) = a(x,0) .

On considère l'équation

 \frac{\partial a(x,t)}{\partial t} = -v(x)\cdot\nabla a(x,t) . 

et on définit  p(x,t) comme l'ensemble de courbes solution du champ  v(x,t)  :

 \left \{ \begin{array}{ll} 
\frac{d}{dt}P(x_0,t) = v(P(x_0,t),t)\\
P(x_0,0)=0 \end{array} \right . 

Soit

 \bar{a}(x_0,t) = a(P(x_0,t),t)  la valeur du champ  a(x,t)  au cours de la caractéristique
 P(x_0,t) .

Alors :

 \frac{d\bar{a}}{dt} = \frac{\partial a}{\partial t} + v \cdot \nabla a = 0

La valeur du champ scalaire ne change donc pas le long de chaque caractéristique.

En particulier :

 \bar{a}(x_0,t) = \bar{a}(x_0,0) = a(x_0,0) = a_0(x_0)
Application de la méthode

D'après les calculs précédents :

  w_2 (x) = w_1(p(x,-\Delta t)) 

On peut alors résoudre le problème numériquement :

Pour chaque point (s,t) il s'agit de revenir en arrière en passant par la caractéristique correspondante à son point afin de récupérer la valeur de la densité et des champs de vitesses.

Or la discrétisation, nous oblige à traiter les cas aux bords où nous devons aller récupérer des valeurs en dehors la limite de notre maillage. Dans ce cas précis, nous faisons une restriction afin de rester à l'intérieur du maillage.

Carac.png

Diffusion

La diffusion mesure la capacité du fluide à se répandre dans son milieu sans apport de forces extérieures. Dans l'équation de Navier-Stokes, le terme "diffusif" de l'équation est :

  \frac{\partial u}{\partial t} = \nu \nabla^2 u 

Pour des raisons de stabilité des méthodes de résolution, il a été montré qu'il est préférable d'utiliser une méthode implicite :

  \left( I - \nu \Delta t \nabla^2 \right) w_3(x) = w_2(x) 

La discrétisation de l'opérateur linéaire :  I - \nu \Delta t \nabla^2 donne une matrice creuse, avec 5 diagonales. En effet, le Laplacien se discrétise en la matrice  n^2 \times n^2 comportant 5 diagonales, avec des 4 sur la diagonale principale, et des -1 sur les autres.

Cette particularité permet d'écrire un solveur très efficace, en l'occurence un solveur de type Gauss Seidel (donc par relaxation), adapté à notre système creux, car nous n'avons pas à avoir de structure de données pour stocker la matrice de l'opérateur Laplacien.

Projection

La projection consistait à utiliser la décomposition de Helmholtz-Hodge du vecteur vitesse, comme si dessous:

w = u + \nabla q 

où la divergence de u est nulle et q est un champ scalaire.


En prenant le produit scalaire par \nabla, on obtient :

\nabla \cdot w_3 = \nabla ^2 q

, que l'on résout avec un solveur spécifique aux équations de Poisson.

De même que pour la diffusion, nous avons choisi un solveur de Gauss-Seidel pour les mêmes raisons (spécificités de la matrice, vitesse de convergence, ...).

Une fois cette équation résolue, nous obtenons q et ainsi nous en déduisons la projection :

 Pw_3 = w_4 = w_3 - \nabla q.
Method.png

Conditions aux bords périodiques

Si on considère des conditions aux bords périodiques, c'est à dire en ayant un domaine de topologie torique, on peut utiliser la transformée de Fourier pour la partie diffusion et projection. La méthode est simple et très élégante car les étapes de projection et de diffusion se réduisent à des opérations simple dans l'espace de Fourier. En effet, l'opérateur \nabla est équivalent à une multiplication par ik

La diffusion devient :

1 + \nu \Delta t k^2

La projection devient :

\hat w(k) - \frac{1}{k^2} (k .\cdot \hat w(k)) k

Dans, ce cas, nous n'avons besoin pour notre solveur que d'un schéma numérique pour l'advection et de la fft.

Toutefois, la complexité de l'algorithme devient O(n^2 log n^2) avec n la taille de la discrétisation.

Nous avons utilisée la librairie fftw pour calculer la fft.

Visualisation de la simulation

Visualisation de la densité

Voici des captures d'écrans de notre programme OpenGL réalisant la simulation. Dans cet exemple, Nous avons placé une source "infinie" de particules au centre du domaine. Comme l'on peut le voir, les conditions aux bords empêchent le fluide de sortir du domaine.

0f2d.png
1f2d.png
2f2d.png
3f2d.png
4f2d.png
5f2d.jpg
6f2d.png
7f2d.png

Visualisation du champ de vitesse avec la méthode de convolution de ligne intégrale (LIC)

La technique du LIC permet de visualiser un champ vectoriel grâce à une image quelconque (mais la visualisation est meilleure avec du bruit blanc) dont on trace le chemin dans l'image en fonction du champ vectoriel. Une fois le chemin calculé, on lui applique une convolution afin d'uniformiser la couleur et de voir la ligne de courant. Le noyau de convolution choisi est un vecteur dont les valeur sont  \frac{1}{n} avec  n la valeur de discrétisation.

Bruit Blanc
LIC du bruit blanc en fonction du champ vectoriel

Visualisation en couleur avec conditions aux bords périodiques

Fluid1.jpg
Fluid2.jpg

Version GPU

L'algorithme de simulation 2D étant implicitement parallèle, nous avons décidé de le porter en CUDA.

Pour ce faire, nous avons choisi le découpage suivant :

  • 1 thread = 16 cellules verticales
  • 1 bloc = 64 threads horizontaux x 4 threads verticaux = 4096 cellules

Le découpage de la grille étant fait automatiquement selon la taille du problème.

GPU bloc.png

Un point critique a été la résolution des systèmes linéaires. Notre version CPU utilisait l'algorithme de Gauss Seidel (adaptée spécifiquement à notre problème creux), mais cet algorithme est intrinsèquement séquentiel. Nous avons donc utilisé une version parallèle de Gauss Seidel, appelée "Red Black Gauss-Seidel", afin d'augmenter les performances. Pour aller plus loin, il existait aussi des méthodes beaucoup plus poussées (multi grilles), mais nous ne sommes pas allés aussi loin.

Pour le rendu OpenGL, nous avons utilisé un VBO partagé entre CUDA et OpenGL, nous ne sommes donc pas limités par la bande passante du PCI Express.

Les tests ont été faits sur une NVidia Geforce GTX 275. Le gain mesuré par rapport à une version séquentielle CPU sur Intel Core i7 920 a été de l'ordre de 15 quant à la taille de la discrétisation, pour un rendu avec framerate acceptable.

Faute de temps, la parallélisation a été relativement naïve et pourrait être améliorée : par exemple, pour certains opérateurs linéaires, chaque thread doit accéder à ses 4 voisins, chaque valeur est donc lue 4 fois dans la mémoire globale, ce qui pourrait être ramené à 1 en utilisant un tableau intermédiaire en mémoire partagée.

Résolution sur des surfaces

La deuxième étape de notre projet s'agit de simuler des fluides sur surfaces tridimensionnelles de topologie arbitraire (on suppose toujours que les surfaces sont fermées).

Choix du type de surface

Surfaces de Catmull Clark

Les surfaces de Catmull Clark sont des surfaces de subdivisions très utilisées dans les modeleur 3D pour l'efficacité de leur algorithme de construction. En effet avec peut de patchs et un nombre faible d'itérations, on peut avoir des surfaces assez lisses. Voici des captures d'écran du programme que nous avons réalisé :

Itération 1
Itération 2
Itération 3
Itération 4

Les surfaces de Catmull-Clark tendent vers des surfaces BSpline bicubiques dont avec un raccordement G1 ce qui convient pour notre problème. Toutefois, l'évaluation pour un point quelconque d'une Surface de Catmull Clark est fastidieuse.

Surfaces de Bezier

Les surfaces de Bézier sont une méthode de définition d'une surface grâce aux courbes de Bézier, avantageuses pour définir une courbe par la donnée de points de contrôle.

Étant donnée une matrice [M] de points de l'espace A_{i,j}, la surface de Bézier correspondante est l'ensemble des points M généré par les valeurs comprises entre 0 et 1 des variables u et v du polynôme :

\overrightarrow {O M} = 	\sum_{i=o}^n 	\sum_{j=o}^m C_n^i v^i (1-v)^{(n-i)} C_m^j u^j (1-u)^{(m-j)} \overrightarrow {O A_{i j}}

On peut créer des surfaces de Bezier en raccordant des patchs de Bezier de façon à avoir une surface G1. De plus, l'évaluation d'un point sur un patch de Bezier tout comme ses dérivées partielles sont très rapides grâce à l'algorithme de De Casteljau.

Placement des points de contrôle pour avoir un raccordement G1
Raccordement G1 de deux patchs de Bezier

Les problèmes principaux

Problème de comparaison de vecteurs dans différents plans tangents

Simuler un fluide sur une surface quelconque pose différents problèmes qui sont :

1. Les conditions aux bords : Lorsque notre surface est composée de différents patchs, il est nécessaire d'avoir un raccord entre ces différentes parties afin d'avoir une continuité du fluide lorsque celui ci les traverse.

2. L'usage de la métrique : nous devons simuler sur des surfaces de forme quelconque, or il est évident que le fluide différemment sur une surface plane et une surface déformée. La prise en compte de la métrique est donc nécessaire afin d'avoir des effets réalistes.

3. La performance de l'algorithme : Le nombre de points avec lesquels nous travaillons étant important, puisque dans certains cas nous devons traiter avec des centaines de patchs qui peuvent avoir une discrétisation très fine. Il est donc absolument nécessaire de faire un code performant et optimisé. De plus, les algorithmes de résolutions doivent être adaptés à nos nouveaux systèmes d'équations.

Modélisation de la surface

Raccord des patchs

Les surfaces sur lesquelles nous avons choisi de faire la simulation sont les surfaces de Bézier, obtenues par recollement G1 de patchs de Bézier. L'évaluation de la surface et de ses dérivées partielles étant très simple (algorithme de Casteljau), c'était un avantage par rapport aux surfaces de Catmull Clark, ou l'évaluation est plus complexe.

A chaque patch de Bézier est associé une discrétisation pour la résolution de l'équation. Chaque patch dispose alors d'un système de coordonnées propres. Pour des raisons topologiques, il n'est pas possible d'étendre le système de coordonnées d'un patch pour inclure tous les patchs de la surface (plusieurs cellules pourraient partager les mêmes coordonnées).

Nous avons donc gardé un système de coordonnées locales pour chaque patch, et programmé des fonctions de transition pour transporter le système de coordonnées d'un patch vers un autre, notamment pour accéder aux cellules voisines.

Résolution de l'équation de fluides

L'objectif est ici de transporter la résolution de l'équation de Navier Stokes sur une surface de Bézier.

La transition n'est pas évidente : projeter directement la solution de l'équation 2D sur la surface peut créer des effets indésirables selon la surface. En effet, l'image de la grille de discrétisation peut être fortement distordue sur la surface comme nous le constatons sur le patch de Bézier suivant.

Grille distordue

Pour résoudre ce problème, nous devons prendre en compte la métrique de la surface. Pour cela, nous utilisons la première forme fondamentale, qui est l'image du produit scalaire usuel du plan tangent à la surface dans l'espace des paramètres (par le difféomorphisme de Bézier inverse). Cette image est aussi un produit scalaire, car les surfaces de Bézier considérées sont suffisamment régulières par le choix des points de contrôle.

Nous obtenons ainsi une norme dans l'espace des paramètres, ce qui nous permet de discrétiser les opérateurs linéaires dont nous avons besoin (gradient, divergence, laplacien et rotationnel). Ainsi, nous pouvons traduire le solveur 2D vers un solveur 3D. Les systèmes considérés sont toujours creux, mais plus complexes qu'en 2D : 9 diagonales, avec des valeurs non répétées le long de celles-ci. Nous avons implémenté une version creuse du gradient conjugué pour résoudre le problème associé.

Un autre point important est la communication entre patchs : on ne peut pas se permettre de résoudre l'équation séparément sur tous les patchs, puisque le fluide doit s'écouler surt toute la surface. La grille de discrétisation du problème est donc plus large que la grille de discrétisation de la surface de Bézier. Ainsi, nous ajoutons une frontière de largeur 1 autour de chaque grille. Ces cellules "fantôme" sont en fait la valeur des patchs voisons aux frontières.

Prise en compte des patchs voisins

Addition de force

Projection de la force sur le plan tangent

L'addition de la force est plus complexe que dans la version 2D du problème :

Nous travaillons sur une surface de l'espace où le fluide est en mouvement et nous souhaitons lui ajouter une force externe. Pour cela nous ne pouvons pas ajouter directement la force aux champs de vitesse de la même manière que pour la 2D, puisque les résultats obtenus sur l'espace des paramètres seront différents sur la surface selon sa forme.

La bonne façon de faire est de projeter la force sur le plan tangent en chaque point, le vecteur résultant étant le vecteur avec lequel on va mettre à jour le champ de vitesse.

Nous supposons que nous pouvons écrire en tout point de l'espace la force sous la forme :

  f(x,y,z) = (f_1(x,y,z),f_2(x,y,z), f_3(x,y,z)) \quad 

Afin de calculer la projection de ce vecteur sur le plan tangent en un certain point P(x,y,z) de la surface, nous devons faire le produit scalaire entre les vecteurs  \frac{\partial}{\partial u} et  \frac{\partial}{\partial v} selon l'endroit où le point est définit :

  p_0 = \phi(u_0,v_0) \quad 
 
  \Pi_{Tp_{0}S}(f(p_0)) = \frac{<f(p_0),\frac{\partial}{\partial u} \phi(p_0)>}{|\frac{\partial}{\partial u} \phi(p_0)|^2} \frac{\partial}{\partial u} \phi(p_0) + \frac{<f(p_0),\frac{\partial}{\partial v} \phi(p_0)>}{|\frac{\partial}{\partial v} \phi(p_0)|^2} \frac{\partial}{\partial v} \phi(p_0)  \quad 

Nous allons ainsi modifier les champs de vitesse avec les coefficients de chaque composante du vecteur précédent normalisé.

 for all patch of the surface
   for all (s,t) in [0,1]x[0,1]
     u(s,t) = (u(s,t) + <force(S(s,t)), dS(s,t,u)>)/dS(s,t,u).norm();
     v(s,t) = (v(s,t) + <force(S(s,t)), dS(s,t,v)>)/dS(s,t,v).norm();
   end;
 end;

Nous appliquons cet algorithme sur chaque patch de la surface.

Advection

La prise en compte de l'orientation relative des patchs est nécessaire pour le codage de l'advection, de même que les conditions aux bords (interaction, continuité entre les bords de chaque patch).

Nous nous sommes alors inspirés de l'advection en 2D pour ce nouvel algorithme en appliquant des modifications simples.

En 2D, lorsque dépassions les bords, nous introduisions une restriction. Ici, il est nécessaire de remarquer que si lors du calcul, le point à aller chercher est en dehors du patch sur lequel nous travaillons, nous calculons alors les nouvelles coordonnées de ce point par rapport au patch auquel il appartient et nous récupérons la valeur afin d'affecter le point initial correctement. De plus, nous devons appliquer une rotation à ce vecteur pour qu'il correspondent aux coordonnées du patch

Rotations.jpg

Diffusion

Comme en 2D, nous devons résoudre l'équation :

  \left( I - \nu \Delta t \nabla^2 \right) u_1 = u_0 

Mais l'expression du Laplacien selon la métrique de la surface est beaucoup plus complexe. Nous ne donnons pas ici les formules utilisées, nous renvoyons à l'article de Jos Stam pour les détails.

Projection

De même que pour la diffusion, l'équation afin de résoudre la projection est totalement différente lorsque la métrique est à gérer. Nous utilisons une nouvelle fois la décomposition de Helmholtz-Hodge :

w = u + \nabla q

\nabla . u = 0

Nous devons résoudre l'équation suivante :

\nabla ^2q = \nabla . w 

Ce qui donne lorsque nous introduisons la métrique :

\frac{\partial}{\partial t}(\sqrt{g}g^{i,j}\frac{\partial q}{\partial x^j} ) = \frac{\partial}{\partial x^i}(\sqrt{g}u_3^i)

Nous avons ainsi :

(\nabla . u)_{i,j} = \frac{1}{(\sqrt{g})_{i,j}h}(D^1_{i,j}+D^2_{i,j})

avec

D^1_{i,j} = (\sqrt{g})_{i+\frac{1}{2},j}(u^1)_{i+\frac{1}{2},j}- (\sqrt{g})_{i-\frac{1}{2},j}(u^1)_{i-\frac{1}{2},j}

D^2_{i,j} = (\sqrt{g})_{i,j+\frac{1}{2}}(u^1)_{i,j+\frac{1}{2}}- (\sqrt{g})_{i,j-\frac{1}{2}}(u^1)_{i,j-\frac{1}{2}}

De manière identique à la partie 2D, nous avons utilisé le même solveur que pour la diffusion, à savoir le gradient conjugué. Enfin après avoir trouvé q, nous calculons son gradient et nous en déduisons la projection comme ci-dessous :

u_4^k = u_3^k - g ^{k,j} \frac{\partial q}{\partial x^j}

et la convention suivante ("à la Einstein")

g^{i,j} \frac {\partial}{\partial x^j} = g^{i,1} \frac {\partial}{\partial x^1}+g^{i,2} \frac {\partial}{\partial x^2}

Visualisation de la simulation

Voici des captures d'écran de notre programme réalisant la simulation de fluides sur des surfaces de Bézier - ici, une boule (24 patchs de Bézier) et un tore (64 patchs)

Simulation de fluide sur une "sphère" formée de patchs de Bézier
Simulation de fluide sur un "tore" formé de patchs de Bézier

Une vidéo !

Voici une vidéo temps réel de notre simulation de fluides, ici, sur un tore

http://vic.michel.perso.sfr.fr/surface.mp4 (~3 Mo, codec H.264)

Conclusions

Comme nous l'avons expliqué lors de l'introduction, le sujet est vaste, et nombreuses sont les pistes à explorer ou à continuer de développer :

Continuité entre patchs : Dans notre programme 3D, la continuité entre nos patchs subissaient à certains moments quelques bugs. Un travail à faire serait donc de trouver ces bugs et de les rectifier.

Utilisation de la métrique : La résolution complète en utilisant la métrique de la surface est à finir

Ajout d'autres forces : Pour l'instant seule la force de gravité est prise en compte, l'ajout d'autres forces comme celle de Coriolis serait bénéfique dans la visualisation finale.

Amélioration de l'interactivité et rendu OpenGL : Certaines améliorations au niveau du résultat et de l'interactivité avec l'utilisateur sont à faire, de même que certains codes OpenGL.

En termes personnels, ce projet fût très enrichissant pour nous tous, le thème touchant directement nos filières et à nos intérêts.

Nous sommes aussi très satisfaits de nos résultats de simulation, au niveau des résultats en 2D (notamment notre simulation avec conditions de Fourier), comme des résultats en 3D sur des patchs et surfaces de Bézier.

Bibliographie

- Jos Stam, "Stable Fluids", In SIGGRAPH 99 Conference Proceedings, Annual Conference Series, August 1999, 121-128

- Jos Stam, "Flows on Surfaces of Arbitrary Topology",ACM Transactions On Graphics (TOG), Volume 22, Issue 3 (July 2003) : Proceedings of SIGGRAPH 2003, 724-731

- Brian Cabral, "Imaging Vector Fields Using Line Integral Convolution", Lawrence Livermore National Laboratory

- Applet de simulation de fluides

- Instants Fun

- Wololo Blues