Programmation en Matlab

Simulation de la propagation d'une onde 2D

Quelques éléments de programmation pour le TP 9: Simulation bidimensionnelle de la propagation d'une onde.

Initialisation: paramètres et constantes physiques

clear all;close all % Pour bien commencer...
tic % Si on veut mesurer par la suite le temps d'exécution (avec alors toc)
c=10; % célérité de l'onde

% Domaine et grille de calcul
Lx=100;Ly=100;
T=10;

deltax=Lx/200;
deltay=Ly/200;
deltat=sqrt(deltax^2+deltay^2)/(2*c);

gax=(c*deltat)^2/(deltax^2);
gay=(c*deltat)^2/(deltay^2);

% Grille de calcul et nombre de points
xx=[0:deltax:Lx];nx=length(xx);
yy=[0:deltay:Ly];ny=length(yy);
tt=[0:deltat:T];nt=length(tt);

% Fréquence: 2Hz, de la source
nu=2;

Initialisation de la matrice

Il s'agit de la discrétisation de la relation (3):

% Pour initialiser une matrice en Matlab, on la remplit de zeros
Mat=zeros(nx,ny,nt);

for i=2:nx-1
    for j=2:ny-1
        Mat(i,j,2)=Mat(i,j,1)+...
            gax/2*(Mat(i-1,j,1)+Mat(i+1,j,1)-2*Mat(i,j,1))+...
            gay/2*(Mat(i,j-1,1)+Mat(i,j+1,1)-2*Mat(i,j,1));
    end
end

Equation des ondes discrétisée

Enfin, on s'attaque à la propagation à proprement parler: la discrétisation (2) de l'équation des ondes (1):

for k=3:nt-1
    for i=2:nx-1
        for j=2:ny-1
            Mat(i,j,k+1)=2*Mat(i,j,k)-Mat(i,j,k-1)+...
                gax*(Mat(i-1,j,k)+Mat(i+1,j,k)-2*Mat(i,j,k))+...
                gay*(Mat(i,j-1,k)+Mat(i,j+1,k)-2*Mat(i,j,k));
        end
    end
    % et on impose la source à chaque instant:
      % sx et xy désignent l'abscisse et l'ordonnée dans la grille 
      % de calcul de la source
    Mat(sx,sy,k+1)=sin(2*pi*sfrq*k*deltat);
    % et on impose aussi ici les conditions aux limites: 
      % conditions sur les bords du domaine 
      % et sur l'éventelle obstacle 
      % (fente simple ou double dans le TP)
end

Il ne manque plus donc que les conditions aux limites (CLA: Conditions aux Limites Absorbantes) et sur un éventuel obstacle.

Conditions aux limites

Les conditions aux limites absorbantes (5), sur les bords du domaine, se discrétisent et s'imposent simplement, à chaque instant k, par:

Mat(1,:,k+1)=Mat(2,:,k);
Mat(nx,:,k+1)=Mat(nx-1,:,k);
Mat(:,1,k+1)=Mat(:,2,k);
Mat(:,ny,k+1)=Mat(:,ny-1,k);

Avec ces conditions, tout se passe (presque, cela reste une approximation numérique…) comme si l'onde se propageait en espace libre.
Sans imposer de CLA (en pratique, en n'imposant rien sur les premières et dernières lignes et colonnes de la matrice), l'onde est réfléchie sur les bords: on simule ainsi le comportement d'une onde dans une cavité.

Obstacle

La présence d'un obstacle se simule simplement en imposant, à chaque instant k,

Mat(Obstaclex,Obstacley,k)=0;

où il faut définir les indices dans notre grille de l'obstacle.

Affichage

Pour l'affichage (et uniquement, aucun sens dans la simulation physique), on peut, pour faire apparaître l'obstacle, imposer

Mat(Obstaclex,Obstacley;:)=1;

Toutes les valeurs de la matrice localisées sur l'obstacle auront la valeur maximale et fixée.

Ensuite, il reste à afficher toutes les valeurs de la matrice Mat à chaque instant k, c'est-à-dire Mat(:,:,k).
L'affichage du contenu d'une matrice peut se faire par exemple avec pcolor:

figure(1);clf;whitebg('w');
for k=1:nt
   pcolor(Mat(:,:,k));
    % Les axes n'ont ici aucun sens physique, on les retire
    axis off;
    % On fixe l'échelle des couleurs
    caxis([-0.5 0.5]);
    % On lisse la représentation par interpolation des données
    shading interp;
    % Et on laisse un peu de temps entre chaque affichage
    pause(0.05)
end