Fluid Silhouettes

De Ensiwiki
Aller à : navigation, rechercher

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


Project schedule.png
Titre du projet Fluid Silhouettes
Cadre Projets de spécialité
Page principale Fluid Silhouettes

Encadrants Stéfanie Hahmann


Présentation du projet

Sujet proposé

Étudiants

Clémence MOUTERDE (MMIS-MCS)

Jérémy ROUOT (MMIS-MCS)

Romain VAVASSORI (MMIS-IRVM)

Fluid Silhouettes

Ce projet se base sur la simulation du comportement d'un suminagashi. Ce dernier est un art dans lequel on dépose des tâches sur un liquide. On fait ensuite bouger le liquide, en soufflant dessus par exemple, pour déformer les tâches.

Dans ce projet, nous cherchons donc à représenter des tâches dans un fluide en mouvement. L'objectif premier est de réaliser une interface graphique 2D permettant de représenter la déformation de la tâche dans le fluide de sorte qu'elle adopte des formes artistiques de type marbré.

Prérequis

Nous avons travaillé à partir des articles [1] et [2] qui proposent une simulation basée sur la résolution de l'équation de Navier-Stokes par la méthode des différences finies.

Équation de Navier-Stokes

Nous résolvons l'équation de Navier-Stokes suivante :


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

Avec :

  • u étant la vitesse du fluide
  • \nu la viscosité
  • p la pression du fluide

On se place dans le cas d'un fluide incompressible. En considérant P, un projecteur qui projette tout vecteur en un vecteur de divergence nulle, nous avons :

P \nabla p = 0 et Pu=u

L'équation devient alors:


\frac{\partial u}{\partial t} = P(-(u.\nabla) u+ \nu \nabla^2 u + f)

Résolution par différences finies

Nous pouvons résoudre cette équation grâce aux différences finies en quatre étapes. Ces quatre étapes correspondent aux quatre termes de l'équation ci-dessus. Nous partons de la vitesse initiale w_0 et nous lui appliquons ces quatre étapes qui sont l'addition de force (w_1), l'advection (w_2), la diffusion (w_3), et enfin la projection (w_4 vitesse finale).

Addition de force

Nous résolvons : \frac{w_1 - w_0}{\delta t}= f

Advection

Nous résolvons :  \frac{\partial w_2}{\partial t}= -(w_1.\nabla)w_2

Diffusion

Nous résolvons :  \frac{\partial w_3}{\partial t} = \nu \nabla^2 w_3

Nous utilisons un schéma d'Euler implicite et la méthode de Gauss-Seidel pour résoudre cette équation.

Projection

D'après la décomposition de Helmoltz-Hodge, w_3 s'écrit de manière unique sous la forme : w_3 = w_4 + \nabla q

avec \nabla. w_4 = 0 et q un champ scalaire.

Ainsi, nous avons \nabla.w_3 = \nabla^2 q et w_4 = w_3 - \nabla q


Comme à l'étape précédente, nous utilisons un schéma d'Euler implicite et la méthode de Gauss-Seidel pour obtenir q.

Affinage des points de contour

Les points de contour vérifient l'équation suivante :

 \frac{\partial v}{\partial t}= -(u.\nabla)v

Avec u la vitesse du fluide et v la position du point de contour.

Afin d'éviter les collisions, nous utilisons un schéma de Runge-Kutta à l'ordre 4 que nous notons \phi.

Ainsi, nous avons : v_{t+\Delta t} = \phi(\Delta t, v_t, u_t)

avec  v_t la position d'un point de contour à l'instant t et  u_t la vitesse d'un point de contour à l'instant t.

Ensuite, nous calculons la distance entre les sommets du contour et si cette distance est supérieure à une valeur seuil notée d, nous ajoutons un point dont nous calculons les coordonnées de la manière suivante :

  • On calcule le centre du segment au pas de temps précédent auquel on applique Runge-Kutta (voir l'équation ci-dessous). Ce point à lui seul a l'inconvénient de donner des contour assez anguleux.

v_{new} = \phi(\Delta t, \frac{1}{2}(v_i^{t-\Delta t}+v_{i+1}^{t-\Delta t}), u_i^{t-\Delta t})

  • On calcule le point centre du segment par interpolation cubique pour avoir une meilleure approximation de la courbe. Ce point donne un contour plus courbe mais a tendance à créer des accumulations trop importantes de points dans certaines zones.
  • Au final on prend la moyenne des deux points pour limiter les inconvénients liés aux deux méthodes.


Au contraire, nous supprimons les points dont la distance est inférieure à d/2.

Cependant, cette méthode conduit à une augmentation très rapide du nombre de point. C'est pourquoi, nous calculons pour chaque point de contour sa propre distance seuil d_i en fonction de la courbure, de la turbulence. Ceci nous permet d'éliminer beaucoup de points dans les zones moins importantes.

d_i = d_{max} c_i^{courbure} c_i^{turbulence} c_i^{proximite}+ \epsilon

Avec

c_i^{courbure} = exp(-|K|) (K étant le rayon de courbure)

c_i^{turbulence} = exp(-|\nabla \times u(v_k)|) avec u(v_k) la vitesse du point de contour

 c_i^{proximite} = 1-exp(- d_i^{proximite}) avec d_i^{proximity} la distance entre le point de contour i et son plus proche voisin

Travail réalisé

Choix d'implémentation

Nous avons choisi de réaliser ce projet en C++ pour pouvoir profiter d'un langage objet.

Nous utilisons aussi plusieurs librairies :

  • Qt pour l'interface graphique.
  • Eigen permettant la gestion des vecteurs, matrices et la partie mathématique du projet.
  • OFELI pour la simulation du fluide par éléments finis.

Nous avons organisé le code de façon modulaire. Nous pouvons donc utiliser plusieurs types de solveurs pour la simulation du fluide ainsi que créer plusieurs types de formes pour les tâches.

Digramme des classes

Diag classe.jpeg

Interface graphique

L'interface graphique se compose de plusieurs outils permettant à l'utilisateur de contrôler le fluide ainsi que de pouvoir créer des formes colorées :

Aperçu de l'interface

Les différents outils mis à disposition sont :

  • La création de formes avec le clic gauche. Elles peuvent être prédéfinies (cercles, carrés, triangles) ou libres, dessinées à la main.
  • Le choix des couleurs pour le fond et pour les tâches.
  • La perturbation du fluide avec le clic droit.

On peut se déplacer dans la zone de travail par un clic de la molette de la souris et aussi zoomer avec la molette. De plus il est possible d'afficher le champ de vitesse calculé de même que le contour des tâches, avec des densités de points variantes suivant les zones de l'image.

Les différents paramètres régissant le comportement du solveur et de l'affinage adaptatif sont réglables dans une fenêtre d'options.

Finalement on peut sauvegarder les images obtenues au format SVG pour garder tous les détails des contours.

Résolution par éléments finis

Présentation de la méthode

La résolution par éléments finis est moins intuitive que celles des différences finies au sens où l'on discrétise l'espace des solutions cherchées (on a toujours une discrétisation du domaine). L'idée est d'écrire une formulation variationnelle de l'équation sur un espace de fonctions V hilbertien. On considère alors une suite d'espaces vectoriels V_h de dimension finie "convergeant" vers V.

Par exemple, en dimension 1, V = C^{0}([0,1]) et on construit V_h comme l'ensemble des fonctions continues, polynomiales de degré inférieur à 1 sur chaque intervalle de la subdivision considérée.


Pour résoudre ce nouveau problème, apparemment plus simple (dimension finie), la méthode consiste à définir un des ensembles D, P et \Sigma où :

  • D (compact) est le domaine de résolution.
  • P un espace vectoriel de dimension finie de fonctions de D dans \mathbb{R}.
  • \Sigma un ensemble fini de forme linéaire sur P.

Dans notre cas, D est le pavé unité [0,1]x[0,1]. Les deux autres ensemble sont modulables et vont influer la performance de la méthode :

on peut considérer P=Vect(1, x, y, xy) (fonctions polynomiales de deux variables de degré inférieur 2

et \Sigma = \{ p \longrightarrow p(0,0), p \longrightarrow p(1,0), p \longrightarrow p(0,1), p \longrightarrow p(1,1) \}.


Dans cet exemple, la triangulation du domaine est constituée de rectangles R_i. L'espace V_h représente l'ensemble des fonctions continues sur [0,1]x[0,1] dont la restriction sur tous les rectangles R_i est dans P. Les fonctions de \Sigma permettent d'obtenir les valeurs aux noeuds du maillage (les sommets des rectangles R_i).


En projetant la formulation variationnelle dans la base de V_h, on aboutit à un système linéaire Au=b(1) d'inconnue u. Plus précisément, l'assemblage de (1) (c'est à dire la construction de A et b), consiste à intégrer sur chaque rectangle les fonctions de base de l'espace V_h. Ces calculs d'intégrales peuvent s'approcher par des formules de quadrature faisant intervenir les valeurs de ces fonctions aux noeuds du maillage (formule du trapèze, formule de Simpson ...).

On note ici une force majeure de cette méthode au niveau algorithmique. Les rendus graphiques sont calculés sur des architectures multicoeurs (GPU) qui bénéficient d'une puissance de calcul bien supérieure à celle des CPU. Néanmoins, pour utiliser le plus efficacement cette puissance de calcul des cartes graphiques, on doit avoir une implémentation parallélisable. C'est le cas pour cette méthode de résolution puisque le calcul d'intégrale sur chaque éléments de la maille peut se faire en parallèle. En outre, la résolution du système linéaire de complexité algorithmique cubique en nombre d'opération peut se paralléliser efficacement.

De surcroît, on cherche à avoir une structure simple pour A (tridiagonale par exemple) qui est directement dû au choix de P et \Sigma. On peut ainsi manipuler des structures creuses et donc diminuer l'espace mémoire.

Application

La formulation variationnelle de l'équation de Navier-Stokes fait intervenir des termes non linéaires ce qui ne permet pas d'aboutir directement à un système linéaire.

En mettant la formulation variationnelle sous la forme f(x)=0, on peut appliquer la méthode de Newton. D'autres méthodes consistent à linéariser les termes non linéaires.

Par souci de simplicité, nous avons utilisé des méthodes de résolution données de la librairie OFELI.

Résultats (démonstrations vidéos) et limites

La résolution des différences finies en procédant suivant l'article [1] permet d'obtenir une simulation réaliste du phénomène étudié (cf vidéo 1 et vidéo 2) ainsi que des rendus proches de l'art suminagashi :

Rendu1.png


Rendu2.png

Au niveau des performances, elles restent limitées pour des ordinateurs personnels (typiquement double coeur cadencé à 2GHz avec carte graphique GeForce 8600M) où la limite de 5 frames par seconde est atteinte avec 3000 points de contours. Néanmoins, la taille de la grille des vitesses est influençable et peut être modulée dans notre application. Il convient de choisir une grille adaptée à la taille de l'affichage.


Nous avons rencontrés des difficultés sur la résolution par la méthode des éléments finis en utilisant des routines de résolution de la librairie OFELI : nous n'arrivions pas à simuler un comportement similaire du fluide, ceci exclue toute comparaison.

Conclusion et perspectives

Nous avons implémenté la simulation planaire d'une tache dans un fluide suivant l'article [1] et proposé une amélioration sur l'affinage des contours. On fournit en outre une interface graphique qui permet à l'utilisateur d'obtenir des rendus proche de l'art suminagashi.

Pour l'optimisation algorithmique, il serait intéressant d'analyser la parallélisation du code dans l'optique de pouvoir l'implémenter sur GPU pour obtenir des rendus plus complexes.

Références

[1] Stable fluids, Stam, J., Proceedings of the 26th annual conference on Computer graphics and interactive techniques (p.121-128), 1999.

[2] Vector Fluids, Ryoichi Ando and Reiji Tsuruno. Proc. Non-Photorealistic Animation and Rendering, 2010.

[3] Analyse numérique et optimisation: une introduction à la modélisation mathématique et à la simulation numérique, Allaire, G., Ecole polytechnique, 2005.