Simulation de fluides
Tuteur: Colin Weill-Duflos
Etudiant : Rousseau Maxime
Lorsque l'on veut afficher de l'eau sur un écran, on veut généralement représenter son mouvement sans qu'un animateur n'ait à préciser à la main le déplacement de chaque polygone représentant la surface de l'eau. Pour cela, on emploie des modélisation physiques du fluide pour simuler son comportement, puis on effectue un rendu du fluide ainsi modélisé. Les modèles physiques employés pour simuler un fluide sont complexes et font apparaître des phénomènes comme des vortex, des vagues... Il faut donc des méthodes de simulation fines capturant la complexité de ces phénomènes.
Introduction: affichage, schéma d’Euler
On utilisera python, avec la bilbiothèque numpy pour créer et manipuler des tableaux, et matplotlib pour afficher des choses.
Commencer par afficher une grille 30x30 en 2d de particules (des points), occupant un coin de l’écran.
Simulation d’un modèle physique simple
Commençons par simuler la chute de chaque particule comme si celle-ci étaient juste des masses en chute libre.
Construire un tableau pour la position de chaque particule (si ce n’est pas déjà fait), la vitesse et l’accélération. Initialiser la position pour avoir la même grille qu’avant, et la vitesse ainsi que l’accélération à 0. Puisque l’on est en 2 dimensions, on peut choisir de représenter position, vitesse, accélération comme des tableaux à 2 dimensions ou bien deux tableaux chacun.
On suppose que chaque particule est soumise à un la gravité.
Rajouter une boucle dans laquelle:
- on fait avancer une variable de temps t d’un pas de temps dt
- on actualise la valeur de l’accélération, puis la vitesse à partir de l’accélération à partir de la position
- on actualise l’affichage avec les nouvelles positions des particules
Pour calculer la vitesse en fonction de l’accélération, et la position en fonction de la vitesse, on se rappelle d’abord le lien entre chaque grandeur, puis on utilise la méthode d’Euler explicite. ( https://fr.wikipedia.org/wiki/Méthode_d%27Euler )
Conditions aux bords
Rajouter des contraintes au bord, simulant un rebond. Chaque particule devrait rebondir jusqu’à sa position initialie. Rajouter une vitesse initiale horizontale non nulle, afin de vérifier que le rebond fonctionne bien sur tous les bords.
On a à présent les conditions aux bords ainsi que la boucle d’actualisation accélération + vitesse + position qui nous serviront par la suite. La difficulté réside à présent dans le calcul des forces afin de déterminer l’accélération.
Partie 1 : le kernel, calcul de quantité physique
Utilisation d’un kernel
On va employer une fonction de kernel appelée W. Cette fonction possède plusieurs propriétés :
- elle vaut 0 au delà d’un rayon h
- son intégrale vaut 1
- elle a une allure en cloche, est lisse et ses “dérivées” sont lisses
- elle ne dépend que de la distance entre deux points
On utilise cette fonction pour passer de la description “discrète” du milieu (seulement défini sur quelques points) à une description permettant d’obtenir une valeur dans n’importe quel point du plan.
On utilise cette fonction est utilisée pour faire une moyenne des valeurs des particules voisines. Avoir une fonction lisse permet d’obtenir un résultat lisse également, avec des variations lisses.
Ici, pour calculer la valeur au point p, on fait la moyenne des valeurs aux particules à une distance ≤h de p. La valeur du kernel est indiquée en rouge: les particules les plus proches vont avoir le plus d’influence sur la valeur totale en p.
On veut également faire contribuer les particules selon leur volume: les particules avec le plus de volume vont contribuer le plus.
Ainsi, on obtient une formule pour obtenir en un point p quelconque la valeur d’une grandeur physique quelconque A :
A(p)=∑(j) A(pj) W(||p−pj||) V(j)
Avec:
- A(j) la valeur de la grandeur physique A associée à la particule j, que l’on suppose connue
- V(j) le volume occupé par la particule j
- ||p−p(j)|| la distance entre le point p et la particule p(j)
Plutôt que de manipulé V(j), on passe par la masse et la masse volumique ρ(j) :
V(j)=m(j)/ρ(j)
Ici, on suppose que chaque particule à la même masse m qui ne varie pas au cours du temps. La formule devient :
A(p)=∑(j) A(p(j)) m/ρ(j) W(||p−p(j)||)
Pour appliquer cette formule, on a besoin de connaître la masse volumique ρj de chaque particule. Pour cela, il suffit d’employet la formule! La grandeur physique A devient ρ, ce qui nous donne :
ρ(p)=∑(j) ρ(j) m/ρ(j) W(||p−p(j)||)
ρ(p)=∑(j) m W(||p−p(j)||)
Cette formule, bien qu’apparemment ayant besoin de connaître ρ, n’en a pas besoin ici: on peut donc l’utiliser pour calculer ρ.
Au début de chaque boucle, on actualisera la valeur de ρ pour chaque particule.
Calcul du kernel
Le kernel est une fonction qui “mange” une distance d, et doit avoir les propriétés définies précédemment. Il y a plusieurs choix de fonctions, nous n’irons pas dans les détails de chacune.
On définit d’abord une distance h correspondant au rayon de notre kernel: ici, on peut essayer h=0.02.
On peut ensuite calculer q=d/h : q sera égal à 1 pour d=h (rappel: d correspond à une distance).
On utilise ensuite le cubic spline kernel :
W(q)=σ(1−1.5q2+0.75q3) pour 0≤q≤1=σ4(2−q)3 pour 1≤q≤2=0 sinon
La constante σ dépend de la dimension: ici, en dimension 2, on a σ=107πh2
On va également avoir besoin de la dérivée W′ :
W′(q)=σ(−3q+2.25q2) pour 0≤q≤1=−0.75σ(2−q)2 pour 1≤q≤2=0 sinon
Partie 2 : calcul des forces
Force de pression
On veut calculer une force de pression, principale cause de mouvement dans notre fluide.
Cette force est donnée par : fp= −1/ρ(i) ∇ p(i), avec rho(i) la densité de la particule i et ∇p(i) le gradient de la pression de la particule i (le gradient est un objet proche de la dérivée, généralisale à plusieurs dimension).
Calcul de gradient
Pour calculer le gradient de la pression, on peut reprendre la formule précédemment:
A(x)= ∑(j) A(x(j)) m/ρ(j) W(||x−x(j)||)
(on note x la position au lieu de p précédent afin de ne pas confondre avec la pression)
On veut calculer le gradient de A:
∇A(x) = ∇(∑(j) A(x(j)) m/ρ(j) W(||x−x(j)||))
On peut suivre les mêmes règles que pour la dérivée: les termes constants (ne dépendant pas de x) sortent:
∇A(x) = ∑(j) A(x(j)) m/ρ(j) ∇W(||x−x(j)||))
Il nous faut donc la possibilité de calculer le gradient de notre fonction “kernel”.
En suivant des règles similaires à celle de la dérivée d’une fonction composée ((f∘g)′(x)=g′(x)∗f′(g(x))) on obtient:
∇W(||x−x(j)||)) = W′(||x−x(j)||) ∗ ∇(||x−x(j)||)
||x−x(j)|| correspondant à d, le premier terme est donnée par la dérivée du kernel. On admet que :
∇(||x−x(j)||) = x−x(j) / ||x−x(j)||
Calcul de la pression
Il nous reste à connaitre la pression en chaque particule. La physique nous donne une relation:
p(i) = k((ρ(i) ρ(0))7−1) avec k et ρ0 deux constantes.
Ce calcul peut donner une pression négative, ce qui est à éviter: si le calcul donne une pression négative, on utilise 0 à la place.
Puisque l’on connait ρ(i) pour chaque particule, il suffit d’appliquer la formule pour connaître p(i).
On peut donc calculer la force de pression:
−1/ρ(i)∇p(x(i)) = −1/ρ(i) ∑(j) p(j) m/ρ(j) ∇W(||x(i)−x(j)||))
Cette formule ne respecte pas la 3e loi de Newton (principe de réciprocité): on lui préfère une version plus verbeuse:
−1/ρ(i)∇p(x(i)) = −∑(j) m(p(i)/ρ(i²) + p(j)/ρ(j²)) ∇W(||x(i)−x(j)||))
On a ainsi une nouvelle force. Suivant la 2e loi de Newton, cette force divisée par m donne un terme de accélération (auquel on ajoute le g de la gravité que l’on utilisait déjà).
Cette force toute seule doit amener les particules à se repousser, et elle peut donner des résultats un peu explosifs.
Viscosité
On rajoute une nouvelle force, qui va “freiner” les particules et va permettre de corriger le côté explosif de la pression. Cette force tend à faire en sorte que toutes les particules aient la même vitesse.
Elle est normalement donnée par la formule: fvisco= νΔv
Où ν est une constante appelée viscosité du fluide et Δv est le “laplacien” de v (v est la vitesse). On ne s’attardera pas sur le laplacien, qui peut être vu comme étant à la dérivée seconde ce que le gradient est à la dérivée.
On pourrait reproduire les calculs pour le gradient avec le laplacien et en déduire qu’il faut calculer le laplacien de notre kernel, mais pour corriger certains effets indésirables on préfère la formulation suivante de la viscosité:
f(i)viscosity = 2νm ∑(j) m/ρ(j) (v(j)−v(i)) [(x(j)−x(i)).∇W(||x(i)−x(j)||) / ||x(i)−x(j)||2+0.01h2]
Où h est le rayon de notre kernel (donc constant et déjà défini), x(i) désigne la position de la particule i et v(i) sa vitesse.
On peut rajouter cette force (divisée par m encore) pour obtenir notre dernier terme de calcul l’accélération.
Les particules devraient avoir un mouvement plus “unifié”, et on doit pouvoir voir se dessiner une sorte de vague si l’on laisse se dérouler la simulation assez longtemps.
Partie 3 : optimisation du calcul des voisins
Optimisation des calculs
La principale raison des faibles performances de la méthode vient actuellement du nombre de comparaisons à effectuer pour le calcul de chaque quantité : pour chaque particule, on prend en compte toutes les autres particules, même les plus éloignées. Beaucoup de temps de calcul est consacré à tester si la distance entre deux particules est inférieure ou supérieure à h.
On se retrouve avec une complexité en O(n2).
On peut tenter de diminuer ceci. Pour cela, on s’appuie sur une grille avec des cases de taille h×h. En récupérant la case à laquelle appartient une particule, on peut ne faire les calculs qu’avec les autres particules de la case et celles des 8 cases voisines.
Il suffit donc de fabriquer cette grille au début de chaque itération, en associant à chaque case une liste de particule contenue dedans. Pour obtenir la case dans laquelle est contenue chaque particule, il faut obtenir l’index (i,j) de la case (deux entiers). On peut les obtenir ainsi: i=⌊xh⌋, où x est la position en x de la particule (on utilise y pour obtenir j).
On peut aussi obtenir ainsi les dimensions de la grille (plus pratique pour l’initialiser): c’est un tableau de m×m, où m=⌊xlim/h⌋ (xlim est la bordure du cadre).
Optimisation des forces
Pour l'optimisation des forces on sait que rhô prend beaucoup de temps de calcul : l'opération du calcul de rho(j) se fait avec une boucle sur le nombre de particule, donc n itérations (où n est le nombre de particule).
Sachant qu’on recalcule rho(j) à chaque calcul de force, et que des forces il y en a entre chaque paire de particule ,on fait appel à cette fonction n² fois. Or rho ne change pas au cours du calcul, en théorie on peut ne le calculer qu'une fois par particule, donc ne faire appel à la fonction que n fois.
On peut donc faire une première boucle avant la boucle des forces où on calcule tous les rho, et d'après utiliser les valeurs stockées dans le tableau rho_ pour le calcul des forces.
Partie 4 : affichage
Rendu
On veut pouvoir afficher la surface du fluide simulé par les particules.
Il faut donc représenter ce fluide d’une manière compatible à un moteur de rendu. En général, les éléments affichés par un moteur sont des modèles 3d consitués de triangles: on voudrait donc une liste de triangles représentant la surface.
On va appliquer la méthode dites des “marching squares” (https://en.wikipedia.org/wiki/Marching_squares) (une variante “marching cubes” existe en 3d). Cette méthode nécessite une fonction de “distance signée” à la surface. Le principe de cette foncton est de donner la distance entre un point quelconque du sens et le point de la surface le plus proche de ce point. Cette fonction est “signée”, c’est à dire qu’elle est positive à l’exterieur de la surface et négative à l’intérieur du volume délimité par la surface: la surface correspond aux points où cette fonction vaut 0.
Ici, pour construire la distance signée d’un point à notre surface, on peut prendre la distance à la particule la plus proche de ce point à laquelle on soustrait r le rayon de la particule (pouvant valoir 2*h comme pour la simulation physique, ou const*2*h avec const une nouvelle constante).
On construit ensuite une grille de taille arbitraire, plus cette grilel contenant de point et plus la surface étant fine mais chère à calculer. En chaque point de cette grille on évalue notre distance signée, puis on considère chaque “petit carré” formé par cette grille. En utilisant le signe des valeurs de la fonction, on peut savoir si le carré se situe entièrement à l’intérieur ou à l’exterieur du volume définir par la surface. S’il se situe sur la surface elle même, alors il y a des coins du carré ayant des signes différents.
On trace alors des segments dans le carré pour représenter la surface. Pour savoir comment les tracer, on utilise cette disjonction de cas :
Pour placer les sommets des segments à tracer, on peut soit les placer au milieu de l’arête du carré, soit on utilise une interpolation linéaire : on considère que la fonction de distance est une fonction linéaire le long de l’arête du carré contenant le sommet, et vu que l’on connaît ses valeurs aux deux extrémité on peut savoir où cette fonction vaut 0 : c’est là qu’on place le sommet de notre segment.
Autrement dit, si l’on connaît f(0), f(1) et que l’on suppose que f est linéaire (de la forme f(x) = a∗x+b) :
f(0) = a∗0+b =b f(1) = a∗1+b =a+f(0) donc a =f(1)−f(0) f(x) = 0 donne x =−f(0) / (f(1)−f(0))
On peut donc obtenir x indiquant où sur l’arête la fonction s’annule.
En répétant ce processus sur tous les carrés de la grille, on obtient une liste de segments représentants notre surface.
En 3D, on utiliserait des cubes pour produire des triangles.
Ce processus peut se faire à toutes les étapes de la simulation, où bien seulement à certaines (toutes les 2, 5, 10 étapes…).
