2.7 Analyse de Fourier

L'analyse de Fourier permet une étude des propriétés d'un schéma aux différences finies, en particulier dans le cadre une équation parabolique, comme l'équation de convection diffusion:

$$\frac{\partial u}{\partial t}+V\frac{\partial u}{\partial x}-\nu\frac{\partial^{2}u}{\partial x^{2}}=0$$ (2.55)

2.7.1 analyse de Fourier de la solution exacte

On cherche une solution sous la forme d'une serie de Fourier en espace

$$u(x,t)=\sum_{k}U_{k}(t)e^{I k x}$$

$k$ est le nombre d'onde: $k=\frac{2\pi}{\lambda}$ (où $\lambda$ est la longueur d'onde)

L'équation d'évolution du mode de Fourier $U_{k}$ s'écrit:

$$\frac{dU_{k}}{dt}+I V k U_{k}+\nu k^{2} U_{k}=0$$ (2.56)

2.7.2 Symbol G de l'équation

Cette équation s'intégre entre $t$ et $t+\Delta t$ , en supposant $V$ et $\nu$ constant sur l'intervalle de temps:

$$U_{k}(t+\Delta t)=U_{k}(t) e^{-(I V k+\nu k^{2})\Delta t}$$

On peut donc écrire la relation entre le mode de Fourier $U_{k}$ à l'instant $t$ et $t+\Delta t$ :

$$U_{k}(t+\Delta t)=\mathcal{G}_{k}(\Delta t) U_{k}(t)$$  $${avec } G_{k}(\Delta t)=e^{-(I V k+\nu k^{2})\Delta t}$$

$\mathcal{G}_{k}$ est appelé le symbol de l'opérateur (ou équation) (2.55).

Ce symbol s'écrit sous la forme générale

$$\mathcal{G}_{k}(\Delta t)=e^{-I\Omega_{k}\Delta t-\Phi_{k}\Delta t}$$ (2.57)

et permet de calculer la solution à l'étape $t+\Delta t$ en fonction de la solution à l'étape $t$

$$U_{k}(t+\Delta t)=\mathcal{G}_{k}(\Delta t) U_{k}(t)$$

• $\Phi_{k}$ est le coefficient de dissipation du schéma et est donc strictement positif si l'équation est dissipative. Pour l'équation (2.55) $\Phi_{k}=\nu k^{2}$ , et l'équation est dissipative si $\nu>0$ .
• $\Omega_{k}$ définie la vitesse de propagation $C$ de l'onde de nombre d'onde $k$ : $\Omega_{k}=C k$ . L'équation de transport est dispersive si la vitesse de propagation $C$ dépend de $k$ . Pour l'équation (2.55) $\Omega_{k}=k V$ , donc $C=V$ et cette équation est donc non dispersive.

2.7.3 Interprétation

Si on cherche la solution de (2.55) sous la forme d'une onde progressive:

$$u(x,t)=A(t) e^{I(kx-\omega t)}$$soit  $$U_{k}=A(t)e^{-I\omega t}$$

pour un symbol $G_{k}$ donnée par la relation (2.57), on en déduit

$$\omega=\Omega_{k}, A(t)=A(0)e^{-\Phi_{k}t}$$

L'onde est donc une onde progressive qui se propage avec une vitesse $\Omega_{k}/k$ et qui se trouve amortie exponentiellement:

$$u(x,t)=A_{0}e^{-\Phi_{k}t}e^{I(kx-\Omega_{k}t)}$$

englishANIMATION onde progressive

2.7.4 Analyse de Fourier du schéma D.F.

Pour analyser le schéma D.F., on éffectue la transformée de Fourier des valeurs $u_{i}^{n}$ dans l'équation aux différences:

$$u_{i}^{n}=\sum_{k}U_{k}^{h}(n) e^{I k i\Delta x}$$

d'où l'on déduit une relation pour le mode $k$

$$U_{k}^{h}(n+1)=\mathcal{G}_{k}^{h}(\Delta t,\Delta x) U_{k}^{h}(n)$$ (2.58)

$\mathcal{G}_{k}^{h}$ est appelée le symbol du schéma D.F. (équivalent discret de $\mathcal{G}_{k}$ ).

Ce symbol $\mathcal{G}_{k}^{h}$ peut sécrire sous la forme

$$\mathcal{G}_{k}^{h}=e^{-I\Omega_{k}^{h}\Delta t-\Phi_{k}^{h}\Delta t}$$

avec

$$\Omega_{k}^{h}=-\frac{\arg(\mathcal{G}_{k}^{h... ...hi_{k}^{h}=-\frac{\ln\left\vert\mathcal{G}_{k}^{h}\right\vert}{\Delta t}$$

Si le schéma D.F. approchait exactement l'équation, on aurait $\mathcal{G}_{k}^{h}=\mathcal{G}_{k}$ , ce qui n'est évidemment pas le cas. On définit donc

• l'erreur d'amplitude $E_{\Phi}$ du schéma D.F.:
  
  $$E_{\Phi}=\left\Vert \mathcal{G}_{k}^{h}\right\Vert -\left\Vert \mathcal{G}_{k}\right\Vert =e^{-\Delta t \Phi_{k}^{h}}-e^{-\Delta t \Phi_{k}}$$

• l'erreur de phase $E_{\Omega}$ du schéma D.F.:
  
  $$E_{\Omega}=\arg(\mathcal{G}_{k}^{h})-\arg(\mathcal{G}_{k})=-\Delta t (\Omega_{k}^{h}-\Omega_{k})$$

La solution d'onde progressive du schéma D.F. est donc une onde qui se propage avec une vitesse $C=\frac{\Omega_{k}^{h}}{k}$ et avec un amortissement exponentielle de coefficient $\Phi_{k}^{h}$ .

Si à l'instant $t^{n}$ , la solution approchée coincide avec la solution exacte. A l'instant $t^{n+1}$ , l'onde calculée aura subi un déphase $E_{\Omega}$ par rapport à la solution exacte, ainsi qu'un amortissement de $E_{\Phi}$ par rapport à cette même onde. La solution calculée sera donc en retard par rapport à la solution exacte si $E_{\Omega}$ est positif, et en avance sinon.

On note aussi, que la relation d'ondes (2.58) est aussi valable pour une perturbation $\epsilon_{i}^{n}=\Psi_{k}^{n}e^{I k i\Delta x}$ , et fournit donc l'évolution de l'amplitude et la phase des perturbations:

$$\Psi_{k}^{n+1}=\mathcal{G}_{k}^{h}(\Delta t,\Delta x) \Psi_{k}^{n}$$

Le module du symbol $\mathcal{G}_{k}^{h}$ de l'équation D.F. est donc le facteur d'amplification du schéma

$$G=\left\Vert \mathcal{G}_{k}^{h}\right\Vert$$

stabilité:

le schéma au D.F. de symbol $\mathcal{G}_{k}^{h}$ est stable (au sens de Von Neumann), si

$$\left\Vert \mathcal{G}_{k}^{h}\right\Vert <1 \forall k$$

consistance:

le schéma au D.F. de symbol $\mathcal{G}_{k}^{h}$ est consistant à l'ordre (n,m) avec l'équation exacte de symbol $\mathcal{G}_{k}$ si:

$$\left\Vert \frac{\mathcal{G}_{k}^{h}-\mathcal{G}_{k}}{\Delta t}\right\Vert =\theta(\Delta x^{n},\Delta t^{m}) \forall k$$

l'ordre du schéma est alors d'ordre $n$ en espace et $m$ en temps.

2.7.5 Analyse d'un schéma explicite centré

Une approximation D.F. de l'équation de convection-diffusion (2.55) est le schéma D.F. centré suivant:

$$\frac{\Delta_{t}u_{i}^{n}}{dt}+V\frac{\overline{\delta}_{x}u_{i}^{n}}{h}-\nu\frac{\delta_{x}^{2}u_{i}^{n}}{h^{2}}=0$$

soit

$$\frac{u_{i}^{n+1}-u_{i}^{n}}{dt}+V\frac{u_{i+1}^{n}-u_{i-1}^{n}}{2h}-\nu\frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{h^{2}}=0$$ (2.59)

On introduit les 2 nombres sans dimensions suivants:

$$\alpha=\frac{V dt}{h}, \beta=\nu\frac{dt}{h^{2}}$$

2.7.5.1 stabilité

Le symbol du schéma s'écrit:

$$\mathcal{G}_{k}^{h}=1+2\beta(\cos kh-1)-I\alpha\sin kh$$

Son module s'écrit en fonction de $y=\cos kh$ :

$$G^{2}=\left\Vert \mathcal{G}_{k}^{h}\right\Vert ^{2}=(1+2\beta(y-1))^{2}+\alpha^{2}(1-y^{2})$$

C'est un polynôme de degré 2 en $y$ que l'on étudie entre $[-1,1]$ où il est positif. Il vérifie:

$$G^{2}(-1)=(4\beta-1)^{2}, G(1)=1$$

Sa dérivée s'écrit:

$$\frac{dG^{2}}{dy}=(8\beta^{2}-2\alpha^{2})y-4\beta(2\beta-1)$$

et vaut en $y=1$

$$\frac{dG^{2}}{dy}(1)=4\beta-2\alpha^{2}$$

Pour que $G^{2}(y)<1$ pour $y\in[-1,1]$ , il est donc suffisant que:

$$G^{2}(-1)\leq1$$et $$\frac{dG^{2}}{dy}(1)\ge0$$

soit

$$2\beta\ge\alpha^{2}$$  et $$2\beta\le1$$

La condition de stabilité conduit à

$$\alpha^{2}\leq2\beta\leq1$$ (2.60)

On retrouve la condition de stabilité pour l'équation de la chaleur:

$$\beta\leq\frac{1}{2} \Rightarrow \frac{\nu dt}{h^{2}}\leq\frac{1}{2}$$ (2.61)

avec une condition supplémentaire sur le pas en temps $dt$ :

$$\alpha^{2}\leq2\beta \Rightarrow \frac{V^{2} dt}{\nu}\leq2$$

que l'on peut écrire en fonction du nombre de Courant $\alpha$ et d'un nombre de Pechlet $P_{e}^{h}=\frac{V h}{\nu}=\frac{\alpha}{\beta}$ :

$$\frac{V^{2} dt}{\nu}=\frac{V dt}{h} \frac{V h}{\nu}=\alpha P_{e}^{h}\leq2$$ (2.62)

2.7.5.2 consistance

En posant $\theta=kh,$ on peut calculer un développement en série de Taylor de $\mathcal{G}_{k}^{h}(\theta)$ en $\theta$ :

$$\Omega_{k}^{h} dt=-\arg(\mathcal{G}_{h}^{k})=\alpha\theta+O(\alpha\theta^{3},\alpha\beta\theta^{3})$$

$$\Phi_{k}^{h} dt=-\ln\left\Vert \mathcal{G}_... ...\frac{\alpha^{2}}{2})\theta^{2}+O(\beta\theta^{4},\beta^{2}\theta^{4})$$

soit en fonction de $h$ et $dt$ :

$$\Omega_{k}^{h}=kV+O(h^{2},dt)$$

$$\Phi_{k}^{h}=(\nu-\frac{V^{2}dt}{2})k^{2}+O(h^{2},dt)$$

d'où l'expression de $\mathcal{G}_{k}^{h}$

$$\mathcal{G}_{k}^{h}=e^{-IkVdt}e^{-(\nu-\frac{V^{2}dt}{2})k^{2}dt}+dt O(h^{2},dt)$$

Le schéma est donc consistant et d'ordre $O(h^{2},dt)$ puisque:

$$\frac{\mathcal{G}_{k}^{h}-\mathcal{G}_{k}}{dt}=O(h^{2},dt)$$

L solution approchée d'onde progressive se déplace (pour les petits nombres d'ondes $k$ , i.e. pour $\theta$ petit) avec une célérité $V$ et est amortie avec une coefficient $\nu-\frac{V^{2}dt}{2}$ , qui doit donc être positif. On retrouve ainsi une des conditions de stabilité (2.62):

$$\beta-\frac{\alpha^{2}}{2}\ge0$$

2.7.5.3 condition de positivité

Si la solution exacte vérifie un principe du maximum:

$$\min_{x\in[a,b]}u(x,t^{n})\leq u(x,t^{n+1})\leq\max_{x\in[a,b]}u(x,t^{n})$$

on souhaite que la solution approchée $u_{i}^{n}$ vérifie un principe identique. Pour cela on peut imposer au schéma aux D.F. d'être positif.

schéma positif:

un schéma au D.F. est dit positif s'il existe des réels positifs $a_{l}$ tels que:

$$u_{i}^{n+1}=\sum_{l}a_{l}u_{i}^{n}$$

Dans ce cas la solution discréte vérifie:

$$\min_{j}u_{j}^{n}\leq u_{i}^{n+1}\leq\max_{j}u_{j}^{n}$$

et est donc $L^{\infty}$ décroissante ( $\left\Vert u_{j}^{n}\right\Vert _{\infty}=\max_{j}(\left\vert u_{j}^{n}\right\vert)$ :

$$\left\Vert u_{i}^{n+1}\right\Vert _{\infty}\l... ...i}^{n}\right\Vert _{\infty}\leq\left\Vert u_{i}^{0}\right\Vert _{\infty}$$

La condition de positivité pour le schéma (2.59) s'écrit:

$$u_{i}^{n+1}=u_{i-1}^{n}(\frac{\alpha}{2}+\beta)+u_{i}^{n}(1-2\beta)+u_{i+1}^{n}(\beta-\frac{\alpha}{2})$$

ce qui implique:

$$\alpha\leq2\beta\leq1$$ (2.63)

Cette condition impose donc les 2 relations suivantes:

$$\beta\leq\frac{1}{2} \Rightarrow \frac{\nu dt}{h^{2}}\leq\frac{1}{2}$$

$$\alpha\leq2\beta \Rightarrow P_{e}^{h}=\frac{V h}{\nu}\leq2$$

La condition

est plus restrictive que la condition de stabilité (2.60), et c'est donc celle que l'on retiendra pour éviter les oscillations dans la solution.

Cette relation impose:

$$\beta=\frac{\nu dt}{h^{2}}\leq\frac{1}{2}$$

et une condition sur le nombre de Pechlet

$$P_{e}^{h}=\frac{V h}{\nu}\leq2$$ (2.64)

Avec la condition (2.61), cette dernière condition vérifie la condition de stabilite (2.61)+(2.62), puisqu'elle impose:

$$\alpha\leq2\beta et 2\beta\leq1 \Rightarrow \alpha\leq1$$

et donc la condition (2.61)+(2.62):

$$\alpha^{2}\leq\alpha\leq2\beta et 2\beta\leq1$$

C'est donc les conditions (2.61)+(2.64) que l'on retiendra pour éviter les oscillations dans la solution:

$$\beta\leq\frac{1}{2}, P_{e}^{h}\leq2$$

On note que si $\nu=0$ , on a $\beta=0$ et le schéma D.F. explicite centré est inconditionnellement instable.

2.7.5.4 Expérimentation Matlab

On a calculé avec Matlab la solution approchée pour une condition initiale

$$u(x,0)=\sin k\pi x$$

et des conditions aux limites périodiques.

La solution exacte est alors une onde progressive amortie:

$$u_{e}(x,t)=e^{-\nu k^{2}\pi^{2}t}\sin(k\pi(x-Vt))$$

Sur un maillage de $N=200$ avec $\nu=1$ , pour le mode $k=6$ et $\beta=1/4$ , on a réalisé différentes simulations pour différentes valeurs de $V$ , i.e. du nombre de Pechlet $P_{e}^{h}$ .

Figure 2.2: solutions numériques pour $P_{e}=0.1$ et $1.0$

Figure 2.3: solutions numériques pour $P_{e}=2$ et $4$

On constate sur la figure (2.2), que pour des faibles valeurs de $P_{e}^{h}$ , la solution numérique approche très bien la solution exacte. L'erreur augmente quand le nombre de Pechlet augmente, la solution numérique étant beaucoup moins diffusive que la solution exacte. On remarque même que pour $P_{e}^{h}=4$ , qui est la limite de stabilité (i.e. $\alpha^{2}=2\beta$ ), la solution numérique ne diffuse plus ( $\Phi_{k}^{h}\approx0$ ).

Ce schéma est donc utilisable uniqument pour des problèmes peu convectifs !

englishANIMATION schéma explicite centré