4.4 Équation de convection-diffusion

4.4.1 Problème physique: dispersion d'un polluant

On considère le problème de la dispersion d'un polluant à la surface d'un liquide en mouvement (figure 4.15). Le polluant est miscible dans le liquide, mais plus léger. On suppose que la vitesse du fluide est horizontale. On peut alors considérer que le polluant diffuse essentiellement à la surface, et négliger les variations suivant la vertical $z$ .

Figure 4.15: diffusion d'une tache de polluant

En notant $C$ la fraction massique de polluant et $\rho$ la densité du fluide, l'équation d'équilibre pour $C$ traduit que la variation temporelle de la quantité de polluant $\rho Cdxdydz$ dans un volume élémentaire $dxdydz$ est égale à un bilan de flux de $C$ sur les facettes du volume. A travers une facette de surface $dS$ et de normale sortante $\vec{n}$ , il existe un flux de diffusion $\lambda\overrightarrow{\nabla}.\vec{n}dS$ ( $\lambda$ est le coefficient de diffusion ) et un flux de convection $-\rho C\overrightarrow{V}.\vec{n} dS$ ( $\overrightarrow{V}$ est la vitesse du fluide). L'équation d'équilibre s'écrit:

$$\frac{\partial\rho C}{\partial t}=div(\lambda\overrightarrow{\nabla}C)-div(\rho C\overrightarrow{V})$$

Compte tenu de l'équation de conservation de la masse du fluide:

$$\frac{\partial\rho}{\partial t}+div(\rho\overrightarrow{V})=0$$

cette équation s'écrit:

$$\rho\frac{\partial C}{\partial t}+\rho\overrightarrow{V}.\overrightarrow{\nabla}C=div(\lambda\overrightarrow{\nabla}C)$$

Compte tenu de l'hypothèse d'indépendance des quantités par rapport à $z$ , et en supposant en outre que la densité $\rho$ et le coefficient de diffusion $\lambda$ sont constants, il vient:

$$\frac{\partial C}{\partial t}+V_{1}\frac{\par... ...tial^{2}C}{\partial x^{2}}+\frac{\partial^{2}C}{\partial y^{2}}\right)=0$$

A cette équation il faut ajouter la condition initiale $C(x,y,t=0)=w(x,y)$ et des conditions aux limites sur les frontières $\Gamma$ du domaine $\Omega$ . On suppose que le polluant se trouve initialement à l'intérieur du domaine $\Omega$ et a pour dimension caractéristique $\delta$ (figure 4.15). On distingue les frontières entrantes $\Gamma_{0}$ , i.e. telle que $\overrightarrow{V}.\overrightarrow{n}<0$ et les frontières sortantes $\Gamma_{1}$ , i.e. telle que $\overrightarrow{V}.\overrightarrow{n}>0$ . Si la convection est prépondérante sur la diffusion, le flux entrant sur les frontières $\Gamma_{0}$ est un flux de convection qui apporte du fluide non pollué dans le domaine $\Omega$ . La condition sur $\Gamma_{0}$ est une condition de Dirichlet $C_{\Gamma_{0}}=0$ . Sur les frontières $\Gamma_{1}$ , le fluide transporte le polluant vers l'extérieur: on impose alors une condition aux limites de Neumann $\frac{\partial C}{\partial n}_{\Gamma_{1}}=0$ , qui autorise la sortie du polluant du domaine $\Omega$ .

Pour un domaine carré de dimension $H$ , et une vitesse $\overrightarrow{V}$ avec des composantes positives $V_{1}>0$ et $V_{2}>0$ , la frontière $\Gamma_{0}$ correspond aux deux cotés $x=0$ et $y=0$ et la frontière $\Gamma_{1}$ aux deux autres cotés opposés $x=H$ et $y=H$ .

En notant $\kappa=\frac{\lambda}{\rho}$ , le problème modèle s'écrit pour un domaine $\Omega$ carré de dimension caractéristique $H$ :

Trouver $u(x,y,t)$ tel que:

$$\frac{\partial u}{\partial t}+V_{1}\frac{\partial u}{\partial x}+... ...partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right)=0 \Omega=[0,H]*[0,H]$$ (4.26)

$$u_{\Gamma_{0}}=0 \Gamma_{0}(x=0, y=0) \„ (\frac{\partial u}{\partial n})_{\Gamma_{1}}=0 \Gamma_{1}(x=H, y=H)$$

$$u(x,y,t=0)=w(x,y)\textrm{}$$

4.4.2 Étude de la solution exacte

Le problème (4.26) est un problème d'évolution parabolique caractéristique des problèmes de mécanique des fluides, avec un terme de convection et un terme de diffusion. Nous allons tout d'abord effectuer une analyse d'ordre de grandeur de chacun de ces termes.

4.4.2.1 Analyse en ordre de grandeur

Soit $\delta$ la dimension caractéristique de la tache initiale, l'ordre de grandeur des différents termes de l'équation (4.26) s'écrit, en notant $\tau$ un temps caractéristique:

$$\frac{\Delta C}{\tau}+V_{1}\frac{\Delta C}{\d... ...da\left(\frac{\Delta C}{\delta^{2}}+\frac{\Delta C}{\delta^{2}}\right)=0$$

Si on considère uniquement la diffusion, le temps caractéristique vaut:

$$\tau_{d}\approx\frac{\delta^{2}}{\lambda}$$

On retrouve le temps caractéristique de diffusion (relation c1tdiff) calculé pour l'équation de la chaleur. Le temps caractéristique $\tau_{d}$ vérifie une relation du type (à rapprocher de la condition de stabilité c1eq36 pour une équation de diffusion):

$$\frac{\lambda\tau_{d}}{\delta^{2}}\approx1$$ (4.27)

C'est le temps caractéristique de décroissance exponentielle des phénomènes de diffusion.

Si on considère uniquement la convection, le temps caractéristique vaut:

$$\tau_{c}\approx\frac{\delta}{V}$$

Ce temps correspond au temps de transport de la tache par le champ de vitesse sur une distance égale à la dimension de la tache. Ce temps caractéristique $\tau_{c}$ vérifie donc une relation du type (à rapprocher de la condition de stabilité de Courant c3eq12):

$$\frac{V \tau_{c}}{\delta}\approx1$$ (4.28)

Pour notre problème, on peut définir un autre temps caractéristique de convection: le temps de sortie $\tau_{s}$ du polluant hors du domaine $\Omega$ . Il est basé sur la dimension caractéristique $L$ du domaine et vérifie:

$$\frac{V \tau_{s}}{L}\approx1$$ (4.29)

Enfin le rapport entre les temps caractéristiques de diffusion $\tau_{d}$ et de convection $\tau_{c}$ est le nombre de Péclet:

$$\frac{\tau_{d}}{\tau_{c}}=\frac{V\delta}{\lambda}=P_{e}$$ (4.30)

qui caractérise l'importance relative du terme de convection par rapport au terme de diffusion.

4.4.2.2 Diffusion convection d'une gaussienne

Considérons la condition initiale suivante:

$$w(x,y)=e^{-\left(\frac{x-x_{0}}{\sigma_{0}}\right)^{2}}e^{-\left(\frac{y-y_{0}}{\sigma_{0}}\right)^{2}}$$ (4.31)

qui décrit une tache gaussienne d'amplitude 1 centrée en $(x_{0},y_{0})$ , de rayon $\sigma_{0}$ .

S'il n'y a pas de convection, cette tache diffuse de façon auto-similaire , i.e. son amplitude diminue et son rayon augmente en conservant une forme gaussienne:

$$u(x,y,t)=A(t)e^{-\left(\frac{x-x_{0}}{\sigma(t)}\right)^{2}}e^{-\left(\frac{y-y_{0}}{\sigma(r)}\right)^{2}}$$

En utilisant la conservation globale de $u$ dans tout le domaine:

$$\int\int u(x,y,t) dxdy=\int\int w(x,y) dxdy=cste \forall t$$

on en déduit la relation entre l'amplitude $A(t)$ et le rayon $\sigma(t)$ :

$$\frac{\pi}{4}A(t)\sigma^{2}(t)=\frac{\pi}{4}\sigma_{0}^{2}$$

(l'intégrale d'une gaussienne vaut: $\int e^{-(x/\sigma)^{2}}dx=\frac{\sqrt{\pi}}{2}\sigma$ ).

On cherche donc une solution de diffusion de l'équation (4.26) (avec $V=0$ ) sous la forme:

$$u(x,y,t)=\left(\frac{\sigma_{0}}{\sigma(t)}\r... ...}}{\sigma(t)}\right)^{2}}e^{-\left(\frac{y-y_{0}}{\sigma(r)}\right)^{2}}$$

En reportant cette relation dans l'équation (4.26) on obtiens l'équation d'évolution de $\sigma(t)$ :

$$\frac{d}{dt}\sigma(t)=-2\frac{\kappa}{\sigma(t)}$$

dont la solution vérifiant $\sigma(0)=\sigma_{0}$ est:

$$\sigma(t)=\sqrt{4\kappa t+\sigma_{0}^{2}}$$ (4.32)

La solution de diffusion de l'équation (4.26) s'écrit donc:

$$u(x,y,t)=\left(\frac{\sigma_{0}}{\sqrt{4\kappa t+\sigma_{0}^{2}}}... ...ight)^{2}}e^{-\left(\frac{y-y_{0}}{\sqrt{4\kappa t+\sigma_{0}^{2}}}\right)^{2}}$$ (4.33)

L'amplitude de cette gaussienne décroît donc suivant la loi:

$$A(t)=\left(\frac{\sigma_{0}}{\sqrt{4\kappa t+\sigma_{0}^{2}}}\right)^{2}$$ (4.34)

En prenant en compte la convection par un champ de vitesse sans cisaillement, cette tache gaussienne est transportée sans déformation et diffuse le long des trajectoires du champ de vitesse comme précédemment. Pour un champ de vitesse constant, les trajectoires sont des droites:

$$x(t)=x_{0}+V_{1}t \„ y(t)=y_{0}+V_{2}t$$

la solution de convection-diffusion de l'équation (4.26) s'écrit donc:

$$u(x,y,t)=\left(\frac{\sigma_{0}}{\sqrt{4\kappa t+\sigma_{0}^{2}}}... ...2}}e^{-\left(\frac{y-y_{0}-V_{2}t}{\sqrt{4\kappa t+\sigma_{0}^{2}}}\right)^{2}}$$ (4.35)

Cette solution est une solution en milieu infini et ne tiens pas compte des conditions aux limites du problème (4.26). Elle constitue cependant une bonne approximation de la solution, si la dimension de la tache $\sigma_{0}$ est petite devant la dimension $L$ du domaine $\Omega$ .

Figure: convection diffusion d'une gaussienne (4.35)

[solution a t=0 et t=0.5]

[amplitude $A(t)$ et taille $\sigma (t)$ ]

On a tracé l'évolution de cette solution (4.35) sur la figure (4.16) pour $\kappa=0.01$ , $\sigma_{0}=0.1$ et $\overrightarrow{V}=[1.\„ 1.]$ . Pour ces valeurs des paramètres, un domaine de longueur $H=1$ , une position initiale $x_{0}=\frac{H}{4}$ et $y_{0}=\frac{H}{4}$ , le temps caractéristique de diffusion vaut $\tau_{d}\approx1.$ , celui de convection $\tau_{c}\approx0.1$ , et le nombre de Péclet $P_{e}\approx10$ . Le problème est donc à convection dominante. Au bout d'un temps $T_{f}=0.5$ $(T_{f}<\tau_{s}=1.0)$ , la tache est à la frontière du domaine avec une amplitude qui a diminuée d'un tiers: $A(T_{f})\approx0.33333$ .

4.4.2.3 Modes propres de diffusion

Pour rechercher des solutions vérifiant les conditions aux limites, on détermine tout d'abord les modes propres de diffusion en utilisant la méthode de séparation de variable décrite au paragraphe c1analytique. Le calcul est identique, et on montre facilement que les modes propres sont les fonctions suivantes:

$$u_{p,q}(x,y,t)=\sin\left(\frac{(2p+1)\pi}{2}\frac{x}{L}\right)\si... ...(2p+1)\pi}{2L}\right)^{2}+\left(\frac{(2q+1)\pi}{2L}\right)^{2}\right)\kappa t}$$ (4.36)

qui vérifient les conditions aux limites:

$$u_{p,q}(0,y,t)=u_{p,q}(x,0,t)=0 (\Gamm... ...frac{d}{dx}u_{p,q}(H,y,t)=\frac{d}{dy}u_{p,q}(x,H,t)=0 (\Gamma_{1})$$

La solution générale de diffusion est alors une combinaison linéaire de ses modes:

$$u(x,y,t)=\sum_{p=0}^{\infty}\sum_{q=0}^{\infty}\alpha_{p,q}u_{p,q}(x,y,t)$$

Figure 4.17: mode propre de diffusion $p=1$ $q=1$

[t=0]

[Amplitude $A(t)$ ]

On a tracé sur la figure (4.17) le mode $p=1$ et $q=1$ , ainsi que l'évolution temporelle de son amplitude pour les mêmes paramètres que précédemment ( $\kappa=0.01$ , $L=1$ ). Sur un temps $T_{f}=0.5$ , l'amplitude de ce mode décroît de $0,8$ .

Si on prend en compte la convection dans le cas d'un champ de vitesse sans cisaillement, la solution initiale est convectée sans déformation et diffusée le long des trajectoires. Par contre, il n'existe pas de solutions analytiques simples qui vérifient les conditions aux limites de (4.26).

En considérant une taille de structure $\delta=\frac{2}{3\pi}L$ , le temps caractéristique de diffusion vaut $\tau_{d}\approx4.5$ , le temps caractéristique de convection $\tau_{c}\approx0,2$ , et le nombre de Péclet $P_{e}\approx21$ .

4.4.3 Discrétisation par différences finies

Nous avons vu dans les chapitres précédents qu'une discrétisation précise de problème parabolique est le schéma de Cranck Nicholson. Appliquée à l'équation (4.26), il s'écrit pour un maillage cartésien de $N_{x}$ points suivant $x$ et $N_{y}$ points suivant $y$ , et de pas $dx$ et $dy$ :

$$\frac{u_{i,j}^{n+1}-u_{i,j}^{n}}{dt} = -\frac{V_{1}}{2}\left(\frac{u_{i+1,j}^{n+1}-u_{i-1,j}^{n+1}}{2dx}+\frac{u_{i+1,j}^{n}-u_{i-1,j}^{n}}{2dx}\right)$$ (4.37)

$$- \frac{V_{2}}{2}\left(\frac{u_{i,j+1}^{n+1}-u_{i,j-1}^{n+1}}{2dy}+\frac{u_{i,j+1}^{n}-u_{i,j-1}^{n}}{2dy}\right)$$

$$+ \frac{\kappa}{2}\left(\frac{u_{i+1,j}^{n+1}-2u_{i,j}^{n+1}+u_{i-1,j}^{n+1}}{dx^{2}}+\frac{u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}}{dx^{2}}\right)$$

$$+ \frac{\kappa}{2}\left(\frac{u_{i,j+1}^{n+1}-2u_{i,j}^{n+1}+u_{i,j-1}^{n+1}}{dy^{2}}+\frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{dy^{2}}\right)$$

C'est un schéma inconditionnellement stable d'ordre 2 en temps et en espace, i.e. en $O(dt^{2},dx^{2},dy^{2})$ .

A chaque itération en temps, on a à résoudre un système d'équations linéaire $\mathcal{A}x=\mathcal{B}$ , de $N=N_{x}N_{y}$ inconnues $\{u_{i,j}^{n+1}\}$ . La matrice $\mathcal{A}$ est une matrice penta-diagonale, qui a la même structure que la matrice du laplacien au paragraphe c4matlap. Pour des très gros maillages, le coût de résolution de ce système linéaire en utilisant les méthodes de résolution du paragraphe c4matlap peut devenir rapidement prohibitif.

On va donc étudier dans le paragraphe suivant une méthode alternative: la méthode des directions alternées implicites.

4.4.4 Méthode des directions alternées implicites

Le principe des méthodes des directions alternées implicites, notées ADI (ADI=Alternated Directions Implicited est le sigle classique des directions alternées en anglais), est de décomposer les opérateurs spatiaux suivant les directions d'espace $x$ et $y$ . On écrit l'équation (4.26) sous la forme symbolique suivante:

$$\frac{\partial u}{\partial t}=L_{1}(u)+L_{2}(u)$$ (4.38)

où $L_{1}$ et $L_{2}$ sont les deux opérateurs suivants:

$$L_{1}=V_{1}\frac{\partial}{\partial x}-\kappa\frac{\partial^{2}}{... ..._{2}=V_{2}\frac{\partial}{\partial y}-\kappa\frac{\partial^{2}}{\partial y^{2}}$$ (4.39)

En notant $u^{n}$ , $u^{n+1}$ , $u^{n+\frac{1}{2}}$ les solutions au temps $t^{n}=ndt$ , $t^{n+1}=(n+1)dt$ et $t^{n+\frac{1}{2}}=(n+\frac{1}{2})dt$ , les développements limitées de $u^{n+\frac{1}{2}}$ peuvent s'écrire de façon symbolique:

$$u^{n+\frac{1}{2}} = \left(e^{\frac{dt}{2}\frac{\partial}{\partial t}}\right)u^{n}=\left(e^{\frac{dt}{2}(L_{1}+L_{2})}\right)u^{n}$$

$$u^{n+\frac{1}{2}} = \left(e^{-\frac{dt}{2}\frac{\partial}{\partial t}}\right)u^{n+1}=\left(e^{-\frac{dt}{2}(L_{1}+L_{2})}\right)u^{n+1}$$

On a utilisé dans ces relations le fait que $u^{n}$ et $u^{n+1}$ sont solutions de l'équation exacte (4.26) pour remplacer $\frac{\partial u}{\partial t}$ en fonction de $L_{1}(u)$ et $L_{2}(u)$ .

En combinant ces deux équations, il vient:

$$\left(e^{-\frac{dt}{2}(L_{1}+L_{2})}\right)u^{n+1}=\left(e^{\frac{dt}{2}(L_{1}+L_{2})}\right)u^{n}$$ (4.40)

En utilisant un développement au premier ordre des exponentielles,

$$\left(1-\frac{dt}{2}(L_{1}+L_{2})\right)u^{n+1}=\left(1+\frac{dt}{2}(L_{1}+L_{2})\right)u^{n}$$

on obtiens le schéma suivant, dans lequel il suffit d'inclure l'approximation spatiale des opérateurs $L_{1}$ et $L_{2}$ :

$$\frac{u^{n+1}-u^{n}}{dt}=\frac{L_{1}(u^{n+1})+L_{1}(u^{n})}{2}+\frac{L_{2}(u^{n+1})+L_{2}(u^{n})}{2}$$ (4.41)

C'est le schéam classique de Cranck Nicholson (4.37).

Pour les schémas ADI, on effectue tout d'abord une factorisation formelle dans (4.40):

$$\left(e^{-\frac{dt}{2}L_{1}}e^{-\frac{dt}{2}L... ...ght)u^{n+1}=\left(e^{\frac{dt}{2}L_{1}}e^{\frac{dt}{2}L_{2}}\right)u^{n}$$

avant le développement limité des exponentielles:

$$\left(1-\frac{dt}{2}L_{1}\right)\left(1-\frac{dt}{2}L_{2}\right)u^{n+1}=\left(1+\frac{dt}{2}L_{1}\right)\left(1+\frac{dt}{2}L_{2}\right)u^{n}$$ (4.42)

Formellement, on a un schéma de type Cranck Nicholson d'ordre 2 mais avec une erreur de troncature différente. Pour résoudre, on introduit la solution intermédiaire $u^{*}$ telle que:

$$\left(1-\frac{dt}{2}L_{1}\right)u^{*} = \left(1+\frac{dt}{2}L_{2}\right)u^{n}$$ (4.43)

$$\left(1-\frac{dt}{2}L_{2}\right)u^{n+1} = \left(1+\frac{dt}{2}L_{1}\right)u^{*}$$

Ces deux équations sont équivalentes à l'équation (4.42). Pour s'en convaincre, il suffit de multiplier la première par $1+\frac{dt}{2}L_{1}$ et la seconde par $1-\frac{dt}{2}L_{1}$ et de les combiner. L'intérêt de cette procédure par rapport à Cranck Nicholson classique est que dans la première équation (4.43), on est implicite suivant $L_{1}$ (i.e. suivant $x$ ) et explicite suivant $L_{2}$ (i.e. suivant $y$ ), et vice-versa dans la seconde. La résolution de ces deux équations sera donc plus facile que la résolution du schéma de Cranck Nicholson (4.41), dans lequel on est implicite suivant les 2 directions $L_{1}$ et $L_{2}$ .

On discrétise ensuite les opérateurs $L_{1}$ et $L_{2}$ avec des différences finies centrées sur un maillage cartésien de $N_{x}$ points suivant $x$ et $N_{y}$ points suivant $y$ , et de pas $dx$ et $dy$ :

$$L_{1}(u_{i,j}) = -V_{1}\frac{u_{i+1,j}-u_{i-1,j}}{2dx}+\kappa\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{dx^{2}}$$

$$L_{2}(u_{i,j}) = -V_{2}\frac{u_{i,j+1}-u_{i,j-1}}{2dy}+\kappa\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{dy^{2}}$$

On obtiens le schéma ADI suivant pour l'équation (4.26):

$$\frac{u_{i,j}^{*}-u_{i,j}^{n}}{dt/2} = -V_{1}\frac{u_{i+1,j}^{*}-u_{i-1,j}^{*}}{2dx}+\kappa\frac{u_{i+1,j}^{*}-2u_{i,j}^{*}+u_{i-1,j}^{*}}{dx^{2}}$$ (4.44)

$$-V_{2}\frac{u_{i,j+1}^{n}-u_{i,j-1}^{n}}{2dy}+\kappa\frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{dy^{2}}$$

$$\frac{u_{i,j}^{n+1}-u_{i,j}^{*}}{dt/2} = -V_{1}\frac{u_{i+1,j}^{*}-u_{i-1,j}^{*}}{2dx}+\kappa\frac{u_{i+1,j}^{*}-2u_{i,j}^{*}+u_{i-1,j}^{*}}{dx^{2}}$$ (4.45)

$$-V_{2}\frac{u_{i,j+1}^{n+1}-u_{i,j-1}^{n+1}}{2dy}+\kappa\frac{u_{i,j+1}^{n+1}-2u_{i,j}^{n+1}+u_{i,j-1}^{n+1}}{dy^{2}}$$

Figure 4.18: schéma ADI4.44 et 4.45

Ces équations ADI correspondent à deux discrétisations de l'équation (4.26) à $t^{n+\frac{1}{2}}$ avec un pas en temps $\frac{dt}{2}$ . La première correspond à une discrétisation entre $t^{n}$ et $t^{n+\frac{1}{2}}$ avec une implicitation dans la direction $x$ et la seconde à une discrétisation entre $t^{n+\frac{1}{2}}$ et $t^{n}$ avec une implicitation dans la direction $y$ (figure 4.18). La variable intermédiaire $u_{i,j}^{*}$ corresponds donc à une approximation de $u_{i,j}^{n+\frac{1}{2}}$ .

La première équation (4.44) couple les valeurs inconnues $u_{i,j}^{*}$ par ligne (figure 4.18), i.e. les valeurs inconnues sur la ligne $j$ $\{u_{i,j}^{*}\}_{i=1,N_{x}}$ sont indépendantes des valeurs sur les autres lignes. Elles sont solutions du système linéaire tri-diagonal $\mathcal{A}^{1}x=\mathcal{B}^{1}$ de dimension $N_{x}$ suivant:

$$\left[\begin{array}{ccccc} \ddots & \ddots & 0 & \cdots & 0 \dd... ...c_{2}^{1}u_{i,j}^{n}+c_{3}^{1}u_{i,j+1}^{n} \vdots \vdots\end{array}\right]$$ (4.46)

avec $a_{1}^{1}=\frac{-\kappa}{dx^{2}}-\frac{V_{1}}{2dx}$ , $a_{2}^{1}=\frac{2}{dt}+\frac{2\kappa}{dx^{2}}$ , $a_{3}^{1}=\frac{-\kappa}{dx^{2}}+\frac{V_{1}}{2dx}$ , et $c_{1}^{1}=\frac{\kappa}{dy^{2}}+\frac{V_{2}}{2dy}$ , $c_{2}^{1}=\frac{2}{dt}-\frac{2\kappa}{dy^{2}}$ , $c_{3}^{1}=\frac{\kappa}{dy^{2}}-\frac{V_{2}}{2dy}$ .

Le second membre $\mathcal{B}^{1}$ peut s'écrire sous la forme d'un produit matrice vecteur $\mathcal{B}^{1}=\mathcal{C}^{1}\{u_{i,j}^{n}\}$ avec une matrice $\mathcal{C}^{1}$ tridiagonale:

$$\mathcal{B}^{1}=\left[\begin{array}{ccccc} \d... ..._{i,j-1}^{n}\\ u_{i,j}^{n}\\ u_{i,j+1}^{n}\\ \vdots\end{array}\right]$$

Pour déterminer les valeurs inconnues $u_{i,j}^{*}$ , il faut donc résoudre $N_{y}$ systèmes linéaires tri-diagonaux de dimension $N_{x}$ .

De même la seconde équation (4.45) couple les valeurs inconnues $u_{i,j}^{n+1}$ par colonne (figure 4.18), i.e. les valeurs inconnues sur la colonne $i$ $\{u_{i,j}^{n+1}\}_{j=1,N_{y}}$ sont indépendantes des valeurs sur les autres colonnes. Elles sont solutions du système linéaire tri-diagonal $\mathcal{A}^{2}x=\mathcal{B}^{2}$ de dimension $N_{y}$ suivant:

$$\left[\begin{array}{ccccc} \ddots & \ddots & 0 & \cdots & 0 \dd... ...c_{2}^{2}u_{i,j}^{*}+c_{3}^{2}u_{i+1,j}^{*} \vdots \vdots\end{array}\right]$$ (4.47)

avec $a_{1}^{2}=\frac{-\kappa}{dy}-\frac{V_{2}}{2dy}$ , $a_{2}^{2}=\frac{2}{dt}+\frac{2\kappa}{dy^{2}}$ , $a_{3}^{2}=\frac{-\kappa}{dy^{2}}+\frac{V_{2}}{2dy}$ , et $c_{1}^{1}=\frac{\kappa}{dx^{2}}+\frac{V_{1}}{2dx}$ , $c_{2}^{1}=\frac{2}{dt}-\frac{2\kappa}{dx^{2}}$ , $c_{3}^{1}=\frac{\kappa}{dx^{2}}-\frac{V_{1}}{2dx}$ .

Le second membre $\mathcal{B}^{2}$ peut encore s'écrire sous la forme d'un produit matrice vecteur $\mathcal{B}^{2}=\mathcal{C}^{2}\{u_{i,j}^{*}\}$ :

$$\mathcal{B}^{2}=\left[\begin{array}{ccccc} \d... ..._{i-1,j}^{*}\\ u_{i,j}^{*}\\ u_{i+1,j}^{*}\\ \vdots\end{array}\right]$$

Pour déterminer valeurs inconnues $u_{i,j}^{n+1}$ , il faut donc résoudre $N_{x}$ systèmes linéaires tri-diagonaux de dimension $N_{y}$ .

A chaque itération en temps, on résoud $N_{y}$ systèmes tri-diagonaux de rang $N_{x}$ et $N_{x}$ systèmes tri-diagonaux de rang $N_{y}$ , ce qui est beaucoup plus efficace que la résolution d'un seul système linéaire de rang $N_{x}N_{y}$ .

Pour les conditions aux limites, il faut modifier la première et la dernière ligne de ces systèmes linéaires. Pour les conditions de Dirichlet sur $\Gamma_{0}$ , on modifie la première ligne de $\mathcal{A}^{1}$ et $\mathcal{A}^{2}$ :

$$\mathcal{A}_{1,j\neq1}^{1}=0 \„ \mathcal{A}_{1,1}^{1}=1\;, \mathcal{B}_{1}^{1}=0 \mathcal{A}_{1,j\neq1}^{2}=0 \„ \mathcal{A}_{1,1}^{2}=1\;, \mathcal{B}_{1}^{2}=0$$ (4.48)

Cette condition fixe la valeur de la première ligne de $u^{*}$ : $\{u_{i,1}^{*}\}_{i=1,N_{x}}=0$ et la première colonne de $u^{n+1}$ : $\{u_{^{1,j}}^{n+1}\}_{j=1,N_{y}}=0$ .

Pour les conditions de Neumann sur $\Gamma_{1}$ , on utilise une condition miroir, qui modifie la dernière ligne de $\mathcal{A}^{1}$ et $\mathcal{A}^{2}$ :

$$\mathcal{A}_{N_{x},N_{x-1}}^{1}=a_{1}^{1}+a_{3}^{1} \mathcal{A}_{N_{y},N_{y-1}}^{2}=a_{1}^{2}+a_{3}^{2}$$ (4.49)

De même le second membre $\mathcal{B}^{1}$ de l'étape 1 pour la ligne $j=N_{y}$ est modifié:

$$\mathcal{B}_{i}^{1}=(c_{1}^{1}+c_{3}^{1})u_{i,N_{y}-1}^{n}+c_{2}^{1}u_{i,N_{y}}^{n}$$ (4.50)

ainsi que le second membre $\mathcal{B}^{2}$ de l'étape 2 pour la colonne $i=N_{x}$ :

$$\mathcal{B}_{j}^{2}=(c_{1}^{2}+c_{3}^{2})u_{N_{x}-1,j}^{*}+c_{2}^{2}u_{N_{x},j}^{*}$$ (4.51)

4.4.4.1 Stabilité et précision du schéma ADI

4.4.4.1.1 Étude de la stabilité:

L'étude de la stabilité utilise le programme Maple 4.4.4.

programme Maple 4.4.4: Etude de la stabilité du schèma ADI (4.44)

> restart:
# Equation de convection-diffusion
> diff(U(x,y,t),t)+V1*diff(U(x,y,t),x)+V2*diff(U(x,y,t),y)=
> kappa*(diff(U(x,y,t),x$2)+diff(U(x,y,t),y$2));eq:=%:
# Schema ADI
> (U[ns,i,j]-U[n,i,j])/(dt/2)+V1*(U[ns,i+1,j]-U[ns,i-1,j])/
 (2*dx)+V2*(U[n,i,j+1]-U[n,i,j-1])/(2*dy)=kappa*((U[ns,i+1,j]
 -2*U[ns,i,j]+U[ns,i-1,j])/dx^2+(U[n,i,j+1]-2*U[n,i,j]+
 U[n,i,j-1])/dy^2);   eqh1:=%:
> (U[n+1,i,j]-U[ns,i,j])/(dt/2)+V1*(U[ns,i+1,j]-U[ns,i-1,j])/
 (2*dx)+V2*(U[n+1,i,j+1]-U[n+1,i,j-1])/(2*dy)=kappa*((U[ns,i+1,j]
 -2*U[ns,i,j]+U[ns,i-1,j])/dx^2+(U[n+1,i,j+1]-2*U[n+1,i,j]+
 U[n+1,i,j-1])/dy^2);  eqh2:=%:
# Etude de la Stabilite
> Up:=(n,i,j)->Psi[n]*exp(I*omega[1]*i*dx)*exp(I*omega[2]*j*dy);
> subs(U[ns,i,j]=Up(ns,i,j),U[n,i,j]=Up(n,i,j),
      U[ns,i+1,j]=Up(ns,i+1,j),U[ns,i-1,j]=Up(ns,i-1,j),
      U[n,i,j+1]=Up(n,i,j+1),U[n,i,j-1]=Up(n,i,j-1),eqh1):
> rel1:=simplify(%*exp(-I*omega[1]*i*dx)*exp(-I*omega[2]*j*dy)):
> Psi[ns]/Psi[n]=solve(subs(Psi[ns]=G*Psi[n],rel1),G);rel11:=%:
#
> subs(U[ns,i,j]=Up(ns,i,j),U[n+1,i,j]=Up(n+1,i,j),
    U[ns,i+1,j]=Up(ns,i+1,j),U[ns,i-1,j]=Up(ns,i-1,j),
    U[n+1,i,j+1]=Up(n+1,i,j+1),U[n+1,i,j-1]=Up(n+1,i,j-1),eqh2):
> rel2:=simplify(%*exp(-I*omega[1]*i*dx)*exp(-I*omega[2]*j*dy)):
> Psi[n+1]/Psi[ns]=solve(subs(Psi[n+1]=G*Psi[ns],rel2),G);
  rel22:=%:
# Facteur d'amplification
> G:=rhs(rel11)*rhs(rel22);
# Etude de chaque terme
> A1=r1*(1-cos(omega[1]*dx));A2=r2*(1-cos(omega[2]*dy));
> rel1:={%,%%}: rel11:={r1=kappa*dt/dx^2,r2=kappa*dt/dy^2}:
> B1=CFL1*sin(omega[1]*dx)/2;B2=CFL2*sin(omega[2]*dy)/2;
> rel2:={%,%%}: rel22:={CFL1=V1*dt/dx,CFL2=V2*dt/dy}:
> 'G'=((A1-1)+I*B1)/((A1+1)+I*B1)*((A2-1)+I*B2)/((A2+1)+I*B2);
  GG:=rhs(%):
# Verification
> G=subs(rel1,rel2,rel11,rel22,GG):simplify(%):rhs(%)-lhs(%);
# Calcul du carre du module de G
> GG;
> G1:=((A1-1)+I*B1)/((A1+1)+I*B1);
  G1M=((A1-1)^2+B1^2)/((A1+1)^2+B1^2);
> G2:=((A2-1)+I*B2)/((A2+1)+I*B2);
  G2M=((A2-1)^2+B2^2)/((A2+1)^2+B2^2);
# G1M et G2M sont donc <1 ==> donc stabilite

On définit les deux étapes (4.44 et 4.45) du schéma ADI (lignes 6 à 13), puis on introduit une perturbation $U_{p}$ décomposée en mode de Fourier (ligne 15), que l'on introduit dans les deux équations (lignes 16 et 22). On calcul l'amplification du mode pour chacune de ces équations: i.e. $\frac{\psi^{*}}{\psi^{n}}$ pour la première (ligne 20) et $\frac{\psi^{n+1}}{\psi^{*}}$ pour la seconde (ligne 22). D'où l'on déduit le facteur d'amplification global du schéma $G=\frac{\psi^{n+1}}{\psi^{n}}$ (ligne 29):

$$G=\frac{a_{1}-1+Ib_{1}}{a_{1}+1+Ib_{1}} \frac{a_{2}-1+Ib_{2}}{a_{2}+1+Ib_{2}}$$

avec $a_{1}=\frac{\kappa dt}{dx^{2}}(1-\cos\omega_{1}dx)$ , $a_{2}=\frac{\kappa dt}{dy^{2}}(1-\cos\omega_{2}dy)$ , $b_{1}=\frac{V_{1}dt}{dx}$ , $b_{2}=\frac{V_{2}dt}{dx}$ .

Le carré du module de $G$ vaut:

$$G^{2}=\frac{(a_{1}-1)^{2}+b_{1}^{2}}{(a_{1}+1... ...{2}}{(a_{2}+1)^{2}+b_{2}^{2}}\leq1 \forall\omega_{1},\omega_{2}$$

Il est plus petit que 1 puisque $a_{1}$ et $a_{2}$ sont positifs.

Le schéma ADI (4.44 et 4.45) est donc inconditionnellement stable.

4.4.4.1.2 Étude de la consistance:

Pour étudier la consistance, on utilise le programme Maple 4.4.4, qui est la suite du programme précédent 4.4.4.

programme Maple 4.4.4: Etude de la consitance du schèma ADI (4.44)

# Erreur de troncature dvt de taylor autour de t=n+1/2
# Equation equivalente
> 1/2*eqh1+1/2*eqh2:eqh12:=lhs(%)-rhs(%):
# subsitution de la solution exacte dans cette equation
> Uex:=(r,p,q)->U(x+(p-i)*dx,y+(q-j)*dy,t+(r-n-1/2)*dt);
> subs(U[n,i,j]=Uex(n,i,j),
     U[n,i+1,j]=Uex(n,i+1,j),U[n,i-1,j]=Uex(n,i-1,j),
     U[n,i,j+1]=Uex(n,i,j+1),U[n,i,j-1]=Uex(n,i,j-1),
     U[n+1,i,j]=Uex(n+1,i,j),
     U[n+1,i+1,j]=Uex(n+1,i+1,j),U[n+1,i-1,j]=Uex(n+1,i-1,j),
     U[n+1,i,j+1]=Uex(n+1,i,j+1),U[n+1,i,j-1]=Uex(n+1,i,j-1),
     U[ns,i,j]=Uex(n+1/2,i,j),
     U[ns,i+1,j]=Uex(n+1/2,i+1,j),U[ns,i-1,j]=Uex(n+1/2,i-1,j),
     U[ns,i,j+1]=Uex(n+1/2,i,j+1),U[ns,i,j-1]=Uex(n+1/2,i,j-1),
     eqh12);   eqh3:=%:
# Developpement de Taylor autour de t=n+1/2
> k:=6:
> U(x,y,t-dt/2)=convert(mtaylor(U(x,y,t-dt/2),[dt],k),diff):
  S1:=%:
> U(x,y+dy,t-dt/2)=convert(mtaylor(U(x,y+dy,t-dt/2),[dy,dt],
  k),diff): S2:=%:
> U(x,y-dy,t-dt/2)=convert(mtaylor(U(x,y-dy,t-dt/2),[dy,dt],
  k),diff):S3:=%:
> U(x+dx,y,t-dt/2)=convert(mtaylor(U(x+dx,y,t-dt/2),[dx,dt],
  k),diff):S4:=%:
> U(x-dx,y,t-dt/2)=convert(mtaylor(U(x-dx,y,t-dt/2),[dx,dt],
  k),diff):S5:=%:
> U(x,y+dy,t)=convert(mtaylor(U(x,y+dy,t),[dy],k),diff):
  S6:=%:
> U(x,y-dy,t)=convert(mtaylor(U(x,y-dy,t),[dy],k),diff):
  S7:=%:
> U(x+dx,y,t)=convert(mtaylor(U(x+dx,y,t),[dx],k),diff):
  S8:=%:
> U(x-dx,y,t)=convert(mtaylor(U(x-dx,y,t),[dx],k),diff):
  S9:=%:
> U(x,y,t+dt/2)=convert(mtaylor(U(x,y,t+dt/2),[dt],k),diff):
  S10:=%:
> U(x,y+dy,t+dt/2)=convert(mtaylor(U(x,y+dy,t+dt/2),[dy,dt],
  k),diff):S11:=%:
> U(x,y-dy,t+dt/2)=convert(mtaylor(U(x,y-dy,t+dt/2),[dy,dt],
  k),diff):S12:=%:
> U(x+dx,y,t+dt/2)=convert(mtaylor(U(x+dx,y,t+dt/2),[dx,dt],
  k),diff):S13:=%:
> U(x-dx,y,t+dt/2)=convert(mtaylor(U(x-dx,y,t+dt/2),[dx,dt],
  k),diff): S14:=%:
# substitution dans l'équation discrete - equation exacte
> subs(S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,eqh3)-
  (lhs(eq)-rhs(eq)):
> simplify(%):collect(%,{dt,dx,dy});
# Schema d'ordre 2 en dt,dx et dy

Pour cela on fait la demi somme des deux équations (4.44 et 4.45) pour obtenir une équation discrète équivalente à l'équation exacte (4.26) (ligne 3):

$$\frac{u_{i,j}^{n+1}-u_{i,j}^{n}}{dt} = -V_{1}\left(\frac{u_{i+1,j}^{*}-u_{i-1,j}^{*}}{2dx}\right)$$ (4.52)

$$- \frac{V_{2}}{2}\left(\frac{u_{i,j+1}^{n}-u_{i,j-1}^{n}}{2dy}+\frac{u_{i,j+1}^{n+1}-u_{i,j-1}^{n+1}}{2dy}\right)$$

$$+ \kappa\left(\frac{u_{i+1,j}^{*}-2u_{i,j}^{*}+u_{i-1,j}^{*}}{dx^{2}}\right)$$

$$+ \frac{\kappa}{2}\left(\frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{dy^{2}}+\frac{u_{i,j+1}^{n+1}-2u_{i,j}^{n+1}+u_{i,j-1}^{n+1}}{dy^{2}}\right)$$

En comparant cette équation (4.52) au schéma de Cranck Nicholson (4.37), on peut retrouver ce dernier en remplaçant simplement dans (4.52) les valeurs de $u_{i,j}^{*}$ par la moyenne $\frac{1}{2}(u_{i„j}^{n}+u_{i,j}^{n+1})$ . Ce schéma ADI est donc bien équivalent au schèma de Cranck Nicolson. Les valeurs intermédiaires $u_{i,j}^{*}$ sont des approximations de la solution à $t=(n+\frac{1}{2})dt$ .

Dans l'équation équivalente (4.52), on substitue la solution approchée par la solution exacte (lignes 6 à 15), et on effectue des développements limités autour de $u_{i,j}^{*}=u(idx,jdy,(n+\frac{1}{2}))$ (lignes 18 à 45). Après soustraction de l'équation exacte, on obtient l'erreur de troncature $E_{t}$ , qui est en $O(dt^{2},dx^{2},dy^{2})$ .

Le schéma ADI (4.44 et 4.45) est donc consistant avec l'équation exacte (4.26), et est d'ordre 2 en temps et en espace, i.e. en $O(dt^{2},dx^{2},dy^{2})$ .

4.4.5 Expérimentation numérique avec Matlab

Le programme Matlab (4.4.5) implémente le schéma ADI en utilisant une programmation matricielle.

programme Matlab 4.4.5: Résolution numérique de l'équation (4.27)

% resolution ADI
clear
% maillage
L=1; H=1;
Nx=51; dx=L/(Nx-1); X=[0:dx:L];
Ny=51; dy=H/(Ny-1); Y=[0:dy:H];
% champ vitesse
v1=1; v2=1;
% parametres
kappa=0.01; dt=0.001; Tf=0.5; nit=round(Tf/dt)
% champ initial
delta=(0.1)^2; x0=L/4; y0=H/4;
Ui=exp(-(X-x0).^2/delta)'*exp(-(Y-y0).^2/delta);
Un=Ui; Us=Ui; Un1=Ui;
% matrices 3D
I1=ones(1,Nx); I2=ones(1,Ny);
A1 =[-kappa/dx^2-v1/(2*dx);  2/dt+2*kappa/dx^2; ...
     -kappa/dx^2+v1/(2*dx)]*I1;
C1=[  kappa/dy^2+v2/(2*dy);  2/dt-2*kappa/dy^2; ...
      kappa/dy^2-v2/(2*dy)]*I1;
A2 =[-kappa/dy^2-v2/(2*dy);  2/dt+2*kappa/dy^2; ...
     -kappa/dy^2+v2/(2*dy)]*I2;
C2=[  kappa/dx^2+v1/(2*dx);  2/dt-2*kappa/dx^2; ...
      kappa/dx^2-v1/(2*dx)]*I2;
% C.L de Dirichlet en i=1 j=1
A1(:,1)=0; A1(2,1)=1; C1(:,1)=0;
A2(:,1)=0; A2(2,1)=1; C2(:,1)=0;
% C.L de Neumann sur les frontieres i=Nx j=Ny
A1(1,Nx)=A1(1,Nx)+A1(3,Nx);  A1(3,Nx)=0;
A2(1,Ny)=A2(1,Ny)+A2(3,Ny);  A2(3,Ny)=0;
% iterations en temps
for it=1:nit
    % 1ere etape ADI
    Us(1:Nx,1)=0;
    for j=2:Ny-1
      B1=C1(1,1:Nx)'.*Un(1:Nx,j-1)+C1(2,1:Nx)'.*Un(1:Nx,j)+...
         C1(3,1:Nx)'.*Un(1:Nx,j+1);
      Us(1:Nx,j)=tridiag(A1,B1);
    end;
    B1=C1(2,1:Nx)'.*Un(1:Nx,Ny)+(C1(1,1:Nx)+...
       C1(3,1:Nx))'.*Un(1:Nx,Ny-1);
    Us(1:Nx,Ny)=tridiag(A1,B1);
    % 2nd etape ADI
    Un1(1,1:Ny)=0;
    for i=2:Nx-1
      B2=C2(1,1:Ny)'.*Us(i-1,1:Ny)'+C2(2,1:Ny)'.*Us(i,1:Ny)'+...
         C2(3,1:Ny)'.*Us(i+1,1:Ny)';
      Un1(i,1:Ny)=tridiag(A2,B2)';
    end;
    B2=C2(2,1:Ny)'.*Us(Nx,1:Ny)'+...
       (C2(1,1:Ny)+C2(3,1:Ny))'.*Us(Nx-1,1:Ny)';
    Un1(Nx,1:Ny)=tridiag(A2,B2)';
    % iteration suivante
    Un=Un1;
end;

Les paramètres du calcul sont définis aux lignes 4 à 10. Les matrices tridiagonales sont construites sur les lignes 17 à 24, puis on applique les conditions aux limites (lignes 26 à 30).

La boucle en temps (lignes 32 à 55) inclus les deux étapes ADI et utilise la fonction tridiag (lin3D) pour la résolution des systèmes linéaires tri-diagonaux.

4.4.5.0.1 Mode propre de diffusion:

Pour valider ce programme, nous avons tout d'abord simuler la diffusion du mode propre (4.36) $p=1$ et $q=1$ avec un maillage de $N_{x}=N_{y}=51$ points dans chaque direction et un paramètre $\kappa=0.01$ . La solution calculée au bout d'un temps $T_{f}=0.5$ avec $dt=0.01$ est tracée sur la figure (4.19). L'allure de la solution (figure 4.19a) coïncide bien avec la solution exacte (figure 4.17a) , ce que confirme le tracé de l'évolution temporelle de la solution au point $x=1$ , $y=1$ comparée à la solution exacte (figure 4.19b).

Figure 4.19: solution de diffusion avec le schéma ADI (mode propre $p=1$ , $q=1$ )

[$t=0.5$ ]

[Amplitude en x=1 et y=1]

Pour tester la précision d'intégration en temps du schéma, nous avons calculer l'erreur au point ( $x=1$ , $y=1$ ) au bout du temps $T_{f}=4$ (de l'ordre du temps caractéristique de diffusion), en fonction du pas d'intégration en temps $dt$ . Le résultat de la figure (4.20a) montre que pour les pas en temps choisis l'erreur est quasiment indépendante du pas en temps $dt$ , et est donc essentiellement une erreur de discrétisation spatiale. On note que les pas en temps choisis sont tels que le pas en temps est beaucoup plus faible que le temps caractéristique de diffusion: $dt\ll\tau_{d}\approx4.5$

Figure 4.20: Erreur numérique du schéma ADI (cas de diffusion)

[Erreur temporelle $N=51$ ]

[erreur spatiale $dt=0.1$ ]

Nous l'avons vérifié en faisant varier le nombre de points du maillage de $N_{x}=N_{y}=11$ à $N_{x}=N_{y}=101$ avec un pas en temps $dt=0.1$ fixé. L'évolution de l'erreur en fonction du pas de discrétisation spatiale $h=dx=dy$ est tracée sur la figure (4.20b), et on constate que l'erreur décroît en $O(h^{2})$ .

Pour cette condition initiale, nous avons aussi fait une simulation avec une vitesse de convection non nulle $(V_{1}=V_{2}=1)$ . Ce cas correspond à un nombre de Péclet $P_{e}=200$ . Les iso-valeurs de la solution sont tracées sur la figure (4.21).

Figure 4.21: Iso-valeurs de la solution (cas de convection-diffusion)

[$t=0$ ]

[$t=0.25$ ]

[$t=0.5$ ]

On note la convection sans déformation de la solution initiale, ce qui confirme que la condition aux limites sur $\Gamma_{1}$ autorise la sortie des structures hors du domaine. On a aussi comparé l'évolution temporelle du maximum de la solution exacte de diffusion et du maximum de la solution calculée, que l'on a tracé sur la figure (4.22a). On vérifie ainsi que la décroissance de la solution est une décroissance visqueuse.

Figure 4.22: Erreur numérique du schéma ADI

[Evolution temporelle]

[Erreur fonction de $dt$ ]

Nous avons ensuite étudié l'influence du pas d'intégration en temps $dt$ , en traçant sur la figure (4.22b) l'écart en fonction de $dt$ entre le maximum de la solution exacte de diffusion et le maximum de la solution calculée. On constate que cet écart croît rapidement en dessus d'une valeur $dt=0.1$ . Cette valeur est justement de l'ordre de grandeur du temps caractéristique de convection $\tau_{c}\approx0.2$ .

En conclusion sur cette simulation, on note que le choix des paramètres numériques a été fixé par la physique du problème, et non par des conditions numériques de stabilité:

1. la maillage $N_{x}=N_{y}=51$ permet de décrire finement la condition initiale
2. le pas en temps $dt=0.01$ vérifie $dt\ll\tau_{d}$ et $dt\ll\tau_{c}$ .

4.4.5.0.2 Convection d'une gaussienne:

Le second cas de calcul correspond à la condition initiale gaussienne (4.31), avec les mêmes paramètres qu'au paragraphe c4gpar.

Figure 4.23: Convection d'une gaussienne

[solution numérique à $t=0$ et $t=0.5$ ]

[Amplitude]

Pour un maillage $N_{x}=N_{y}=51$ et un pas en temps $dt=0.01$ , on a tracé sur la figure (4.23a) la solution à $t=0$ et à $t=0.5$ . Cette solution se compare très bien avec la solution exacte tracée sur la figure (4.16). On a aussi comparé l'évolution temporelle de l'amplitude de la tache gaussienne calculée avec le schéma ADI avec l'expression analytique (4.34). Ces deux courbes coïncident, ainsi que le montre la figure (4.23b).

Figure 4.24: Convection d'une gaussienne $(P_{e}^{h}=20)$

[isovaleurs à $P_e^h=20$ ]

[isovaleurs à $P_e^h=2$ ]

[profil à $x=0.75$ ]

Nous avons effectué une seconde simulation avec un coefficient de diffusion plus petit $\kappa=0.001$ . Dans ce cas, la solution numérique présente des oscillations (figure 4.24). Pour cette valeur de $\kappa$ , le nombre de Péclet de maille:

$$P_{e}^{h}=\frac{V h}{\kappa}$$

vaut $P_{e}^{h}=20$ au lieu de $P_{e}^{h}=2$ avec la valeur de $\kappa$ précédente. Sur le tracé des iso-valeurs (4.24a), on constate l'apparition de légères oscillations, caractérisées par la présence de nombreuses lignes iso-valeurs $u=0$ , qui n'existent pas à $P_{e}^{h}=2$ (figure 4.24b). Le tracé d'un profil à $x=0.75$ (figure 4.24c) montre bien l'apparition d'une oscillation numérique au pied de la tache gaussienne. Ces oscillations numériques sont de même nature que celles étudiées au paragraphe c3centre du chapitre précédent. Elles apparaissent dès que le Péclet de maille $P_{e}^{h}$ devient plus grand que 2 et indiquent que le maillage n'est plus suffisamment fin pour capter la solution de convection avec ce schéma centré.