MATLABで学ぶ音響工学　球面波2

スポンサーリンク

前回の記事に引き続き、今回は2次元でのプログラム方法について説明致します。

ちなみに、前回の記事で1次元での球面波のプログラム方法を紹介しています。

前回の記事の繰り返しになりますが、下記の物性記号を用いて、球面波の圧力pは下式で表せます。
ρ：空気の密度
j：√-1
ω：角周波数
k：音波の波数（2π／波長）

上式の距離rをxとyの直交座標系で表すと下のようになります。

$$
r = \sqrt{x^2+y^2}
$$

プログラム上では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