Mise en œuvre de la réalisation diffusive de l’opérateur Z

 Les simulations ont été effectuées sur un maillage de 128 (r) ×1024 (θ) points dans le disque unité privé d’un disque centré à l’origine, de rayon 0.1 avec conditions de Dirichlet nulles. La condition d’impédance est quand à elle discrétisée au moyen de 20 points en la variable ξ. La transformation de Fourier sur le bord transparent est la FFT, fonction disponible en standard sous Matlab. Les résultats de simulation confirment dans tous les cas la parfaite absorption de la frontière, tant pour les incidences orthogonales aux bord qu’obliques ou rasantes. Il faut noter que des fluctuations trop rapides engendrent des erreurs numériques qui se propagent et dégradent la solution: il faut en pratique respecter un minimum de 10 points par période pour une onde sinusoïdale. En Figure 1 est visible l’évolution d’un ébranlement provoqué par une source ponctuelle. La neutralité du bord aux ondes rasantes est clairement mise en évidence. Après passage de l’ébranlement, le milieu reste au repos. Une onde rampante de la forme ϕ = Re ¡ k(r) ei(ωt+nθ) ¢ , provoquée par une source sur le cercle de rayon 0.375 animée d’un mouvement de rotation uniforme, est montrée en Figure 2. N.B.: la transition entre la zone de champ proche (c’est-à-dire lorsque la vitesse de l’onde est inférieure à la célérité du milieu, la décroissance étant alors exponentielle) et la zone de rayonnement par propagation (en spirale, à décroissance en √ 1 r ) est clairement visible. Enfin, on peut voir en Figure 3 l’évolution provoquée par une source ponctuelle constante. En simulant une évolution en espace libre, la condition d’impédance permet ainsi de résoudre au moyen de l’équation des ondes dans un domaine limité un problème d’équilibre statique du type: ∆Φ = f dans R2.

Listing des programmes 

ondepol_imp_RE_stock.m % résolution de l’équation d’ondes 2D pour un domaine circulaire absorbant. % Programme unifié qui permet d’utiliser indifféremment imp pour frontière droite sur un cercle ou % imp adaptée pour le cercle en choisissant à l’appel 1 (pour adaptée) autre chose (pour frontière plane). % En commentaires : possibilité de mettre une source haute fréquence, une Gaussienne, une source temporaire % sur tout le bord interieur, ondes rampantes… % rm et rs sont les rayons min et max du domaine % Nr et Nth sont les nombres de points de discrétisation selon le rayon et l’angle % A0, A1 e A2 sont les matrices Nr*Nth-1 des valeurs de Phi en tn-1, tn et tn+1 % T est le nombre d’itérations en temps, donc à quelque chose près (multiplication par dt) le temps de simulation % X et Y maillage du domaine, M matrice des valeurs de phi. function [X,Y,M]=ondepol_imp_RE_stock(rm,rs,Nr,Nth,T,choix) dim=20; M=zeros(Nr,Nth,20); % vecteur de matrices Nr*Nth qui servira a stocker de temps en temps l’état du système V=zeros(20,1); % vecteur qui servira a stocker les valeurs en un point pour explosion dr=(rs-rm)/(Nr-1); dth=2*pi/(Nth-1); discretr=linspace(rm,rs,Nr); discrett=linspace(0,2*pi,Nth); X=discretr’*cos(discrett); % matrice des abscisses de discrétisation Y=discretr’*sin(discrett); % matrice des ordonnées de discrétisation % Impulsion de départ : Gaussienne %M(:,:,1)=gaussiennepol(X,Y,0,0.8,0.001); % CI pour la simulation de la source ponctuelle M(round(8*Nr/10),1,1)=1; A0=M(:,1:Nth-1,1); A1=A0; A2=zeros(Nr,Nth-1); dt=0.25*min(dr,rm*dth) temps_de_simu=dt*T c1=( (1-dt^2/dr^2)*ones(1,Nr-2)-dt^2./(discretr(2:Nr-1)*dth).^2 )’; c2=(dt^2./(2*discretr(2:Nr-1)*dr))’; c3=(dt^2./(discretr(2:Nr-1)*dth).^2)’; C1=c1*ones(1,Nth-3); C2=c2*ones(1,Nth-3); C3=c3*ones(1,Nth-3); mu=4; % fréquence de la source omega=2*pi*mu; % pulsation de la source psi=zeros(dim,Nth-2); % une colonne de psi est le vecteur à k composante psi_n, n balayant les fréquences. psi0=zeros(dim,1); fftphi=zeros(Nth-1,1); % vecteur sortie de la RE (vecteur Fourrier(phi)) F=zeros(dim,round((Nth-1)/2));