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.

Simulation de fluides : partie 1

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

Simulation de fluides : partie 2

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.