Stabilité du schéma MAC

On ne peut pas analyser la stabilité d'un schéma pour une équation non linéaire, il est nécessaire d'en faire une approximation. On linéarise l'équation sous forme non-conservative ( $u u_x + v u_y = u_0 u_x + v_0 u_y$) avec $u_0$ et $v_0$ des constantes, on élimine le décalage des cellules et on oublie la pression :

$$\underbrace{\frac{u_{i,j}^{n+1} - u_{i,j}^{n}}{\Delta \tau}}... ...frac{1}{Re} (\delta_x^2 u_{i,j}^{n} + \delta_y^2 u_{i,j}^{n} )$$

Tout comme dans le chapitre $5$, $\delta_x$ représente les différences centrées. On suppose l'erreur de la forme :

$$\epsilon_{i,j}^n = G^n e^{ipx + iqy}$$

alors :

$$\begin{array}{rcl} \delta_x \epsilon_{i,j}^n&=&\frac{\epsilo... ...big( -4 \sin^2(\frac{p \Delta x}{2}) \big)\\ &&\\ \end{array}$$

et :

$$\frac{\epsilon_{i,j}^{n+1}-\epsilon_{i,j}^{n}}{\Delta \tau} = \frac{G-1}{\Delta \tau} G^n e^{i p x_i + i q y_j}$$

Finalement, avec l'hypothèse $\Delta x = \Delta y$ et en prenant $\xi = p \Delta x$ et $\eta = q \Delta y$, on a :

$$\frac{G-1}{\Delta \tau} + \underbrace{\frac{u_0}{ \Delta x}... ...\frac{\xi}{2}+\sin^2 \frac{\eta}{2} \right)}_{\textrm{R\'eel}}$$

donc :

$$\vert G\vert^2 = \left\vert 1- \frac{4}{Re \Delta x^2} \Del... ...ght)^2 \left\vert u_0 \sin{\xi} + v_0 \sin{\eta} \right\vert^2$$

Il y a évidement une limite sur la longueur d'onde que l'on peut capturer numériquement avec un maillage donné ; il faut au moins deux intervales $\Delta x$ pour capturer une onde (th. de Shannon).

Soit $k$ le nombre d'onde, on a alors comme longueur d'onde $\frac{2 \pi}{k}$ et :

$$2 \Delta x \leqslant \frac{2 \pi }{k} \leqslant L$$

avec $L$ la longueur du domaine selon $x$. Ainsi :

$$k \Delta x = \xi \leqslant \pi$$

On choisi donc la valeur maximale de $\xi$ et $\eta$ comme :

$$\xi = \eta = \pi$$

dans ce cas :

$$\vert G\vert^2 = \left\vert 1 - \frac{4}{Re} \frac{\Delta \tau}{\Delta x^2}.2 \right\vert^2$$

Le schéma est alors stable si :

$$\begin{array}{c} -1 \leqslant 1 -\frac{8}{Re} \frac{\Delta \... ...c{8}{Re} \frac{\Delta \tau}{\Delta x^2} \leqslant 0 \end{array}$$

soit :

$$\frac{4}{Re} \frac{\Delta \tau}{\Delta x^2} \leqslant 1$$

Si on avait choisi pour $\xi$ et $\eta$ :

$$\xi = \eta = \frac{\pi}{2}$$

on aurait eu :

$$\vert G\vert^2 = \left\vert 1 - \frac{4}{Re} \frac{\Delta \t... ...lta \tau}{\Delta x}\right)^2 \vert u_0+v_0\vert^2 \leqslant 1$$

pour la stabilité, soit :

$$1 - \frac{4}{Re} \frac{\Delta \tau}{\Delta x^2} \left( 2 - \... ...a \tau}{\Delta x} \right)^2 \vert u_0 + v_0\vert^2 \leqslant 1$$

En majorant le terme à gauche et en simplifiant, on a :

$$2-\frac{4}{Re} \frac{\Delta \tau}{\Delta x^2} + \big( \frac{\Delta \tau}{\Delta x} \big)^2 \vert u_0+v_0\vert^2 \geqslant 0$$

donc :

$$\frac{1}{4} Re \Delta \tau \vert u_0 + v_0\vert^2 \leqslant 1$$

Cette condition donne la limite de stabilité convective.

Les deux limites trouvées doivent être satisfaites. Une troisième analyse de stabilité peut être effectuée en éliminant le terme non-linéaire et en cherchant à calculer la stabilité de l'écoulement de Stokes. Dans ce cas, on peut trouver :

$$4 \frac{\Delta \tau}{\Delta x^2} \left( \frac{1}{Re} + \frac{\Delta \tau c^2}{2} \right) \leqslant 1$$

et on trouve une méthode sur $c^2$.

Une petite remarque avant de terminer cette section. Dans la section suivante, on trouve la méthode fonction de courant, vorticité qui peut s'appliquer pour des écoulements stationnaires avec une méthode similaire à celle de la compressiblité artificielle.