Utilisateur:Vingtsens/brouillon
===========================================================================
modifierStable-Fluid
modifierLe Stable-Fluid est une méthode de simulation de fluide (CFD) décrite par Jos Stam en 1999 dans le cadre de la conférence SIGGRAPH99 [1].
Elle utilise une grille eulerienne: un espace composé exclusivement de fluide où le mouvement du fluide est décrit par un champs de vecteurs vitesse (aussi nommé vélocité).
Formulation mathématique
modifierEquations de Navier-Stokes
modifierLe modèle du Stable-Fluid utilise les équations de Navier-Stokes dans le cas des fluides incompressibles.
- Équation d'incompressibilité
- Équation de bilan de la quantité de mouvement
Ces équations peuvent être interprétées de la manière suivante:
- Le terme représente le vecteur vitesse à un point donné du fluide.
- Les termes d'accélérations correspondent à:
Où le second membre est nommé advection.
- Les forces de viscosité correspondent au terme:
- Les forces de pressions correspondent au terme:
- Les forces extérieures correspondent au terme .
- L’incompressibilité se traduit par la relation:
Théorème de décomposition de Helmholtz-Hodge
modifierA partir de là, on cherche à simplifier les équations de Navier et Stokes en utilisant le théorème de décomposition de Helmholtz-Hodge. De cette manière on inclus le terme de projection .
De plus, ce théorème nous permet de supprimer plusieurs éléments de l'équations, comme l'équation de l'incompressibilité et le terme relatif à la pression:
Ainsi on arrive à l'équation suivante:
Le modèle numérique
modifierMaintenant, au niveau de l'implémentation, les différents composants (diffusion, advection, projection, forces extérieures) ne sont pas calculés puis additionnés ensemble contrairement à sa forme mathématique. On considère plutôt que chaque composant, à partir d'un champs de vecteurs d'entrée, effectue des modifications sur le champs et donne le résultat au composant suivant.
On définit des opérateurs pour l'advection (), la diffusion (), l'application des forces extérieures () et la projection () et opérateur équivalent à la formulation mathématique pour un pas-de-temps.
Ainsi, à chaque pas-de-temps, le champs de vecteur u passe successivement par l'opérateur d'advection puis ceux de diffusion, des forces extérieures et en dernier par l'opérateur de projection.
Implémentation
modifierPour l’implémentation de base 2D (peu de différence pour la 3D), considérons une grille régulière de N*N cellules, où chaque cellule est définit par un vecteur vélocité (au centre de la cellule) et une pression (pour la partie projection).[2]
Les données seront représentées dans des tableaux ou i et j correspondent à la position de la cellule courante dans la grille.
Diffusion
modifierL'étape de diffusion correspond à la résolution de la viscosité. La vitesse du fluide va impacter les cellules voisines. Pour cela, un solveur linéaire sera utilisé.
Si correspond à la vélocité d'une des cellules et , , et aux les cellules voisines. Alors le solveur doit résoudre, le système où pour chaque cellule:
où
: le pas de temps.
: coefficient de diffusion/viscosité.
: nombre de cellules => correspond à .
Pour la résolution du système, le solveur du Gauss-Seidel peut être utilisé. Dans cet exemple, on travaille dans un tableau où seul les voisins directs sont utilisés, les autres cellules n'apparaissent pas et ont comme coefficient 0.
int N = 100 // largeur et hauteur du tableau int kmax = 100 // nombre d'itérations à effectuer (solveur) float diff = 0.01 // coefficient de diffusion float u[N,N] // composante horizontale des vecteurs vélocité float v[N,N] // composante verticale des vecteurs vélocité float u0[N,N], v0[N,N] // tableaux temporaires relatifs à u et v // solveur linéaire Gauss-Seidel void linear_solver (int b, float * x, float * x0, float a, float c, int kmax ) { int i, j, k; for ( k=0 ; k<kmax ; k++ ) { FOR_EACH_CELL x[IX(i,j)] = (x0[i,j] + a*(x[i-1,j]+x[i+1,j]+x[i,j-1]+x[i,j+1]))/c; END_FOR set_bnd (b, x ); // conditions limites } } // opérateur de diffusion void diffuse_step() { SWAP ( u0, u ); SWAP ( v0, v ); float a=_fm->dt()*diff*N*N; lin_solve (1, u, u0, a, 1+4*a,kmax ); lin_solve (2, v, v0, a, 1+4*a,kmax ); }
Advection
modifierDurant l'étape d'advection, on considère donc qu'il y a une particule au niveau du vecteur vélocité, on effectue un backtracking pour déterminer la position de la particule au pas-de-temps précédent, enfin on effectue une interpolation pour connaitre la vélocité de la particule. Cette dernière remplacera la valeur du vecteur vélocité. On cherche donc à garder l'énergie contenue dans la particule.
advection (int b, float * d, float * d0, float * u, float * v) { int i, j, i0, j0, i1, j1; float x, y, s0, t0, s1, t1, dt0; dt0 = _fm->dt()*_n; //coefficient FOR_EACH_CELL // position de la particule au pas de temps précédent x = i-dt0*u[IX(i,j)]; y = j-dt0*v[IX(i,j)]; // condition limites, les bords de la grille if (x<0.5f) x=0.5f; if (x>_n+0.5f) x=_n+0.5f; i0=(int)x; i1=i0+1; if (y<0.5f) y=0.5f; if (y>_n+0.5f) y=_n+0.5f; j0=(int)y; j1=j0+1; // interpolation s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1; d[IX(i,j)] = s0*(t0*d0[IX(i0,j0)]+t1*d0[IX(i0,j1)])+ s1*(t0*d0[IX(i1,j0)]+t1*d0[IX(i1,j1)]); END_FOR // conditions limites, bord de la grilles set_bnd (b, d ); } advection_step() { SWAP ( u0, u ); SWAP ( v0, v ); diffuse_vel_x (1, u, u0, _visc,this->_algoVel_diffuse_k); diffuse_vel_y (2, v, v0, _visc,this->_algoVel_diffuse_k); }
Projection
modifierLa projection permet de rétablir la conservation de masse et donc l'incompressibilité du fluide. Pour cela, on doit déterminer les quantités de matière entrante et sortante au niveau de chaque cellule, puis corriger les vitesses de telle façon à avoir autant de matière entrante que sortante.
La divergence correspond à la différence entre les flux entrants et sortants. Pour cela on prend la composante horizontale des vecteurs voisins situés à l'horizontal et la composante verticale des vecteurs voisins situés à la verticale.
: composante horizontale du vecteur vélocité
: composante verticale du vecteur vélocité
project (float * u, float * v, float * p, float * div, int kmax ) { int i, j; //Calcul de la divergence FOR_EACH_CELL div[i,j] = -0.5f*(u[i+1,j]-u[i-1,j]+v[i,j+1]-v[i,j-1])/_n; p[i,j] = 0; END_FOR set_bnd (0, div ); set_bnd (0, p ); //Résolution du système lin_solve (0, p, div, 1, 4, kmax); //Mise à jour des vélocités FOR_EACH_CELL u[i,j] -= 0.5f*_n*(p[i+1,j]-p[i-1,j]); v[i,j] -= 0.5f*_n*(p[i,j+)]-p[i,j-1]); END_FOR set_bnd (1, u ); set_bnd (2, v ); }
Conditions aux limites
modifierExtensions de la méthode
modifierStaggered-Grid
modifierDans une staggered-grid[3] les variables scalaires (comme la pression) sont stockées aux centre des cellules, tandis que les vélocités sont localisées sur les bords (ou faces) de la cellule.
Utiliser une Staggered Grid évite un "odd-even decoupling"[4] entre la pression et la vélocité. On peut connaitre exactement le volume de fluide entrant et sortant de la cellule, ce qui n'est pas le cas si les vélocités sont au centre. Donc la Staggered-Grid permet des calculs plus précis et moins de dissipation numérique.
Le désavantage est que les valeurs sont stockées à différents emplacements, ce qui rend plus difficile le contrôle de volume.
Backtraking
modifierLe backtracking est la méthode qui permet de trouver la position d'une particule dans un fluide eulérien au pas-de-temps précédent. Dans l'advection de base, il s'agit simplement de prendre la vélocité, multiplier par le pas-de-temps et par les dimensions d'une cellule pour retrouver la position précédente de la particule.
Cela dit, le fluide suit un parcours curviligne et la simple utilisation de la vélocité ne suffit pas pour trouver précisément l'ancienne position de la particule. Il peut être effectué des corrections afin d'améliorer la précision. Il y a le Runge-Kutta, le MacCormack[5] et le BFECC[6] (Back and Forth Error Compensation and Correction)
La première figure du schéma ci-contre montre le cas du backtracking simple: on peut remarquer que le backtracking récupère une vélocité (en gris) qui n'appartient pas à la bonne ligne de caractéristique. Le MacCormack permet de se rapprocher de la ligne de caractéristique (dernière figure du schéma, vecteur en gris).
Flux sur une surface arbitraire et Distorsion
modifierDans le cas où on choisit de travailler sur un flux de fluide parcourant une surface, on se retrouve dans le cas où les cellules formant le système ne sont pas des carrés. Ces déformations doivent être prises en compte pour éviter des phénomènes irréalistes[7].
Surface Libre
modifierLa notion de Surface Libre concerne la modélisation de la surface de l'eau. Il existe différentes méthodes comme le level-set[8] ou le fast and robust tracking of fluid surface[9].
Cellules triangulaires et tétraèdres
modifierDans les cas précédents, ce sont des cellules carrés (ou au moins des quadrilatères) qui ont été utilisés pour les cas 2D et des cubes pour le cas 3D. Il est néanmoins possible d'utiliser des cellules triangulaire[10] ou des tétraèdres.
Vorticité
modifierIl s'agit de remplacer l'étape de projection par une étape où l'on cherche à préserver la vorticité[11]. La vorticité correspond aux phénomènes de tourbillons. Ainsi au lieu de vérifier qu'une cellule contient la bonne quantité de fluide, on vérifie l'incompressibilité par les flux autour de points.
Liens externes
modifierSimulation de l'écoulement d'un fluide - Humbert Florent - 2008
Fast Fluid Dynamics Simulation on the GPU - GPU Gems - NVIDIA
Simulation de fluide sur les surfaces - Ensiwiki
Fluid Simulation for Visual Effects - Magnus Wrenninge - 2003
Notes
modifier- "Stable Fluid" - Jos Stam - 1999
- Real-Time Fluid Dynamics for Games
- "Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface" - Harlow, F. H. and Welch, J. E. - 1965,
- "Numerical Heat Transfer" - SV Patankar
- [http://physbam.stanford.edu/~fedkiw/papers/stanford2006-09.pdf "An Unconditionally Stable MacCormack Method"]
- "Using BFECC for Fluid Simulation"
- "Flows on Surfaces of Arbitrary Topology"
- "level-set method for fluid interface"
- "Fast and Robust Tracking of Fluid Surfaces"
- "Inviscid and Incompressible Fluid Simulation on Triangle Meshes"
- "Stable, Circulation-Preserving, Simplicial Fluids"