Simulation et rendu d'un fluide

De Ensiwiki
Aller à : navigation, rechercher


Pres.jpg Simulation et rendu d'un fluide
Projet Projet de spécialité 2A
Sous-projet Image
Étudiants Émilie BARRE (MMIS-MCS), Pierre JOLIVET (MMIS-MCS), Antoine MILLIEZ (MMIS-IRVM)
Promo 2011
Encadrants Marie-Paule CANI Marie-Paule.Cani@grenoble-inp.fr, Christophe PICARD Christophe.Picard@grenoble-inp.fr

Introduction

De nos jours, les problèmes liés à la mécanique des fluides sont observés de manière très pointue par les chercheurs en mathématiques appliquées et par les physiciens. Ceci s'explique par le fait que ces mêmes problèmes interviennent dans nombre de phénomènes de la vie courante ou du monde industriel :

  • impact du profil d'un avion sur sa consommation
  • impact des lignes d'une voiture sur ses performances
  • ou bien encore, impact d'une fuite dans une canalisation haute pression en mer

Comme nous allons le voir par la suite, de nombreuses méthodes existent pour résoudre des problèmes fluides, mais de ces différentes méthodes émerge un problème commun, l'équation de Navier-Stokes (et ses possibles simplifications dans le cadre d'un fluide incompressible par exemple) :

\frac{\partial \left( \rho \vec{v} \right)}{\partial t} + \overrightarrow{\nabla} \cdot \left(\rho \vec{v} \otimes \vec{v} \right) = - \overrightarrow{\nabla} p + \overrightarrow{\nabla} \cdot \overrightarrow{\overrightarrow {\tau}} + \rho \vec{f}

  • t représente le temps
  • \rho désigne la masse volumique du fluide
  • \vec{v} désigne la vitesse eulérienne d'une particule fluide
  • p désigne la pression
  • \overrightarrow{\overrightarrow{\tau}} = \left( \tau_{i,j} \right)_{i,j} est le tenseur des contraintes visqueuses (en \rm Pa)
  • \vec{f} désigne la résultante des forces massiques s'exerçant dans le fluide (en \rm N.kg^{-1})
Photo d'une fuite de pétrole en mer

Notre but pour ce projet a été de modéliser une fuite de pétrole en mer, à l'instar des ingénieurs [BP].

Les différentes approches

Il existe principalement deux approches pour évaluer les propriétés physiques telles que la vitesse d'un fluide :

  • l'approche lagrangienne qui consiste à suivre l'évolution des trajectoires de chaque particule. Pour toute particule située initialement en (x_0, y_0, z_0), on connaît pour tout t les nouvelles coordonnées (x,y,z) = f(x_0, y_0, z_0, t)
  • l'approche eulérienne qui ne tient compte que du champ de vecteurs des vitesses, en tout point (x,y,z), et pour tout instant t, on connait le vecteur vitesse V(x,y,z,t)


Les méthodes numériques adaptés

Il semble au premier abord difficile d'élaborer une méthode pour l'approche lagrangienne étant donné que cela impliquerait de stocker la position et la vitesse de chaque élément du fluide, par exemple des particules infinitésimales. Nous allons donc dans un premier temps discuter des méthodes pour l'approche eulérienne.

La méthode la plus simple reste l'utilisation d'une grille eulérienne, i.e. on discrétise l'espace et les opérateurs utilisés dans l'équation de Navier-Stokes puis grâce à des schémas numériques on se ramène à la résolution d'un grand système linéaire. Une méthode plus sophistiquée à savoir la méthode Volume of Fluid permet de résoudre les problèmes à surface libre. Elle a fait ses preuve et peu s'adapter à beaucoup de contextes. Le but est de localiser et observer la déformation de la surface du fluide. Pour ce faire, on doit encore une fois discrétiser le volume d'étude selon une grille. On définit alors une fonction "rapport" \nu qui est égale à l'intégrale de la fonction caractéristique du fluide sur une unité élémentaire du volume. Par exemple, si dans un voxel de la grille, il n'y pas de fluide, on a alors \nu = 0, si par contre le voxel est plein de fluide, \nu = 1.

Si on appelle V le champs des vitesses, alors on obtient grâce à la dérivée matérielle de \nu l'équation d'Hamilton-Jacobi suivante :
\dfrac{\partial \nu}{\partial t} + V \cdot \nabla \nu = 0

qui ressemble fortement à une équation de type level-set (le lecteur averti aura remarqué la similitude entre les deux problèmes). Malheureusement, \nu n'est pas continue, on ne peut donc pas appliquer des méthodes numériques du type schéma implicite Euler d'ordre 1. Les méthodes les plus efficaces pour ce genre de problème sont celles liées à la reconstruction géométrique.


Après avoir observé quelques méthodes eulériennes, nous allons maintenant parler d'une méthode lagrangienne qui a été l'objet d'une étude plus approfondie de notre part. Il s'agit de la méthode Smoothed-Particle Hydrodynamics (SPH). Il est nécessaire pour l'appliquer de discrétiser le fluide avec des éléments que l'on appelle particules, ayant un rayon h. On a besoin de fonctions dites noyaux qui dépendent de la distance entre une certaine particule et ses voisines. Lorsque l'on s'écarte trop d'une particule, ces fonctions sont nulles. Elles vont permettre de calculer les grandeurs physiques au voisinage d'une particule, comme la densité, ou la pression, pour en déduire par l'évolution en position des particules. On gagne ainsi un temps précieux en ne calculant pas les contributions pour ces mêmes grandeurs de toutes les particules, sous reserve d'utiliser un algorithme de recherche de voisin adapté. Mathématiquement, si on note W la fonction noyau, on peut estimer la quantité F en un point M_i = (x_i, y_i, z_i) par
F(M_i) = \sum_{j\mbox{ }\mathrm{une}\mbox{ }\mathrm{particule}} m_j \dfrac{F(M_j)}{\rho_j} W(M_j - M_i, h)

En pratique, cette somme se limite à un nombre de particules inférieur à 100.

Si on cherche la masse volumique \rho, la formule précédente devient
\rho(M_i) = \sum_{j\mbox{ }\mathrm{une}\mbox{ }\mathrm{particule}} m_j  W(M_j - M_i, h)
L'équation continue à discrétiser pour trouver le déplacement des particules est
P = k \cdot (\rho - \rho_0)
P est la pression, k une constante et \rho la masse volumique. On cherche à trouver une pression constante dans tout le fluide, il s'agit donc d'une relaxation des particules en fonction de leur densité.

Visualisation 3D

Étude des SPH

Nous allons dans cette partie detaillé le modèle des Smoothed-Particle Hydrodynamics. C'est ce modèle que nous avons implementé en utilisant comme référence l'article de Simon Clavet, Philippe Beaudoin, et Pierre Poulin

Description détaillée du modèle

Comme vu en introduction, les SPH permettent d'utiliser un modèle particulaire pour modéliser des problèmes fluides. L'algorithme principale pour avancer d'un pas de temps \Deltat dans le modèle est le suivant :

pour toutes particules i
    // Appliquer la gravité
    vi += \Deltat.g
// Modifier la viscosité avec les impulsions visqueuses mutuelles
AppliquerViscosité
pour toutes particules i
    // Sauvegarder la position précédente
    xiprev = xi
    // Modifier les positions via un schéma eulérien d'ordre 1
    xi += \Deltat.vi
// Relaxation en fonction des densités et des collisions
RelaxationDoubleDensité
CalculCollisions
pour toutes particules i
    // On utilise l'ancienne position pour calculer la nouvelle vitesse
    vi = (xi - xiprev) / \Deltat

Où chaque sous procédures AppliquerViscosité, RelaxationDoubleDensité et CalculCollisions peuvent encore une fois s'écrire en pseudo-code sous les formes suivantes, en ayant au préalable défini u_{ij} comme le vecteur unitaire partant de la particule i et pointant vers j :

AppliquerViscosité RelaxationDoubleDensité CalculCollisions
pour chaque pair de particules i < j voisines
    q = rij / h
    si q < 1
        // Vitesse radiale entrante
            u = (vi - vj) . uij
                si u > 0
                    I = \Deltat (1 - q)(\sigmau + \betau2) uij
                    vi -= I / 2
                    vj += I / 2
pour toutes particules i
    \rho = 0
    \rhovoisinage = 0
    // Mise à jour des deux densités
    pour toutes particules j voisine de i
        q = rij / h
        si q < 1
            \rho += (1-q)2
            \rhovoisinage += (1-q)3
    // Mise à jour des deux pressions
    P  = k (\rho - \rho_0)
    Pvoisinage = kvoisinage.\rhovoisinage
    dx = 0
    pour toutes particules j voisine de i
        q = rij / h
        si q < 1
            // Déplacement des particules
            D = \Deltat2(P.(1-q) + Pvoisinage.(1-q)2)uij
            xj += D/2
            dx -= D/2
    xi += dx
pour toutes particules i dans un solide rigide 
avec un coefficient de friction \mu
    calculer I = vnormal - \muvtangentielle
    // Appliquer \mu à la particule
    // Modifier la position de la particule
       pour l'extraire du solide

Ajustement du modèle

La différence majeure entre le papier cité ci-dessus et l'approche générale des SPH vient de la relaxation par double densité. Nous verrons dans la partie suivante la différence de rendu obtenue en ajoutant cette seconde densité mais nous allons dans un premier temps détailler d'un point de vue théorique la nécessité d'un tel ajout.

Experience réalisée par [Khayyer et al.]

La méthode SPH n'est pas rigoureusement incompressible, on parle en faite de Weakly Compressible SPH (WCSPH), mais on peut cependant, comme l'ont montré Lee, Moulinec et al., l'adapter pour qu'elle le devienne, on parle alors de Incompressible SPH (ISPH), que l'on peut encore affiner pour obtenir des simulations numériques plus justes, on parle alors de Corrected ISPH (CISPH). Sans cet ajout, le risque de clustering devient d'autant plus important que le nombre de particules augmente. Clustering est un terme de Computational Fluid Dynamics (CFD) qui signifie que les particules ont tendance à s'amasser en quelques grosses gouttes.

L'ajout d'une seconde densité permet d'éviter ce phénomène, en faisant contribuer plus de particules pour les calculs des densités. On a tendance à calculer des déplacements plus importants, et donc à mieux répartir les particules dans l'espace.

On aurait pu penser à rajouter une force qui dépend de la distance entre particule, comme une force de Lennard-Jones, mais cette approche ajoute des erreurs de calcul pour des fluides avec une forte viscosité qui ont tendance à se rigidifier.

La deuxième densité permet également la création d'une tension de surface, est permet ainsi pour les fluides assez visqueux d'obtenir des résultats pratiques proches de la réalité.



Nous avons du modifier le modèle de base où les particules sont uniquement soumises à la gravité pour prendre en compte la poussé d'Archimède. Il a également été nécessaire de rajouter une détection des collisions avec un plan rigide à hauteur fixé pour simuler la surface de la mer (au repos, i.e. sans vague).

Travail réalisé

Implémentation de la méthode SPH

Plusieurs problèmes algorithmiques émergent lors de la création des divers algorithmes pour les méthodes particulaires en général.

Le premier que nous avons rencontré est lié au stockage des positions et des vitesses. Nous avons utilisé des std::vector mais cette structure n'est malheureusement pas simple d'utilisation lorsque l'on cherche à paralléliser le code.

La construction de la fonction de hachage revient à trouver une fonction injective de width × height × depth dans \mathbb{N} pour minimiser les collisions

Le problème suivant concerne la recherche des particules voisines. Il existe divers techniques, citons en particuliers les kd-tree, et nous avons choisis l'approche par table de hachage. Pour se faire il est nécessaire de discrétiser l'espace de travail, qui est un volume fini de \mathbb{R}^3 de taille width × height × depth, en N^3 cubes. On définit alors la fonction de hachage comme la fonction qui prend en paramètre des coordonnées (x, y, z) et qui retourne l'entier hash(x, y, z) = x / N + width / N × y / N + width / N × height / N × z / N. On associe alors à chaque clef un std::set qui regroupe toutes les particules ayant la même valeur de hachage.

Pour chercher les voisins d'une particule située en (x0, y0, z0), il ne reste plus qu'à obtenir de la table de hachage les std::set en hash(x0 + i, y0 + j, z0 + k) avec (i, j, k) \in \{-1, 0, 1\}.

En C++, les tables de hachage sont implémentées dans le Technical Report 1 en tant que std::tr1::unordered_map, encore une fois, l'utilisation de cette structure n'est pas thread-safe.

Comparaison des approches

Les 2 images ci-dessous nous on permis d'étudier les 2 modèles, on voit bien que lorsqu'il n'y a pas de double densité, il y a formation de 2 tas indépendants de particules, tandis qu'avec la double densité, on a une répartition uniforme des particules. Dans les 2 cas, on part de la même condition initiale : les positions des particules sont les mêmes au début.


Visualisation

La visualisation d'un système de particules peut se faire de deux façons différentes :

  • Représentation des points dans l'espace.
  • Représentation d'une surface implicite dont l'ensemble des particules est le squelette.

Comme on peut s'en douter, un simple rendu par points n'est pertinent que si le nombre de particules est fort élevé. Cependant même avec une machine très puissante pouvant représenter plusieurs dizaines de milliers de particules, un tel rendu n'était pas satisfaisant pour représenter du pétrole. Nous avons donc adopté une représentation par surface implicite, décrite par une équation type f(x,y,z) = iso.

La représentation d'une surface implicite peut se faire grâce à deux techniques majeures :

  • Le marching cube : l'espace est discrétisé en une grille 3D, pour chaque cube de cette grille, on calcule les intersections des arêtes avec la surface, et en fonction des différents cas possibles on effectue un rendu dans chaque cube.

Cette méthode est extrêmement gourmande en calculs, et n'est pas très adaptée dans notre cas, car les particules de pétrole vont essentiellement se concentrer en une colonne évasée. Beaucoup de zones de l'espace sont donc vides et il est inutile de faire autant de tests sur chaque case vide.

  • Le lancer de rayon : des rayons lumineux sont lancés depuis l'oeil de l'observateur de la scène, et on étudie les points d'intersection de ces rayons avec la surface.

Nous avons adopté cette méthode, elle sera détaillée plus bas.

Surface implicite

Une surface implicite est une surface décrite par une équation du type f(x,y,z) = iso. Nous avons utilisé pour f la fonction de Wyvill f(d) = -\frac{4}{9}d^6 + \frac{17}{9}d^4 - \frac{22}{9}x^2 + 1 pour 0 < d < 1. Où d = \sqrt{x^2 + y^2 + z^2}.

Représentation de la fonction de Wyvill

A un point de l'espace, cette fonction va associer une valeur proche de 1 si le point est très proche d'une particule, et une valeur plus faible si on s'éloigne de la particule. Pour chaque point de l'espace, on somme les contributions des fonctions de Wyvill en chaque particule, et si la somme dépasse l'iso défini, le point de l'espace est à l'intérieur de la surface implicite.

Cette approche permet de construire des Blobs (pour Blinn Objects, élaborés par Blinn en 1982). Une particule isolée va générer une boule centrée en elle-même, mais deux particules proches vont générer deux sphères "rejointes" comme par capillarité (voir illustration). Ce modèle est donc très bien adapté à la représentation d'un fluide constitué de gouttes.

Raycasting

Vision dans l'espace

Pour représenter ces Blobs, nous avons utilisé la technique du raycasting. Pour faire simple, le raycasting consiste à jeter "à travers chaque pixel" de l'écran un rayon provenant de l'oeil de l'observateur, et d'étudier les objets rencontrés par ce rayon.

Nous définissions tout d'abord l'emplacement de l'oeil (ou de la caméra) dans l'espace, ainsi que la direction dans laquelle elle regarde (en rouge sur l'illustration).

La scène est ensuite observée à travers une fenêtre de l'espace. Il faut définir cette fenêtre. Nous avons décidé de la placer à une distance fixe de la caméra. Afin de jouer sur le zoom, il faut donc définir l'angle de vue représenté sur le schéma, qui définira la taille de cette fenêtre.

Toutes ces informations ne nous permettent jusqu'à présent que de définir le plan de l'espace contenant la fenêtre d'observation. Il est donc également nécessaire de définir le vecteur de l'espace décrivant la direction vers le haut. Ce vecteur est ensuite projeté sur le plan en question, et permet de définir l'orientation de la fenêtre dans le plan.

Tous les éléments sont définis! Il reste à discrétiser cette fenêtre en fonction de la résolution de l'image désirée. Le vecteur oeil->pixel(i,j) est ensuite calculé et définit la direction du rayon.

On avance ensuite à pas constant le long de chaque rayon. Si un point du rayon permet de dépasser l'iso définie, alors c'est qu'on est rentré dans la surface à représenter, et une dichotomie permet de déterminer l'emplacement exact de l'intersection entre la surface et le rayon.

Si le rayon atteint une profondeur limite que l'on a définie, sans avoir rencontré la surface, alors on estime qu'il a atteint le fond de la mer et une couleur bleue est renvoyée.

Afin de définir le sable et la surface, on définit la profondeur du fond de la mer et la hauteur de la surface. Si un rayon atteint l'un de ces deux plans on renvoie : une couleur bleue éclairée pour la surface, et une nuance d'orange pour le sable. Cette nuance d'orange est calculée en fonction de sa localisation afin de donner l'impression d'une texture de sable.

Une fois qu'un rayon atteint le sable, on le fait remonter à la verticale. Si un rencontre un Blob, une couleur assombrie est renvoyée pour le sable. Si en revanche la surface est atteinte sans croiser de blob, on renvoie la couleur du sable éclairé. Ceci nous permet de définir une ombre portée.

Rendu final


Conclusion

Ce projet nous a permis de mettre en application beaucoup de notions très utilisées dans le monde de l'analyse numérique et de la représentation graphique. L'implémentation de A à Z de l'ensemble des deux notions précédentes nous a forcé à mettre en place une démarche scientifique adaptée à un travail de grand envergure et de développer nos connaissances entrevues dans nos cours respectifs en MMIS-MCS et MMIS-IRVM.

Évolution possible

  • Meilleure parallélisation (actuellement, nous travaillons sur CPU à l'aide de QThread, on peut penser à du calcul GPU très adapté pour le raytracer entre autre à l'aide de CUDA)
  • Adaptation du modèle pour être fiable pour les fluides incompressibles

Slides

Disponible ici

Bibliographie

  • Particle-based Viscoelastic Fluid Simulation de Simon Clavet, Philippe Beaudoin et Pierre Poulin
  • Modéliser les évacuateurs de crue avec la méthode numerique SPH de Eun-Sug Lee, Damien Violeau, Réza Issa, Guillaume Thibault, Raphael Marc
  • Visualisation de calculs hydrodynamiques SPH de EDF
  • Analyse structurale d'un barrage flottant anti-hydrocarbure par élèments finis de Frédéric Muttin et Serge Nouchi
  • Corrected SPH for incompressible fuid for accurate water-surface-tracking in plunging breaker par Abbas Khayyer