Stable-Fluids

Le Stable-Fluids est une méthode de simulation de fluide (MFN) 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 champ de vecteurs vitesse (aussi nommé vélocité).

Formulation mathématique

Équations de Navier-Stokes

Le modèle du Stable-Fluid utilise les équations de Navier-Stokes dans le cas des fluides incompressibles.

Équation d'incompressibilité
${\overrightarrow {\nabla }}\cdot {\vec {v}}=0$
Équation de bilan de la quantité de mouvement
${\frac {\partial {\vec {v}}}{\partial t}}=-\left({\vec {v}}\cdot {\overrightarrow {\nabla }}\right){\vec {v}}-{\frac {1}{\rho }}{\overrightarrow {\nabla }}p+\nu \nabla ^{2}{\vec {v}}+{\vec {f}}$

Ces équations peuvent être interprétées de la manière suivante :

Le terme ${\vec {v}}$ représente le vecteur vitesse à un point donné du fluide.
Les termes d'accélérations correspondent à :
${\frac {\partial {\vec {v}}}{\partial t}}+\left({\vec {v}}\cdot {\overrightarrow {\nabla }}\right){\vec {v}}$
où le second membre est nommé advection.
Les forces de viscosité correspondent au terme :
$\nu \nabla ^{2}{\vec {v}}$
Les forces de pression correspondent au terme :
${\frac {1}{\rho }}{\overrightarrow {\nabla }}p$
Les forces extérieures correspondent au terme ${\vec {f}}$ .
L’incompressibilité se traduit par la relation :
${\overrightarrow {\nabla }}\cdot u=0$

Théorème de décomposition de Helmholtz-Hodge

À 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 inclut le terme de projection $\mathbb {P}$ .

\mathbb {P} {\frac {\partial {\vec {v}}}{\partial t}}=\mathbb {P} \left(-\left({\vec {v}}\cdot {\overrightarrow {\nabla }}\right){\vec {v}}-{\frac {1}{\rho }}{\overrightarrow {\nabla }}p+\nu \nabla ^{2}{\vec {v}}+{\vec {f}}\right)

De plus, ce théorème nous permet de supprimer plusieurs éléments de l'équation, comme l'équation de l'incompressibilité et le terme relatif à la pression :

\mathbb {P} ({\overrightarrow {\nabla }}v)=0

\mathbb {P} ({\overrightarrow {\nabla }}p)=0

\mathbb {P} ({\vec {v}})={\vec {v}}

Ainsi on arrive à l'équation suivante :

{\frac {\partial {\vec {v}}}{\partial t}}=\mathbb {P} \left(-\left({\vec {v}}\cdot {\overrightarrow {\nabla }}\right){\vec {v}}+\nu \nabla ^{2}{\vec {v}}+{\vec {f}}\right)

Le modèle numérique

Concernant la résolution numérique, 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 champ de vecteurs d'entrée, effectue des modifications sur le champ et donne le résultat au composant suivant.

On définit des opérateurs pour l'advection ( $\mathbb {A}$ ), la diffusion ( $\mathbb {D}$ ), l'application des forces extérieures ( $\mathbb {F}$ ) et la projection ( $\mathbb {P}$ ) et $\mathbb {S}$ opérateur équivalent à la formulation mathématique pour un pas-de-temps.

\mathbb {S} (u)=\mathbb {P} \circ \mathbb {F} \circ \mathbb {D} \circ \mathbb {A} (u)

Ainsi, à chaque pas-de-temps, le champ 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.

Mise en œuvre de la méthode

Considérons un cas 2D (peu de différence pour la 3D), et considérons une grille régulière de N×N cellules, où chaque cellule est définie 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 $v^{i,j}$ ou i et j correspondent à la position de la cellule courante dans la grille.

Diffusion

L'étape de diffusion correspond à la résolution de la viscosité. La vitesse du fluide va influencer les cellules voisines. Pour cela, un solveur linéaire sera utilisé.

Si $v^{i,j}$ correspond à la vélocité d'une des cellules et $v^{i+1,j}$ , $v^{i-1,j}$ , $v^{i,j+1}$ et $v^{i,j-1}$ aux cellules voisines, alors le solveur doit résoudre le système où pour chaque cellule :

(1-4\mathrm {A} )\cdot {\vec {v^{i,j}}}=\mathrm {A} \cdot {\vec {v^{i+1,j}}}+\mathrm {A} \cdot {\vec {v^{i-1,j}}}+\mathrm {A} \cdot {\vec {v^{i,j+1}}}+\mathrm {A} \cdot {\vec {v^{i,j-1}}}

où

$A = d t \cdotν\cdotN\timesN$
$d t$ : le pas de temps.
$ν$ : coefficient de diffusion/viscosité.
$N\timesN$ : nombre de cellules ⇒ correspond à 1/Aire.

Pour la résolution du système, le solveur du Gauss-Seidel peut être utilisé. Dans cet exemple, on travaille dans un tableau où seuls les voisins directs sont utilisés, les autres cellules n'apparaissent pas et ont comme coefficient 0.

Pseudo-Code, Algorithme diffusion

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[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

Durant l'étape d'advection, on considère donc qu'il y a une particule au niveau du vecteur vélocité, on effectue un retour sur trace (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.

Pseudo-Code, Algorithme advection

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 grille
	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

La 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.

div^{ij}=-0,5(vx^{i+1,j}-vx^{i-1,j})\cdot \mathrm {d} y-0,5(vy^{i,j+1}-vy^{i,j-1})\cdot \mathrm {d} x

vx^{i,j}

: composante horizontale du vecteur vélocité

{\vec {v^{i,j}}}

vy^{i,j}

: composante verticale du vecteur vélocité

{\vec {v^{i,j}}}

\mathrm {d} x

: distance entre les vecteurs

v^{i+1,j}

et

v^{i-1,j}

\mathrm {d} y

: distance entre les vecteurs

v^{i,j+1}

et

v^{i,j-1}

Dans le cas d'une grille de N×N cellules on a $\mathrm {d} x=\mathrm {d} y={\frac {1}{\mathrm {N} }}$

Pseudo-Code, Algorithme projection

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

Extensions de la méthode

Staggered-Grid

Dans 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.

Retour sur trace

Le retour sur trace (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 retour sur trace simple : on peut remarquer que le retour sur trace 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

Dans 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

La 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

Dans les cas précédents, ce sont des cellules carrées (ou au moins des quadrilatères) qui ont été utilisées pour les cas 2D et des cubes pour le cas 3D. Il est néanmoins possible d'utiliser des cellules triangulaires^[10] ou des tétraèdres.

Vorticité

Il 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

Notes

[1] "Stable Fluid" - Jos Stam - 1999

[2] Real-Time Fluid Dynamics for Games

[3] "Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface" - Harlow, F. H. and Welch, J. E. - 1965,

[4] "Numerical Heat Transfer" - SV Patankar

[5] "An Unconditionally Stable MacCormack Method"

[6] "Using BFECC for Fluid Simulation"

[7] "Flows on Surfaces of Arbitrary Topology"

[8] "level-set method for fluid interface"

[9] "Fast and Robust Tracking of Fluid Surfaces"

[10] "Inviscid and Incompressible Fluid Simulation on Triangle Meshes"

[11] "Stable, Circulation-Preserving, Simplicial Fluids"

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]