放射音場計算

無限バッフルに囲まれた円盤(円形)トランスデューサからの放射音場計算

円盤トランスデューサからの放射音場計算

MATLABで書いたプログラムその3です。 無限バッフル(剛体)に囲まれた、円盤(円形)トランスデューサから、三次元空間に放射された音場情報の計算をします。 結果は円盤トランスデューサと平行な二次元平面における複素音圧分布を出力します。 比較的計算時間が早いのが特徴です。

計算手法としては、幾何的な手法であると言えます。 詳細は以下の文献を参照してください。 この手法は、円盤トランスデューサだけでなく、矩形トランスデューサの計算にも応用可能です。

P. R. Stepanishen,"Transient Radiation from Pistons in an Infinite Planar Baffle," J. Acoust. Soc. Am. 49 (1971) pp. 1629-1638

原理の概略

円盤トランスデューサは一様に正弦振幅をしているとします。 しかし、まずはインパルス応答を計算します。 観測面のある一点、観測点でのインパルス応答を計算し、時間応答をシミュレーションします。 このインパルス応答は、円盤トランスデューサと観測点との位置関係で幾何的に計算できます。

このインパルス応答と、正弦振幅による調和振動応答との関係はフーリエ変換になっています。 インパルス応答をフーリエ変換し、円盤トランスデューサの駆動周波数におけるフーリエ係数に密度等の係数をかけたものが調和振動における音場を示します。 インパルス応答とは違い、ある値となります。

この観測点における音場計算を、繰り返し行うことで、二次元音圧分布が得られます。

これまで、我々の研究室ではこの原理に基づき音場計算をするプログラムが使われてきましたが、計算時間が数時間もかかるという問題がありました。

高速化への工夫

そこで、いくつかの解決策を考えました。 まず、使われてきた音場計算プログラムは、インパルス応答のフーリエ変換をFFTにより全周波数に渡って行っていました。 必要であるのは、超音波トランスデューサの駆動周波数だけですので、特定の周波数における積分のみで十分です。 これにより、計算コストが大幅に縮小されました。

次に、その積分にあたり、単純に区分求積法を適応すると十分な計算精度を得るためにインパルス応答を計算する際の時間刻みが十分小さくなければなりません。 そこで、シンプソンの公式を用いることで、時間刻みがそれほど小さくなくても十分な精度を得ることが出来ました。

最後に、円盤トランスデューサによる調和振動応答ですので、超音波トランスデューサからの距離で音場は決まります。 従って、観測面における音軸上の点からの距離に対する計算を行い、二次元音場を計算せずに一次元の音場を計算し、その後に内挿により二次元に戻しました。

以上の工夫をし、matlabで高速化をした結果、音場計算が数秒で計算可能になりました。 このプログラムの長所は前述しましたが、計算速度と計算精度です。

Matlabによる実装

図に示すような座標系を考え、matlabで実装しました。

%円盤トランスデューサから放射された二次元複素音場を求める
clear

%設定が必要な定数
C=1500;%音速(m/s)
Rho=1000;%密度(kg/m3)
D=4e-3;%超音波トランスデューサの半径(m)
F=1.43e6;%超音波の周波数(Hz)
Z=0;%観測面と振動面の距離(m)
r_max=12e-3;%測定半径(m)
x_step=1e-4;%x,y軸の刻み幅(m)
nk=100;%シンプソン積分において、nk等分する(偶数である必要がある)

%変数
r_step=x_step*1e-2;%r軸刻み幅(m) 必要刻み幅の100倍としている
ko=1:2:nk-1;%シンプソン積分のおける奇数番
ke=2:2:nk-2;%シンプソン積分のおける偶数番
p=zeros(size((-r_max:x_step:r_max)',1),size((-r_max:x_step:r_max)',1));%二次元複素音場
int_h=zeros(size((0:r_step:r_max * 1.5)',1),1);%インパルス応答を積分したもの
h=zeros(nk,1);%r1<ct<=r2のインパルス応答

% r=0の時のインパルス(Integrate[c Exp[I w t], {t, z/c, Sqrt[z^2 + d^2]/c}])
jw = j*2*pi*F;
int_h(1) = C * (exp(jw * sqrt(D^2 + Z^2) / C) - exp(jw * Z / C)) / jw;

%0<r<=Dのインパルス
for r=r_step:r_step:r_max * 1.5 % sqrt(2) では1点足りなくなるかもしれない
    r1=(Z^2+(r-D).^2).^0.5;
    r2=(Z^2+(r+D).^2).^0.5;
    rk = round(r ./ r_step) + 1;
    
    %ct=r1の値(r>=Dの場合におけるインパルス応答の立ち上がり)
    ct = r1;
    t = ct ./ C;      
    ctz2 = ct.^2 - Z^2;   
    h0 = (C * exp(jw .* t)).*(r1==Z)+(acos((ctz2+r.^2-D^2)./(2*r.*ctz2.^0.5))*C/pi .* exp(jw .* t)).*(r1~=Z);
    
    %ct=r1以降、インパルス応答が0になるまで
    ctstep = (r2 - r1) / nk;
    ct = r1 + ctstep : ctstep : r2;
    t = ct ./ C;
    ctz2 = ct.^2 - Z^2;
    h = acos((ctz2+r.^2-D^2)./(2*r.*ctz2.^0.5))*C/pi .* exp(jw .* t);
    
    % シンプソンの積分公式
    summ = h0 +2*sum(h(ke))+4*sum(h(ko))+ h(nk);
    int_h(rk) = (C * (exp(jw * r1 / C) - exp(jw * Z / C)) / jw).*(r<D) + summ * ctstep / C / 3;
end

%得られたH(r)から音圧分布p_r(r)を求める
p_r=-j*Rho*2*pi*F*int_h;

%二次元音場(complex)の計算
[x y]=meshgrid(-r_max:x_step:r_max, -r_max:x_step:r_max);
r=sqrt(x.^2+y.^2);
p=reshape(p_r(round(r./r_step)+1), size(x));
最終更新 May 14, 2020: ブログ更新 (e287c4c)