🐡
MATLABで学ぶ音響工学 球面波2
前回の記事に引き続き、今回は2次元でのプログラム方法について説明致します。
ちなみに、前回の記事で1次元での球面波のプログラム方法を紹介しています。
前回の記事の繰り返しになりますが、下記の物性記号を用いて、球面波の圧力pは下式で表せます。
ρ：空気の密度
j：√-1
ω：角周波数
k：音波の波数（2π／波長）
上式の距離rをxとyの直交座標系で表すと下のようになります。
プログラム上ではxとyという2変数に対して、rが決まり、rが決まると圧力pが決まるという流れになります。そして、そのpを行列の要素内に保存して、全てのxとyの組み合わせが計算し終わったら、表示をするという流れです。
例えば、xを列、yを行とする行列Pの中に各pを格納していきます。
プログラムをみた方が簡単に理解できると思うので、下のプログラムをご確認ください。
clear all;close all;clc
x=-2:0.1:2;
y=-2:0.1:2;
rho=1.293;
A=1;
c=340; %音速
freq=100; %周波数
w=2*pi*freq;
k=w/c;
t=(0:0.1:2*pi)/w;
P=zeros(length(y),length(x));
for ii1=1:length(x)
for ii2=1:length(y)
r=sqrt(x(ii1)^2+y(ii2)^2);
if r<0.1
r=0.1; %r=0の時にP(ii1,ii2)が無限大になるのを防ぐ目的
end
P(ii2,ii1)=j*w*rho*A./(4*pi*r).*exp(j*(k*r));
end
end
% x(ii1)およびy(ii2)という位置における振幅abs(P(ii1,ii2))は一意に決まるため
% P(ii1,ii2)を複素数平面上で回転させれば、音波が伝搬するアニメーションが作成できる
% figure(1)
% surf(x,y,abs(P))
% shading interp
% ylabel('X [m]')
% xlabel('Y [m]')
t=(0:0.3:2*pi)/w;
for ii1=1:length(t)
clf
surf(x,y,real( P*exp(j*(-w*t(ii1))) ))
shading interp
ylabel('X [m]')
xlabel('Y [m]')
view([0 90])
% pause(0.1)
saveas(gcf,[num2str(ii1+100) '.png'],'png')
end
1.5MBまでしかアップロードできないので、粗いアニメーションになっていますがご了承ください。
今回はこのへんでGood luck
私の個人ブログもよろしくお願いします。
Discussion