4.3 Équation des ondes

4.3.1 Problème physique: oscillations d'une surface libre

On s'intéresse aux oscillations de la surface d'un liquide contenu dans un réservoir. Au repos le réservoir contient un liquide sur une hauteur $h_{0}$ . La surface du liquide est plane, horizontale et notée $\Omega$ . La répartition de pression dans le liquide est hydrostatique:

$$p_{s}=p_{0}+\rho_{0}g (h_{0}-z)$$

On perturbe la surface $h(x,y,t)$ à l'instant initiale. Celle ci se met alors à osciller de part et d'autre de sa position initiale $h_{0}$ . On néglige les effets de viscosité, et on applique les équations de conservation de la masse et de la quantité de mouvement à un cylindre élémentaire de base $dxdy$ et de hauteur $h(x,y,t)$ (figure 4.7).

Figure 4.7: oscillation de la surface libre

En notant $\overrightarrow{U}=\left[u(x,y,t),v(x,y,t)\right]$ la vitesse moyenne (suivant z) et $\rho$ la densité, le bilan de masse pour le cylindre de volume $hdxdy$ s'écrit:

$$\frac{\partial}{\partial t}(\rho h)+\frac{\partial\rho hu}{\partial x}+\frac{\partial\rho hv}{\partial y}=0$$

De même le bilan de quantité de mouvement s'écrit:

$$\frac{\partial\rho hu}{\partial t}+\frac{\partial\rho hu^{2}}{\partial x}+\frac{\partial\rho huv}{\partial y}+\frac{\partial p}{\partial x} =$$ 0

$$\frac{\partial\rho hv}{\partial t}+\frac{\partial\rho huv}{\partial x}+\frac{\partial\rho hv^{2}}{\partial y}+\frac{\partial p}{\partial y} =$$ 0

On considère que le fluide est un liquide incompressible $(\rho=cste)$ , et que la perturbation de la surface libre $\xi=h-h_{0}$ est faible. On peut alors considérer que la répartition de pression reste hydrostatique:

$$p(x,y,z,t)=p_{0}+\rho g (h-z)$$

En notant que $h=h_{0}+\xi$ , $p=p_{s}+\rho g\xi$ et en linéarisant les équations précédentes ( $\xi\ll h_{0}$ , $\rho g\xi\ll p_{s}$ , $U\ll1$ ), il vient:

$$\frac{\partial\xi}{\partial t}+h_{0}\frac{\partial u}{\partial x}+h_{0}\frac{\partial v}{\partial y} =$$ 0 (4.7)

$$\frac{\partial u}{\partial t}+g\frac{\partial\xi}{\partial x} =$$ 0

$$\frac{\partial v}{\partial t}+g\frac{\partial\xi}{\partial y} =$$ 0

En dérivant la première équation par rapport à $t$ , la seconde par rapport à $x$ et la troisième par rapport à $y$ , on obtiens l'équation de propagation de la perturbation $\xi$ de la surface libre en éliminant $u$ et $v$ :

$$\frac{\partial^{2}\xi}{\partial t^{2}}-gh_{0}\left(\frac{\partial^{2}\xi}{\partial x^{2}}+\frac{\partial^{2}\xi}{\partial y^{2}}\right)=0$$ (4.8)

C'est une équation des ondes qui traduit la propagation d'ondes de surface avec une célérité $c_{0}=\sqrt{gh_{0}}$ .

A cette équation, il faut ajouter une condition aux limites sur la frontière $\Gamma$ du réservoir. La condition physique est la condition de vitesse normale nulle $\overrightarrow{U}.\overrightarrow{n}=0$ sur $\Gamma$ . Pour obtenir une condition sur $\xi$ , on utilise la combinaison suivant $\overrightarrow{n}$ des 2 dernières équations (4.7):

$$\frac{\partial\overrightarrow{U}.\overrightarrow{n}}{\partial t}+g\frac{\partial\xi}{\partial n}=0$$

d'où l'on déduit la condition aux limites sur $\xi$ :

$$\left(\frac{\partial\xi}{\partial n}\right)_{\Gamma}=0$$ (4.9)

Les conditions initiales sont données par la déformation initiale $w(x,y)$ et la vitesse de cette déformation. En supposant la vitesse de déformation initiale nulle, on a donc:

$$\xi(x,y,0)=w(x,y) \„ \frac{\partial\xi}{\partial t}(x,y,0)=0$$ (4.10)

Le problème modèle associé s'écrit:
Trouvez $u(x,y,t)$ tel que

$$\frac{\partial^{2}u}{\partial t^{2}}-c_{0}^{2}(\frac{\partial^{2}u}{\partial^{2}x}+\frac{\partial^{2}u}{\partial^{2}y}) =0 \textrm{\mbox{ dans }}\Omega$$ (4.11)

$$\frac{\partial u}{\partial n}=0 \textrm{sur $\Gamma$}$$

$$u(t=0)=w \„ \frac{\partial u}{\partial t}(t=0)=0$$

Dans le cas d'un domaine $\Omega$ circulaire de rayon $r=1$ , on peut avantageusement passer en coordonnées polaires $(r,\theta)$ . Le problème modèle s'écrit:
Trouver $u(r,\theta,t)$ tel que:

$$\frac{\partial^{2}u}{\partial t^{2}}-c_{0}^{2}\left(\frac{1}{r}\f... ...}\right)+\frac{1}{r^{2}}\frac{\partial^{2}u}{\partial\theta^{2}}\right)=0 [0,1]*[0,2\pi]$$ (4.12)

$$\frac{\partial u}{\partial r}=0\textrm{ \mbox{ en } r=1}$$

$$u(t=0)=w \„ \frac{\partial u}{\partial t}(t=0)=0$$

4.3.2 Étude de solutions analytiques

Dans la cas cylindrique (4.12), on a cherché la solution analytique à l'aide du programme Maple (4.3.2).

programme Maple 4.3.2: Solution analytique de l'équation des ondes (4.12)

> restart;with(plots):
# Equation des ondes en polaire
> diff(u(r,theta,t),t$2)=c0^2/r*diff(r*diff(u(r,theta,t),r),r)
      +c0^2/r^2*diff(u(r,theta,t),theta$2);
> eq:=%:
# Solutions en variables séparées de type onde periodique en theta
> u(r,theta,t)=A(r)*exp(I*k*theta)*exp(I*lambda*c0*t);
> subs(%,eq):simplify(%/exp(I*k*theta)/exp(I*lambda*c0*t)/c0^2);
  eq1:=%:
> rhs(eq1)-lhs(eq1);eq2:=%:
> dsolve(eq2,A(r));
# Modes propres
> Ue:=(r,theta,t)->BesselJ(k,lambda*r)*exp(I*k*theta)*exp(I*lambda*c0*t);
> subs(r=1,diff(Ue(r,theta,t),r));
# Pour chaque valeur de k, on calcule les racines lambda
> -BesselJ(k+1,lambda)+k*BesselJ(k,lambda)/lambda=0; eq3:=lhs(%):
> plot({subs(k=0,eq3),subs(k=1,eq3),subs(k=2,eq3)},lambda=0..20);
> Um:=(k,p)->BesselJ(k,lambda[k,p]*r)*exp(I*k*theta)*exp(I*lambda[k,p]*c0*t);
> Uex=sum(sum(C[k,p]*Um(k,p),p=0..M),k=0..N);
# Conditions initiales particulieres (avec c0=1)
# 1/ mode radial
> lambda1:=fsolve(subs(k=0,eq3),lambda,10..12);
> Uex1:=BesselJ(0,lambda1*r)*cos(lambda1*t);
> animate3d([r*cos(theta),r*sin(theta),Uex1],r=0..1,theta=0..2*Pi,
     t=0..2*(2*Pi/lambda1),frames=50);
# 2/ mode angulaire
> lambda2:=fsolve(subs(k=1,eq3),lambda,6..10);
> Uex2:=BesselJ(1,lambda2*r)*cos(lambda2*t)*cos(theta);;
> animate3d([r*cos(theta),r*sin(theta),Uex2],r=0..1,theta=0..2*Pi,
      t=0..2*(2*Pi/lambda2),frames=50);

On a cherché une solution élémentaire à variables séparées sous la forme d'une onde se propageant à une célérité $\lambda c_{0}$ et périodique en $\theta$ (ligne 7):

$$u(r,\theta,t)=A(r) e^{Ik\theta}e^{I\lambda c_{0}t} k\in\mathbb{Z}$$

En reportant dans l'équation, on obtiens une équation de Bessel (ligne 9) pour $A(r)$ :

$$r^{2}\frac{d^{2}A}{dr^{2}}+r\frac{dA}{dr}+(r^{2}\lambda^{2}-k^{2})A=0$$

dont Maple nous fournit la solution générale avec la fonction dsolve (ligne 11). Cette solution générale est une combinaison linéaire de fonctions de Bessel de première et deuxième espèce: $BesselJ$ et $BesselY$ . De ces deux familles de fonctions de Bessel, on ne retient que la famille $BesselJ$ , qui est la seule a avoir une valeur finie (égale à 1) en $r=0$ . La solution élémentaire s'écrit (ligne 13)

$$U_{e}=BesselJ(k,\lambda r)e^{Ik\theta}e^{I\lambda c_{0}t}$$

Cette solution élémentaire doit vérifier la condition à la limite $\frac{\partial U_{e}}{\partial r}=0$ en $r=1$ (ligne 14), ce qui impose pour chaque valeur de $k$ des valeurs de $\lambda$ particulières. Les valeurs possibles de $\lambda$ sont les racines de la fonction $F_{k}$ (ligne 16):

$$F_{k}(\lambda)=-BesselJ(k+1,\lambda)+\frac{k BesselJ(k+1,\lambda)}{\lambda}$$

On a tracé cette fonction pour différentes valeur de k sur la figure (4.8). Pour une valeur de $k$ fixé, on a une infinité de racines $\lambda_{k,p}$ (c'est l'équivalent des racines $\frac{k\pi}{2}$ des fonctions $\cos\omega x$ en coordonnées cartésiennes). La solution élémentaire dépend donc de deux paramètres entiers $k$ et $p$ :

$$U_{e}(k,p)=BesselJ(k,\lambda_{k,p})e^{Ik\theta}e^{I\lambda_{k,p}c_{0}t}$$

Figure 4.8: fonction $F_{k}(\lambda )$ pour différentes valeurs de k

La solution générale de l'équation des ondes (4.12) est donc une combinaison linéaire de ces solutions élémentaires:

$$U_{ex}(r,\theta,t)=\sum_{k=0}^{\infty}\sum_{p=1}^{\infty}C_{k,p}BesselJ(k,\lambda_{k,p})e^{Ik\theta}e^{I\lambda_{k,p}c_{0}t}$$ (4.13)

où $\lambda_{k,p}$ est la $p^{\mbox{i\\lq {e}me}}$ racine de $F_{k}(\lambda)=0$ . Les coefficients $C_{k,p}$ permettent à $U_{ex}$ de vérifier les conditions initiales. A titre d'exemple, on a déterminé et tracé la solution modale:

$$BesselJ(k,\lambda_{k,p})\cos(k\theta)\cos(\lambda_{k,p}c_{0}t)$$ (4.14)

pour les deux cas particuliers:

1. $3^{\mbox{i\\lq {e}me}}$ mode radiale ($k=0,p=3$ et $\lambda_{0,3}\approx10,173468$ ), qui ne dépend donc pas de $\theta$ (lignes 22 à 24)
2. $3^{\mbox{i\\lq {e}me}}$ mode angulaire ($k=1,p=3$ et $\lambda_{1,3}\approx8,536316$ ) (lignes 27 à 29)

On a tracé ces deux modes propres sur la figure (4.9), et dans le programme Maple on a l'animation de ces modes avec la commande $animate3d$ .

Figure 4.9: modes propres de l'équation des ondes

[k=0]

[k=1]

4.3.3 Schéma explicite

La discrétisation de l'équation des ondes (4.11) avec un schéma explicite s'écrit en coordonnées cartésiennes sur un maillage régulier de pas $dx$ et $dy$ :

$$\frac{u_{i,j}^{n+2}-2u_{i,j}^{n}+u_{i,j}^{n-1}}{dt^{2}}=c_{0}^{2}... ...,j}^{n}}{dx^{2}}+\frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{dy^{2}}\right)$$ (4.15)

C'est l'extension naturelle du schéma explicite 1D (c3eq15) du chapitre précédent.

Pour discrétiser l'équation en coordonnées polaires (4.15), on discrétise le domaine polaire $\Omega=[0,1]*[0,2\pi]$ en $(r,\theta)$ avec un pas $dr$ et $d\theta$ , ce qui correspond à des points sur des rayons et des cercles dans le domaine physique $(x,y)$ (figure 4.10). On note $N_{r}$ et $N_{\theta}$ le nombre de noeuds suivant $r$ et $\theta$

Figure 4.10: maillage D.F. en coordonnées polaires

En notant $u_{i,j}^{n}=u(idr,jd\theta,ndt)$ les valeurs aux noeuds, la discrétisation par différences finies explicites de l'équation (4.12) s'écrit:

$$\frac{u_{i,j}^{n+2}-2u_{i,j}^{n}+u_{i,j}^{n-1}}{dt^{2}} = c_{0}^{2}\left(\frac{u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}}{dr^{2}}+\frac{1}{r_{i}}\frac{u_{i+1,j}^{n}-u_{i-1,j}^{n}}{2dr}\right.$$ (4.16)

$$+\left.\frac{1}{r_{i}^{2}}\frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{d\theta^{2}}\right)$$

4.3.3.0.1 conditions aux limites:

A cette équation, on ajoute la condition aux limites du problème (4.12) en $r=1$ : $\frac{\partial u}{\partial r}=0$ . Pour imposer cette condition, on utilise une condition miroir qui permet de calculer la valeur inconnue $u_{N_{r}+1,j}$ dans l'équation discrète sur $u_{N_{r},j}$ :

$$u_{N_{r}+1,j}=u_{N_{r}-1,j}$$ (4.17)

A cette condition physique, il faut ajouter des conditions numériques liées à la transformation en coordonnées polaires:

1. conditions de périodicité en $\theta$
   
   $$u(r,0)=u(r,2\pi)$$   et  $$\frac{\partial u}{\partial\theta}(r,0)=\frac{\partial u}{\partial\theta}(r,2\pi)$$

2. valeur finie en $r=0$ indépendante de $\theta$
   
   $$\lim_{r\rightarrow0}u(r,\theta)=u_{0}$$

La première condition traduit une périodicité en $\theta$ de période $2\pi$ . Elle se traduit au niveau discret par le fait que les équations pour la ligne $\theta=0$ correspondent aux inconnues $u_{i,1}$ et font intervenir des valeurs en $\theta=-d\theta$ qui ne sont pas définies, i.e. $u_{i,0}^{n}=u(idr,-d\theta)$ . De même les équations pour la dernière ligne $\theta=2\pi-d\theta$ correspondent aux inconnues $u_{i,N_{\theta}}$ et font intervenir des valeurs en $\theta=2\pi$ qui ne sont pas définies, i.e. $u_{i,N_{\theta}+1}^{n}=u(idr,2\pi)$ . Les conditions de périodicité permettent d'imposer :

$$u_{i,0}^{n}=u_{i,N_{\theta}}^{n} u_{i,N_{\theta}+1}^{n}=u_{i,1}^{n}$$ (4.18)

dans les équations pour $j=1$ et $j=N_{\theta}$ .

La seconde condition est nécessaire, car l'équation discrétisée (4.16) dégénère en $r=0$ , de même que l'équation exacte (4.12) à cause des termes en $1/r$ .

Figure 4.11: conditions aux limites en $r=0$

Pour lever cette dégénérescence en $r=0$ , on utilise l'équation discrétisée en coordonnées cartésiennes (4.15). Avec les notations de la figure (4.11) et suivant les axes $(x,y)$ , cette équation s'écrit en $r=0$ :

$$\frac{u_{0}^{n+2}-2u_{0}^{n}+u_{0}^{n-1}}{dt^... ...5}^{n}}{dr^{2}}+\frac{u_{1,3}^{n}-2u_{0}^{n}+u_{1,7}^{n}}{dr^{2}}\right)$$

En effectuant une rotation des axes $(x,y)$ de 45 degrés, on obtiens une autre équation équivalente:

$$\frac{u_{0}^{n+2}-2u_{0}^{n}+u_{0}^{n-1}}{dt^... ...6}^{n}}{dr^{2}}+\frac{u_{1,4}^{n}-2u_{0}^{n}+u_{1,8}^{n}}{dr^{2}}\right)$$

La valeur de $u_{0}$ étant unique, on choisit la moyenne de ces équations:

$$\frac{u_{0}^{n+2}-2u_{0}^{n}+u_{0}^{n-1}}{dt^... ...ft(\frac{1}{8}\sum_{j=1}^{8}u_{1,j}^{n}\right)-u_{0}^{n}}{dr^{2}}\right)$$

soit, de façon générale si on a $N_{\theta}$ noeuds dans la direction $\theta$ :

$$\frac{u_{0}^{n+2}-2u_{0}^{n}+u_{0}^{n-1}}{dt^{2}}=4c_{0}^{2}\left... ...{N_{\theta}}\sum_{j=1}^{N_{\theta}}u_{1,j}^{n}\right)-u_{0}^{n}}{dr^{2}}\right)$$ (4.19)

Cette dernière équation peut s'interpréter comme un bilan de flux sur un disque de rayon $dr/2$ . En intégrant l'équation (4.12) sur ce disque, il vient, après utilisation du théorème de Green:

$$\int_{0}^{\frac{dr}{2}}\frac{\partial^{2}u}{\... ...i}\left(\frac{\partial u}{\partial r}\right)_{dr/2} \frac{dr}{2}d\theta$$

On approxime chacun de ses termes par différences finies. Pour le premier terme, on utilise l'approximation de la dérivée seconde en temps en $r=0$ :

$$\int_{0}^{\frac{dr}{2}}\frac{\partial^{2}u}{\... ...prox\frac{u_{0}^{n+2}-2u_{0}^{n}+u_{0}^{n-1}}{dt^{2}}\pi\frac{dr^{2}}{4}$$

et pour le second l'approximation de la dérivée première en $dr/2$ sur chaque rayon d'angle $\theta_{j}$

$$\left(\frac{\partial u}{\partial r}\right)_{dr/2}\approx\frac{u_{1,j}-u_{0}}{dr}$$

ce qui fournit l'approximation de la seconde intégrale:

$$\int_{0}^{2\pi}\left(\frac{\partial u}{\parti... ...^{N_{\theta}}\frac{u_{1,j}-u_{0}}{dr}\frac{dr}{2}\frac{2\pi}{N_{\theta}}$$

En combinant ces deux approximations, on retrouve l'équation (4.19), qui permet de calculer l'évolution temporelle de la valeur $u_{0}$ .

4.3.3.0.2 condition initiale

La résolution numérique du schéma explicite (4.16) nécessite l'initialisation de la solution $u^{0}$ à $t=0$ et $u^{1}$ à $t=dt$ . On applique la même démarche que pour l'équation des ondes en 1D (paragraphe c3ondes). Disposant des deux conditions initiales du problème (4.12), la valeur $u^{0}$ est donnée par la première condition:

$$u_{i,j}^{0}=w(idr,jd\theta)$$ (4.20)

et la valeur de $u^{1}$ est obtenue à partir d'un développement limité en temps à l'ordre 2 autour de $u^{0}$ :

$$u_{i,j}^{1}=u_{i,j}^{0}+\left(\frac{\partial ... ...+\left(\frac{\partial^{2}u}{\partial t^{2}}\right)_{t=0}\frac{dt^{2}}{2}$$

La valeur de $\frac{\partial u}{\partial t}$ est fournie par la seconde condition initiale, et on utilise l'équation exacte pour calculer $\frac{\partial^{2}u}{\partial t^{2}}$ en fonction du laplacien de $w$ , que l'on discrétise ensuite par différences finies centrées:

$$\left(\frac{\partial^{2}u}{\partial t^{2}}\right)_{t=0}=c_{0}^{2}\Delta w$$

On obtiens ainsi la valeur de $u^{1}$ avec une précision $O(dt^{2},dx^{2})$ , identique à celle du schéma:

$$u_{i,j}^{1} = w_{i,j}+\frac{c_{0}^{2}dt^{2}}{2}\left(\frac{w_{i+1,j}-2w_{i,j}+w_{i-1,j}}{dr^{2}}+\frac{1}{r_{i}}\frac{w_{i+1,j}-w_{i-1,j}}{2dr}\right.$$ (4.21)

$$+\left.\frac{1}{r_{i}^{2}}\frac{w_{i,j+1}-2w_{i,j}+w_{i,j-1}}{d\theta^{2}}\right)$$

4.3.3.1 Stabilité et précision du schéma

L'étude de la stabilité et de la consistance est effectuée tout d'abord sur l'équation discrétisée en coordonnées cartésiennes (4.15), et nous en déduirons ensuite les propriétés pour l'équation discrétisée en coordonnées polaires (4.16).

4.3.3.1.1 Étude de la stabilité:

L'étude de la stabilité utilise le programme Maple 4.3.3. On remplace (ligne 12) dans l'équation discrétisée définie à la ligne 7, la solution approchée $u_{i,j}^{n}$ par une perturbation, que l'on a décomposé en mode de Fourier suivant $x$ et $y$ (ligne 11):

$$Up_{i,j}^{n}=\psi^{n}e^{I\omega_{1}idx}e^{I\omega_{2}jdy}$$

Après simplification, on obtiens une équation du second degré pour le facteur d'amplification $G=\frac{\psi^{n+1}}{\psi^{n}}$ (ligne 19):

$$G^{2}+2bG+1=0$$ (4.22)

Le coefficient $b$ est simplifié (lignes 21 et 22), et s'exprime en fonction de 2 nombres de Courant $C_{FL_{1}}=\frac{c_{0}dt}{dx}$ et $C_{FL_{2}}=\frac{c_{0}dt}{dy}$ et de $y_{1}=\sin\omega_{1}\frac{dx}{2}$ et $y_{2}=\sin\omega_{2}\frac{dy}{2}$ :

$$b=-1+2C_{FL_{1}}^{2}y_{1}^{2}+2C_{FL_{2}}^{2}y_{2}^{2}$$

Le produit des racines de l'équation (4.22) étant égale à 1, la condition de stabilité $\vert G\vert\le1$ impose donc que ces racines soient complexes conjuguées, i.e. que le discriminant soit négatif:

$$\Delta=b^{2}-1\le0 \forall y_{1},y_{2}\in[-1,1]$$

ce qui conduit à la condition (ligne 25 et 26):

$$C_{FL_{1}}^{2}y_{1}^{2}+C_{FL_{2}}^{2}y_{2}^{2}\le1 \forall y_{1},y_{2}\in[-1,1]$$

La condition de stabilité du schéma explicite (4.15) s'écrit donc:

$$C_{FL_{1}}^{2}+C_{FL_{2}}^{2}\le1$$

C'est une condition de Courant:

$$C_{FL}=\frac{c_{0}dt}{h}\le1$$ (4.23)

basé sur une longueur caractéristique $h$ de la maille différence finie définie par:

$$\frac{1}{h}=\frac{1}{dx}+\frac{1}{dy} \leadsto h=\frac{dxdy}{\sqrt{dx^{2}+dy^{2}}}$$ (4.24)

Si les pas de discrétisation sont égaux ( $dx=dy$ ), on a $h=\frac{dx}{\sqrt{2}}$ . La condition de stabilité est donc plus sévère en 2D qu'en 1D.

programme Maple 4.3.3: Etude de stabilité et de consistance de l'équation des ondes

> restart;with(plots):
# Equation des ondes en cartesien
> diff(U(x,y,t),t$2)=c0^2*(diff(U(x,y,t),x$2)+
                           diff(U(x,y,t),y$2));
> eq:=%:
# Equation D.F.
> (U[i,j,n+1]-2*U[i,j,n]+U[i,j,n-1])/dt^2=c0^2*((U[i+1,j,n]-
   2*U[i,j,n]+U[i-1,j,n])/dx^2+(U[i,j+1,n]-
   2*U[i,j,n]+U[i,j-1,n])/dy^2); eqh:=%:
# Etude de stabilite
> Up:=(i,j,n)->Psi[n]*exp(I*omega[1]*i*dx)*exp(I*omega[2]*j*dy);
> subs(U[i,j,n+1]=Up(i,j,n+1),U[i,j,n-1]=Up(i,j,n-1),
  U[i,j,n]=Up(i,j,n),U[i-1,j,n]=Up(i-1,j,n),
  U[i+1,j,n]=Up(i+1,j,n),U[i,j-1,n]=Up(i,j-1,n),
  U[i,j+1,n]=Up(i,j+1,n),eqh);
#
> rel1:=simplify(%*exp(-I*omega[1]*i*dx)*exp(-I*omega[2]*j*dy)):
> simplify(subs(Psi[n+1]=G*Psi[n],Psi[n-1]=Psi[n]/G,rel1*G/Psi[n]));
> collect(dt^2*(lhs(%)-rhs(%)),G)=0;eq3:=lhs(%):
# calcul du coefficient de G
> coeff(eq3,G):expand(%);
> subs(cos(omega[1]*dx)=1-2*y[1]^2,cos(omega[2]*dy)=1-2*y[2]^2,
       dx=dt*c0/CFL[1],dy=dt*c0/CFL[2],%):simplify(%);b:=%/2;
# Racines G complexes conjuguees si  Delta<0
> Delta:=b^2-1;factor(%);
> CFL[1]^2*y[1]^2+CFL[2]^2*y[2]^2<1;
# Condition de stabilite
> CFL[1]^2+CFL[2]^2<1;cdts:=%:
> subs(CFL[1]=c0*dt/dx,CFL[2]=c0*dt/dy,cdts);
> dt*c0*sqrt(1/dx^2+1/dy^2)<1;
# Erreur de troncature
> Uex:=(p,q,r)->U(x+(p-i)*dx,y+(q-j)*dy,t+(r-n)*dt);
> subs(U[i,j,n]=Uex(i,j,n),U[i,j,n-1]=Uex(i,j,n-1),U[i,j,n+1]=
   Uex(i,j,n+1),U[i+1,j,n]=Uex(i+1,j,n),U[i-1,j,n]=Uex(i-1,j,n),
   U[i,j+1,n]=Uex(i,j+1,n),U[i,j-1,n]=Uex(i,j-1,n),
   lhs(eqh)-rhs(eqh)); rel3:=%:
#
> U(x,y,t+dt)=convert(mtaylor(U(x,y,t+dt),[dt],8),diff);S1:=%:
> U(x,y,t-dt)=convert(mtaylor(U(x,y,t-dt),[dt],8),diff):S2:=%:
> U(x+dx,y,t)=convert(mtaylor(U(x+dx,y,t),[dx],8),diff):S3:=%:
> U(x-dx,y,t)=convert(mtaylor(U(x-dx,y,t),[dx],8),diff):S4:=%:
> U(x,y+dy,t)=convert(mtaylor(U(x,y+dy,t),[dy],8),diff):S5:=%:
> U(x,y-dy,t)=convert(mtaylor(U(x,y-dy,t),[dy],8),diff):S6:=%:
> subs(S1,S2,S3,S4,S5,S6,rel3-(lhs(eq)-rhs(eq))):simplify(%);
# Schema d'ordre 2

Par analogie, la condition de stabilité du schéma explicite (4.16) en coordonnée polaire est aussi la condition de Courant (4.23). Il faut cependant définir la longueur caractéristique $h$ de la maille en coordonnée polaire. Cette longueur ne doit pas être basée sur les dimensions $dr$ et $d\theta$ de la maille dans l'espace transformé, mais sur les dimensions de la maille $dr$ et $rd\theta$ dans l'espace physique:

$$h=\frac{rdrd\theta}{\sqrt{dr^{2}+(rd\theta)^{2}}}$$

La maille dans l'espace physique n'étant pas constante (elle dépend de $r$ ), on doit prendre la plus petite valeur de $h$ dans le maillage, qui est obtenue pour $r=dr$

$$h_{min}=\frac{dr^{2}d\theta}{\sqrt{dr^{2}+(drd\theta)^{2}}}\approx drd\theta$$

ce qui fournit la condition de stabilité du schéma (4.16):

$$C_{FL}=\frac{c_{0}dt}{drd\theta}\le1$$ (4.25)

On constate que cette condition en coordonnée polaire est beaucoup plus contraignante que la condition (4.23) en coordonnées cartésiennes. Pour un maillage de $N$ points dans chaque direction, il faut choisir un pas en temps $dt$ tel que $dt\le c_{0}\sqrt{2}/N$ en cartésien et $dt\le c_{0}/N^{2}$ en polaire.

4.3.3.1.2 Erreur de troncature:

Le calcul Maple (lignes 32 à 44 du programme 4.3.3) fournit l'erreur de troncature pour le schéma (4.15) en cartésien:

$$E_{t}=\frac{1}{12}\frac{\partial^{4}u}{\parti... ...}}{12}\frac{\partial^{4}u}{\partial y^{4}}dy^{2}+O(dt^{4},dx^{4},dy^{4})$$

Le schéma explicite (4.15) est donc d'ordre 2 en temps et en espace (i.e. a une précision en $O(dt^{2},dx^{2},dy^{2})$ ).

Par analogie, on en déduit que le schéma explicite (4.16) est aussi d'ordre 2 en temps et en espace (i.e. a une précision en $O(dt^{2},dr^{2},d\theta^{2})$ ).

Enfin on note que ces schémas explicites (4.15) et (4.16) sont des schémas dispersifs, mais non dissipatifs comme en 1D (puisque le facteur d'amplification $G$ à un module égal à 1).

4.3.4 Expérimentation numérique avec Matlab

Le programme Matlab (4.3.4) implémente le schéma explicite en cordonnées polaires (4.16). Les paramètres du calcul sont définies sur les lignes 3 à 13, avec un $C_{FL}$ définit par la relation (4.25). Pour tenir compte simplement des conditions de périodicité, la solution approchée $u_{i,j}^{n}$ est définie comme un tableau de dimension $N_{r}*(N_{\theta}+1)$ , i.e. avec $r$ variant de 0 ( $i=1)$ à $1$ $(i=N_{r})$ et $\theta$ de 0 $(j=1)$ à $2\pi+d\theta$ $(j=N_{\theta}+1)$ . La solution initiale $w$ est définie sur les lignes 14 à 21 comme combinaison linéaire des 2 modes propres (4.15). Cette solution initiale permet l'initialisation des champs $Un_{0}$ (lignes 24 à 26). La seconde condition initiale (4.21) est appliquée lors de la première itération en temps en notant que cette condition est dans notre cas équivalente à l'équation aux différences (4.16) pour $n=1$ avec $u_{i,j}^{-1}=u_{i,j}^{1}$ , i.e. l'équivalent d'une condition miroir à l'instant initiale. Cette condition peut alors être implémentée en initialisant $u^{0}$ et $u^{-1}$ avec $w$ , et en calculant $u^{1}$ à la première itération avec la formule générale (4.16), et un coefficient $c_{0}^{2}$ divisé par deux (ligne 37).

Dans les itérations en temps (lignes 34 à 57), on utilise l'équation aux différences (4.16) pour les noeuds internes $1<i<N$ , $1<j<N_{\theta}$ en utilisant la programmation matricielle Matlab (ligne 40) . La condition aux limites en $r=1$ fournit la valeur aux noeuds frontières $i=N$ (ligne 45). Les conditions de périodicité fournissent les valeurs aux noeuds frontières $j=1$ et $j=N_{\theta}+1$ (lignes 49 à 50). Enfin l'équation pour les noeuds en $r=0$ est écrite aux lignes 52 à 54.

programme Matlab 4.3.4: Résolution numérique de l'équation des ondes (4.11)

% equation des ondes en polaire
% schema explicite
clear;
R1=1;
Ntheta=50; Nr=80;
dtheta=2*pi/(Ntheta-1);dr=R1/(Nr-1);
R=[0:dr:R1]';
Theta=[0:dtheta:2*pi];
% pts du maillage
X=R*cos(Theta); Y=R*sin(Theta);
% parametre
c0=1; CFL=0.9;
dt=CFL*dr*dtheta/c0
%  cdts initial
k1=0; a1=1.0;
lambda1=fzero(inline('0*besselj(0,r)-r*besselj(1,r)'),10)
BJ1=inline('besselj(0,r)','r');
k2=1; a2=0.5;
lambda2=fzero(inline('1*besselj(1,r)-r*besselj(2,r)'),10)
BJ2=inline('besselj(1,r)','r');
W=(a1*BJ1(lambda1*R)*cos(k1*Theta)+...
   a2*BJ2(lambda2*R)*cos(k2*Theta));
% initialisation
Un0=zeros(Nr,Ntheta+1);
Un0(:,1:Ntheta)=W; Un0(:,Ntheta+1)=Un0(:,2);
Un=Un0;Un1=Un;
% noeuds internes
I=[2:Nr-1];J=[2:Ntheta];
RI=R(I)*ones(1,Ntheta-1);
% schema D.F.
Tf=2*(2*pi/c0/lambda1); 
nit=round(Tf/dt)
% iteration
for it=1:nit
   coef=c0^2*dt^2;
   if (it==1) coef=coef/2; end;
   % noeuds internes
   Un1(I,J)=2*Un(I,J)-Un0(I,J) + ...
     (coef/dtheta^2)*(Un(I,J+1)-2*Un(I,J)+Un(I,J-1))./(RI.^2)+ ...
     (coef/(2*dr))*(Un(I+1,J)-Un(I-1,J))./RI + ...
     (coef/dr^2)*(Un(I+1,J)-2*Un(I,J)+Un(I-1,J));
  % C.L. en r=1
   Un1(Nr,J)=2*Un(Nr,J)-Un0(Nr,J) + ...
     (coef/dtheta^2)*(Un(Nr,J+1)-2*Un(Nr,J)+Un(Nr,J-1))/(R(Nr))+ ...
     (coef/dr^2)*(Un(Nr-1,J)-2*Un(Nr,J)+Un(Nr-1,J));
   % periodicite
   Un1(1:Nr,1)=Un1(1:Nr,Ntheta);
   Un1(1:Nr,Ntheta+1)=Un1(1:Nr,2);
   % C.L. en r=0
   Um=sum(Un(2,2:Ntheta))/(Ntheta-1);
   Um1=2*Un(1,1)-Un0(1,1)+(4*coef/dr^2)*(Um-Un(1,1));
   Un1(1,1:Ntheta+1)=Um1;
   % iteration suivante
   Un0=Un; Un=Un1;
end;

Sur la figure (4.12), on a tracé la solution calculée au bout d'une période, avec $C_{FL}=0.9$ , $N_{r}=80$ et $N_{\theta}=50$ , pour les deux conditions initiales étudiées analytiquement au paragraphe sol_ondes2d. Elles se comparent parfaitement aux solutions analytiques de la figure (4.9).

Figure 4.12: solutions numériques de l'équation des ondes

[k=0]

[k=1]

Pour tester la stabilité du schéma, nous avons fait varier le pas en temps pour deux maillages donnés: ( $N_{r}=40$ , $N_{\theta}=25$ ) et ( $N_{r}=80$ , $N_{\theta}=50$ ) avec la condition initiale suivante:

$$BesselJ(0,\lambda_{0,3})\cos(\lambda_{0,3}c_{... ...\frac{1}{2}BesselJ(1,\lambda_{1,3})\cos(\theta)\cos(\lambda_{1,3}c_{0}t)$$

Nous avons ensuite calculé l'erreur en $r=0$ en comparant la solution approchée sur l'axe $u_{1,1}$ et la solution exacte $u_{ex}=\cos(c_{0}\lambda_{0,3}t$ ) sur un temps $\tau$ de l'ordre de deux périodes $\tau\approx1.4$ . On a tracé ces évolutions sur la figure (4.13). On constate bien que la solution diverge dès que le nombre de Courant (4.25) est supérieur ou égale à 1. On peut aussi noter que si l'on choisit une condition initiale ne dépendant pas de $\theta,$ la solution reste stable pour des $C_{FL}$ beaucoup plus grands (i.e. $C_{FL}\le N_{\theta}$ ), ce qui montre que l'instabilité la plus sévére proviens de la discrétisation du terme en $\theta$ . Théoriquement au bout d'un nombre très grand d'itérations, les erreurs d'arrondis devraient pouvoir déstabiliser la solution, mais ici la symétrie des calculs fait que ces erreurs sont indépendantes de $\theta$ et la solution non perturbée reste stable.

Figure 4.13: évolution de l'erreur sur l'axe en fonction du $C_{FL}$

[ $N_r=40, N_\theta =25$ ]

[ $N_r=80, N_\theta =50$ ]

Pour étudier la précision du calcul, nous avons calculer l'erreur au centre pour différents maillages avec un $C_{FL}=0.9$ avec la même condition initiale. La solution au centre étant indépendante de $\theta$ , nous avons uniquement fait varier la discrétisation suivant $r$ en choisissant des valeurs de $N_{r}$ de 10 à 640. La taille caractéristique du maillage est $h\approx\frac{1}{N_{r}}$ , et nous avons tracé l'erreur sur l'axe en fonction de $h$ sur la figure (4.14).

Figure 4.14: erreur sur l'axe en fonction de la taille du maillage

On constate sur cette figure que l'erreur se comporte à la limite en $O(h^{2})$ , ce qui était prévue par la théorie. Cela montre que notre condition en $r=0$ préserve la précision d'ordre 2 du schéma.