Stable-Fluids

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

• Équation d'incompressibilité

  ${\displaystyle {\overrightarrow {\nabla }}\cdot {\vec {v}}=0}$

• Équation de bilan de la quantité de mouvement

  ${\displaystyle {\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 ${\displaystyle {\vec {v}}}$ représente le vecteur vitesse à un point donné du fluide.
• Les termes d'accélérations correspondent à :
  ${\displaystyle {\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 :
  ${\displaystyle \nu \nabla ^{2}{\vec {v}}}$
• Les forces de pression correspondent au terme :
  ${\displaystyle {\frac {1}{\rho }}{\overrightarrow {\nabla }}p}$
• Les forces extérieures correspondent au terme ${\displaystyle {\vec {f}}}$.
• L’incompressibilité se traduit par la relation :
  ${\displaystyle {\overrightarrow {\nabla }}\cdot u=0}$

À 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 ${\displaystyle \mathbb {P} }$.

${\displaystyle \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 :

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

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

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

Ainsi on arrive à l'équation suivante :

${\displaystyle {\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)}$

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 (${\displaystyle \mathbb {A} }$), la diffusion (${\displaystyle \mathbb {D} }$), l'application des forces extérieures (${\displaystyle \mathbb {F} }$) et la projection (${\displaystyle \mathbb {P} }$) et ${\displaystyle \mathbb {S} }$ opérateur équivalent à la formulation mathématique pour un pas-de-temps.

${\displaystyle \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.

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 ${\displaystyle v^{i,j}}$ correspond à la vélocité d'une des cellules et ${\displaystyle v^{i+1,j}}$ , ${\displaystyle v^{i-1,j}}$ , ${\displaystyle v^{i,j+1}}$ et ${\displaystyle v^{i,j-1}}$ aux cellules voisines, alors le solveur doit résoudre le système où pour chaque cellule :

${\displaystyle (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 = dt⋅ν⋅N×N
• dt : le pas de temps.
• ν : coefficient de diffusion/viscosité.
• N×N : 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.

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 );
}

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.

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);
}

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.

${\displaystyle 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}$

${\displaystyle vx^{i,j}}$ : composante horizontale du vecteur vélocité ${\displaystyle {\vec {v^{i,j}}}}$

${\displaystyle vy^{i,j}}$ : composante verticale du vecteur vélocité ${\displaystyle {\vec {v^{i,j}}}}$

${\displaystyle \mathrm {d} x}$ : distance entre les vecteurs ${\displaystyle v^{i+1,j}}$ et ${\displaystyle v^{i-1,j}}$

${\displaystyle \mathrm {d} y}$ : distance entre les vecteurs ${\displaystyle v^{i,j+1}}$ et ${\displaystyle v^{i,j-1}}$

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

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 );
}

simple grille et 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 simple et deux étapes du MacCormack

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

Distorsion appliquée aux cellules afin de gérer des flux sur une surface arbitraire.

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

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

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.

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.