Masassiah Blog

Masassiah のブログです。主に読書で得た気づきをまとめています。

音響FDTD2次元

イメージ 1

以前に書いた音響FDTD1次元バージョンを拡張して、同じくMATLABで音響FDTD(Finite Difference Time Domain)2次元バージョンのプログラムを作ってみた。
音圧と粒子速度に関する運動方程式と連続の式を、空間、時間領域で差分方程式に展開して、逐次計算しています。ある時刻の音圧絶対値の空間分布の画像を1枚だけ掲載していますが、プログラムでは2次元ガウス分布が自由空間に広がっていくアニメーションを表示するようにしています。

MATLABプログラム
% acoustic FDTD 2D ver1.00
clear all;
close all;
%===Parameter===%
rho=1.21;K=1.42e5;c=sqrt(K/rho);
d=2e-3;NX=101;NY=101;
CFL=0.5;
dt=CFL*d/c;NT=200;
x=0:d:d*(NX-1); y=0:d:d*(NY-1);
[Y,X]=meshgrid(y,x);
%===FDTD's coefficient===%
pu=dt/rho/d;
pp=K*dt/d;
pm=(c*dt-d)/(c*dt+d);
%===Memory occupation===%
p=zeros(NY,NX);
ux=zeros(NY,NX);uy=zeros(NY,NX);
%===Initial condition===%
sigma=20e-3;
r=(X-0.1).^2+(Y-0.1).^2;
p=exp(-(r/sigma^2));
%===Stat Time Step===%
for T=0:NT
    %===Update u===%
    uy(2:NY-1,2:NX-1)=uy(2:NY-1,2:NX-1)-pu*(p(3:NY,2:NX-1)-p(2:NY-1,2:NX-1));
    ux(2:NY-1,2:NX-1)=ux(2:NX-1,2:NX-1)-pu*(p(2:NY-1,3:NX)-p(2:NY-1,2:NX-1));
    %===tmp for mur1st===%
    p(1,:)=p(3,:); p(NY,:)=p(NY-2,:);
    p(:,1)=p(:,3); p(:,NX)=p(:,NX-2);
    %===Update p===%
    p(3:NY-2,3:NX-2)=p(3:NY-2,3:NX-2)-pp*(...
        uy(3:NY-2,3:NX-2)-uy(2:NY-3,3:NX-2)...
       +ux(3:NY-2,3:NX-2)-ux(3:NY-2,2:NX-3));
    %===mur1st===%
    p(2,:)=p(1,:)+pm*(p(3,:)-p(2,:));
    p(NY-1,:)=p(NY,:)+pm*(p(NY-2,:)-p(NY-1,:));
    p(:,2)=p(:,1)+pm*(p(:,3)-p(:,2));
    p(:,NX-1)=p(:,NX)+pm*(p(:,NX-2)-p(:,NX-1));
    %===Figure===%
    imagesc(x,y,abs(p));
    axis('equal');
    caxis([0,0.25]);
    xlim([0,max(x)]);ylim([0,max(y)]);
    grid on;
    drawnow
end