MATLABで音響FDTDプログラムを作成しました.とりあえず1次元バージョンです.
*自由音場における音圧p,粒子速度uに関する偏微分方程式を差分化して,逐一音圧と粒子速度を交互に計算していきます.
*密度と体積弾性率は室温20度の値を用いています.
*Murの1次吸収境界を用いています.
*初期波源は音圧に1次元ガウシアンパルスを与え,粒子速度をゼロとしています.
*自由音場における音圧p,粒子速度uに関する偏微分方程式を差分化して,逐一音圧と粒子速度を交互に計算していきます.
*密度と体積弾性率は室温20度の値を用いています.
*Murの1次吸収境界を用いています.
*初期波源は音圧に1次元ガウシアンパルスを与え,粒子速度をゼロとしています.
%===Parameters===% NT=200;NX=101; SIG=10e-3/sqrt(2); RHO=1.21;KKK=1.42e5; cc=sqrt(KKK/RHO); dx=1e-3;dt=1e-6; mur=(cc*dt-dx)/(cc*dt+dx); duu=dt/RHO/dx; dpp=KKK*dt/dx; x=0:dx:(NX-1)*dx; p=exp(-((x-0.05).^2)/2/SIG/SIG); u=zeros(size(x)); for T=0:NT; t=T*dt; %===update u===% u(1:NX-1)=u(1:NX-1)-duu*(p(2:NX)-p(1:NX-1)); %===update p===% tmp1=p(2);tmp2=p(NX-1); p(2:NX-1)=p(2:NX-1)-dpp*(u(2:NX-1)-u(1:NX-2)); %===Mur1st===% p(1)=tmp1-mur*(p(1)-p(2)); p(NX)=tmp2-mur*(p(NX)-p(NX-1)); %===Plot===% plot(x,p); axis([0,max(x),-0.2,1.2]); drawnow; end