Projet image 2012 : Visualisation de fluides

De Ensiwiki
Aller à : navigation, rechercher

Tux.png  Linux  CDROM.png  Informatique  Mathematiques.png  Mathématiques 

Project schedule.png
Titre du projet Projet image 2012 : Visualisation de fluides
Cadre Projets de spécialité
Page principale Projet image 2012 : Visualisation de fluides

Encadrants Stefanie Hahmann


Introduction

Exemple de suminagashi

Cette page a pour objectif de présenter le travail réalisé dans le cadre de notre projet de spécialité 2a. La version originale du sujet est disponible sur le site officiel du sujet.

Ce projet consiste à modéliser et visualiser le comportement d'un solide plongé dans un fluide. Pour illustrer ceci, on peut prendre l'exemple d'un art traditionnel japonais appelé suminagashi. Il consiste à faire flotter de l'encre sur de l'eau et souffler sur le bassin de façon à obtenir un motif aléatoire.
L'encre est ici le solide, et l'eau le fluide.

Équipe

Notre équipe pour ce projet était constituée de deux personnes :

Nous avons également été aidés par notre tuteur :

Objectifs et contexte

La présentation du sujet était très axée sur le domaine artistique. En effet, les images produites sont esthétiques et peuvent ainsi être utilisées en tant que motif ou texture (par exemple en reliure de livre). Cependant, il peut aussi être observé du point de vue simulation. Nous visualisons les effets d'un fluide sur un solide suite à l'application de forces particulières (appliquée par l'utilisateur). Notre application pourrait donc être la base d'un simulateur d'écoulements d'un fluide, utile dans des domaines tels que l'aérodynamique. Elle serait alors un outil de simulation permettant de remplacer des tests coûteux en soufflerie.

Les objectifs de notre projet ont ainsi étés les suivants :

  • S'appuyer sur des articles scientifiques pour réaliser une application permettant une visualisation de fluides réaliste
  • Permettre à l'utilisateur d'interagir avec le simulateur, i.e appliquer des forces, rajouter des éléments, etc
  • Améliorer les méthodes du simulateur pour effectuer les calculs plus rapidement
  • Réfléchir à des méthodes permettant d'améliorer la qualité de l'image affichée, et les essayer

Travail réalisé et livrables

Nous avons ainsi réalisé une application de visualisation de fluides en deux dimension, en implémentant les techniques présentées dans les articles listés dans la bibliographie. Nous avons également ajouté certaines améliorations proposées dans ces articles, que ce soit au niveau du simulateur, ou du rendu visuel.
Notre application est fournie avec un manuel utilisateur permettant de comprendre toutes les fonctionnalités de notre application. Nous avons également réalisé cette page afin de présenter notre travail, les résultats que nous avons obtenus, et dresser un bilan de notre projet.

Modélisation

Dans le cadre de notre sujet, un fluide est caractérisé par un champs de vecteurs vitesse. Le comportement des fluides est aujourd'hui bien connu et peut être décrit à partir d'équation physiques. Nous allons utiliser l'équation principale de la mécanique des fluides pour modéliser le comportement du champs de vecteurs : l'équation de Navier-Stokes.

Équation de Navier Stokes

L'équation de Navier Stokes de base est la suivante :

 
\left\{
\begin{array}{ll}
\nabla \cdot u = 0 & \text{(0)}\\
\frac{\partial u}{\partial t} = (\underbrace{-(u \cdot \nabla) u}_{\text{advection}} - \underbrace{\frac{1}{\rho} \nabla p}_{\text{pression}} + \underbrace{\nu \Delta u}_{\text{diffusion}} + \underbrace{f}_{\text{addition de force}}) & \text{(1)}
\end{array}
\right.

Avec :

  • u : vitesse du fluide
  • p : pression
  • \rho : densité
  • \nu : viscosité

Cependant ce n'est pas l'équation principale que nous utiliserons. Nous allons réduire celle-ci en supprimant le terme relatif à la pression. Pour cela, nous allons projeter l'équation (2) en utilisant la décomposition de Helmholtz-Hodge. Tout vecteur se décompose de manière unique en une somme d'un gradient et d'un champ scalaire. On peut ainsi définir un opérateur de projection P projetant un vecteur sur sa composante de divergence nulle.L'équation devient alors :

  
\left\{
\begin{array}{ll}
\nabla \cdot u = 0 & \text{(0)}\\
\frac{\partial u}{\partial t} = \underbrace{P}_{\text{projection}}(\underbrace{-(u \cdot\nabla) u}_{\text{advection}} + \underbrace{\nu \Delta u}_{\text{diffusion}} + \underbrace{f}_{\text{addition de force}}) & \text{(3)} 
\end{array}
\right.

C'est à partir de cette équation que nous développerons notre solveur.

Résolution de l'équation de Navier-Stokes

La résolution d'une telle équation différentielle se fait par étapes, en ajoutant au fur et à mesure les différents membres de l'équation (3). On commence la résolution de l'équation avec la solution calculée précédemment :

 w_{0}(x) = u(x,t)

La procédure se décompose en quatre étapes :

  w_0 \xrightarrow{\text{addition de force}} w_1 \xrightarrow{\text{advection}} w_2 \xrightarrow{\text{diffusion}} w_3 \xrightarrow{\text{projection}} w_4

La solution au temps t+ \Delta t est alors w_{4}(x).

Étape d'ajout de force

L'équation se limite donc à:

\frac{\partial w_{1}}{\partial t} = f 

En considérant que la force est constante sur la période du pas de temps \Delta t, on obtient la solution w_{1}(x) de cette équation :

 w_{1}(x) = w_{0}(x)+\Delta t f(x,t)

Étape d'advection

Étant donné que w_{1} satisfait l'équation précédente, l'équation à résoudre est la suivante :

 \frac{\partial w_2}{\partial t}= -(w_1 \cdot\nabla)w_2

Pour résoudre celle-ci, nous utilisons une méthode présentée l'article [2] appelée méthode des caractéristiques. Il faut imaginer le fluide comme un ensemble de particules déplacées à différentes vitesses. A chaque pas de temps, ces particules bougent et se retrouvent dans une nouvelle position. Cette méthode consiste à trouver la particule à la date t-\Delta t qui se retrouve à la position étudiée à la date t. Sa vitesse est alors celle recherchée. Ceci nécessite deux étapes :

  • rechercher la position de la particule
  • interpoler la vitesse de la particule à partir de sa position

Le calcul de la vitesse est nécessaire car le fluide est discrétisé, ce qui implique une interpolation pour tout calcul de vitesse. Il existe plusieurs façons d'interpoler les vitesse. Nous rediscuterons des choix effectués dans une autre section.

Étape de diffusion

L'équation est ici :

 \frac{\partial w_3}{\partial t} = \nu \Delta w_3

Cette équation assez classique, peut être résolue de différente manières. Une façon simple de faire, que nous avons adopté, est de discrétiser cette équation. On obtient alors l'équation suivante :

 {w_3}_{i, j} = {w_2}_{i, j} + \nu \Delta t ({w_2}_{i-1, j} + {w_2}_{i+1, j} + {w_2}_{i, j-1} + {w_2}_{i, j+1} - 4  {w_2}_{i, j})/h^2 

Avec h la taille d'un élément de notre grille de discrétisation.
Cette équation est un système linéaire, et peut donc être résolu par une méthode classique telle que Gauss-Seidel.

Étape de projection

Nous cherchons à présent la projection de w_{3} permettant d'avoir un résultat de divergence nulle. L'équation est un problème de Poisson :

w_{4} = P w_3 = w_{3} - \nabla q

Avec q un champ scalaire. L'inconnue est ce champ scalaire. Pour le calculer, on calcule la divergence de w_{3} ce qui nous donne l'équation :

 \nabla \cdot w_{3} = \nabla^2 q = \Delta q 

Une fois encore, on discrétise cette équation pour obtenir :

 \frac{{w_{3x}}_{i+1, j}- {w_{3x}}_{i-1j}}{2h} + \frac{{w_{3y}}_{i, j+1}- {w_{3y}}_{i, j-1}}{2h} = \frac{q_{i-1, j}+ q_{i+1, j}+ q_{i, j-1}+ q_{i, j+1} - 4q_{i, j}}{h^2}

De nouveau, il s'agit d'un système linéaire. On peut utiliser le même solveur que celui de l'étape de diffusion.

Création de Fluid solver

Maintenant que nous disposons d'un solveur permettant de calculer l'évolution du champ de vecteurs vitesse, il faut l'insérer dans notre application et l'utiliser pour visualiser les fluides. Dans le contexte de notre projet, un solide est un contour fermé : c'est à dire que ses frontières peuvent être dessinées sans lever le crayon.
En ce qui concerne les technologies utilisées, nous avons optés pour QGLViewer, interface qt qui utilise OpenGL. Nous allons présenter dans cette section comment nous avons utilisé le solveur construit précédemment pour déformer ces solides. Nous évoquerons également de quelle manière nous affichons le résultat de notre simulation.

Méthode

Comme dit précédemment, les solides qui nous intéressent sont des "contours fermés". Ils sont caractérisés par ce contour. Nous allons donc calculer les nouvelles positions du contour, puis colorer l'intérieur de celui ci.

représentation d'un contour fermé

Ce contour est modélisé par un ensemble de points reliés. Le principe est donc de déplacer ces points en calculant leur nouvelle position à partir du champ de vecteurs vitesse du solveur. Lorsque ces nouvelles positions sont obtenues, on peut chercher à améliorer le rendu (voir section Améliorations apportées).

Déplacement du contour

Notre modèle discrétise l'espace et le temps. En effet, on calcule la nouvelle position des points à intervalle de temps fixé, et une grille discrète est utilisée pour stocker le champ de vecteurs vitesse.
Pour calculer la position d'un point après modification du champ de vecteurs vitesse, on intègre sa position dans le champ par la méthode de Range-Kutta d'ordre 4. La nouvelle position de chaque point dépend de son ancienne position, de la vitesse de ce point, et du pas de temps :

  u^{t+\Delta t} = \phi(\Delta, u^{t}, v^{t})

Avec :

  • u^{t} position au temps t
  • v^{t} vitesse au temps t
  • \Delta t pas de temps utilisé
  • \phi la fonction de calcul du schéma de Range-Kutta

Remplissage du contour

Nous avons utilisé une méthode classique de remplissage utilisée pour les polygones concaves. Cette méthode s'appelle méthode de Stencil. Elle consiste à utiliser un buffer pour stocker les pixels devant être colorés et ceux devant être masqués. Opengl met à disposition un stencil buffer que nous avons utilisé.

Améliorations apportées

Rendu avec un nombre de points fixe

Notre avons cherché à améliorer notre application pour rendre plus rapide les calculs et obtenir un meilleur rendu. Ces améliorations sont particulièrement importantes dans le cadre du rendu. En effet, sans celle-ci, le contour est définit avec un nombre constant de points. Lorsque la forme devient complexe, la distance entre les points augmente et entraîne des collisions (illustré à gauche ci-dessous), ainsi qu'une absence de courbe (figure de droite).

Les méthodes utilisées sont présentées dans les articles [3][4].

Courbes de Bézier

Une première idée est de remplacer les segments entre les points du contour par des courbes. Reste à choisir le type de courbe. Nous avons opté pour une des courbes les plus utilisées dans le domaine de l'infographie : les courbes de Bézier. Celles-ci sont définies à partir de quatre points de contrôle.

points de contrôle

Avec : v_{i} = p_{i+1}-p_{i-1}
Pour calculer la courbe entre p_{i-1} et p_{i} on utilise les points de contrôle (p_{i-1}, p_{i-1} + v_{i-1}/8, p_{i} - v_{i}/8, p_{i}) Le nombre de points définissant le contour est donc toujours stable, mais ceux-ci permettent de définir de nouveaux point pour l'affichage.

Ajout de points

Une seconde amélioration consiste à ajouter des éléments à la liste définissant le contour dès que la distance entre deux points dépasse une certaine distance. Une fusion de points (retirer de la liste) est également prévue lorsque les points sont trop proches. La notion de "trop proche" étant définie à partir de la distance précédente.
La difficulté consiste à choisir cette distance. Dans cette première amélioration, celle-ci est constante. Cela mène rapidement à un nombre de points très important ce qui ralenti l'application.

Rendu avec ajout de points et fusion

Courbure

L'ajout de points donne de bons résultats, mais devient rapidement lourd lorsque les formes deviennent complexes. L'idée est alors de faire varier la distance maximum et minimale entre les points en fonction de la figure. Si celle ci est "droite", peu de points sont nécessaires. Si au contraire il y a une courbe importante, il faut définir plus de points.
La distance d_{i} est donc calculée de façon à ce qu'elle dépende d'un coefficient de courbure :

 d_{i} = d_{max}(\epsilon){c_{i}}^{courbure} + \epsilon

Avec \epsilon une distance minimale et d_{i} =
d_{max}(\epsilon) une distance maximale dépendant de \epsilon.
La formule utilisée pour définir {c_{i}}^{courbure} est la suivante :

 {c_{i}}^{courbure} = exp(-a \times |k|)

Avec k la courbure définie par la formule approchée suivante (où l'on considère que les points du contours sont les points également espacés dans le temps d'une courbe paramétrée):

 k = \frac{\ddot{y}\dot{x} - \dot{y}\ddot{x}}{(\dot{x}^2 + \dot{y}^2)^{\frac{3}{2}}}

En utilisant cette distance variable, nous obtenons une répartition des points plus importante dans les courbes.

Rendu avec distance dépendant de la courbure

Turbulences

L'article [3] propose de prendre en compte les "turbulences" du fluide. L'idée étant que lorsque les vitesses sont élevées et fortement variables dans une région, il y aura besoin de plus de points pour définir le contour. Nous avons ajouté cette méthode en calculant un nouveau coefficient :

 {c_{i}}^{turbulence} = exp(-b \times |\nabla \times u(v_{i})|)

Avec v_{i} la position du vertex i et u(v_{i}) sa vitesse.

Rendu avec distance dépendant des turbulences

Combinaison

On peux combiner les deux coefficients précédents pour définir la distance selon la formule suivante :

 d_{i} = d_{max}{c_{i}}^{courbure}{c_{i}}^{turbulence} + \epsilon

On a alors le rendu suivant :

Rendu avec distance dépendant des deux coefficients

On remarque que l'on a une répartition avec beaucoup de points aux bons endroits mais on se retrouve avec des agrégats. En effet, les points sont mal répartis le long du contour.

Diffusion

Pour homogénéiser la répartition des points, on pondère la distance entre deux points par les autres distances des paires de points proches. On utilise ici une gaussienne qui nous permet de diffuser plus ou moins les distances par un plus ou moins grand nombre de distance voisine.

Rendu avec diffusion

Au final, on obtient une "belle" répartition des points sur le contours en jouant sur les différents paramètres du modèle. En comparant les images avec le raffinement et sans, on observe l'intérêt de la méthode: avoir un contour plus précis, et qui le reste quelque soit la perturbation apportée.

On peut noter qu'il existe encre d'autres méthodes pour améliorer le contour, qui consistent majoritairement en des tests de distances avec les contours opposés pour, par exemple, diviser le contour en deux. Cependant ces méthodes nous ont semblé trop lourdes à implémenter et pas dans le cadre du sujet de contours fermés.

Performances et choix

Nous avons désignés cette application dans l'optique d'être facilement utilisable par l'utilisateur et surtout fluide, ce qui a pour contrepartie d'utiliser des approximations. Les 2 points principaux pour augmenter les performances ont été de limiter le nombre d'itérations lors de la méthode de Gauss-Seidel pour la résolution des systèmes linéaires et la limitation du nombre de points dans un contour.

On peut noter qu'il reste l'option de paralléliser le code sur GPU (avec CUEDA par exemple) pour obtenir de meilleurs performances, ce n'est cependant pas un point que nous avons abordé.

L'interface graphique a été réalisé sous Qt avec Qtdesigner qui nous a permit de créer rapidement une interface fonctionnelle sans connaissances préalable de Qt.

Interface utilisateur

Bilan

Notre application donne des rendus très réalistes et a une prise en main aisée pour un utilisateur lambda. Ceci la rend efficace pour produire des images voulues par l'utilisateur.

Bien que l'application a pour principal but d'être utilisée comme un outils permettant de faire des dessins ressemblant à des Suminagashis, l'ensemble du travail réalisé peut être adapté et complété par des méthodes permettant de charger, par exemple, des simulations et être utilisé en milieu industriel pour simuler des phénomènes réels.

Sources

Installation

Manuel utilisateur

Sources

Bibliographie

[1] Jos Stam Real-Time Fluid Dynamics for Games, Game Developer Conference (2003)<br\> [2] Jos Stam Stable Fluids, SIGGRAPH 99 Conference Proceedings, Annual Conference Series (1999)
[3] Ryoichi Ando et Reiji Tsuruno Vector Graphics Depicting Marbling Flow, Computers & Graphics. Volume 35, (2011)
[4] Ryoichi Ando et Reiji Tsuruno Vector Fluid: A Vector Graphics Depiction of Surface Flow, 8th International Symposium on Non-Photorealistic Animation and Rendering (2010)