MATLAB を使って音波の伝搬を FDTD 法でシミュレーション

MATLAB で音波（sound wave）の伝搬を FDTD（Finite Difference Time Domain）法でシミュレーションしてみた。

はじめに

音波は空気その他の物質，すなわち媒質（medium）中に存在する波動であることから，媒質粒子の振動している状態が媒質の中を伝わっていることになる。

このことから，音波については媒質粒子の振動と，この振動が媒質中を伝搬することとの両方について考えなければならない。

以前に書いた音響 FDTD 1 次元バージョンを拡張して，MATLAB で音響 FDTD 2 次元バージョンのプログラムを作ってみた。

king-masashi.hatenablog.com

FDTD 法による音波の伝搬解析

音圧（sound pressure）と粒子速度（particle velocity）に関する運動方程式と連続の式を，空間，時間領域で差分方程式に展開して，逐次計算している。

ある時刻の音圧絶対値の空間分布の画像を 1 枚だけ掲載しているが，プログラムでは 2 次元ガウス分布が自由空間に広がっていくアニメーションを表示するようにしている。

なお，吸収境界には，Mur の 1 次吸収境界条件（Mur's first order absorbing boundary condition）を採用している。

伝搬解析に用いたパラメータ

音波の伝搬解析に用いたパラメータ(Parameters used in the electromagnetic field analysis.) を示す。

• 空気密度 ρ = 1.21 kg/m³
• 体積弾性率 K = 1.42 × 10⁵ N/m²
• 計算領域 200 mm × 200 mm
• 格子間隔 2 mm
• ステップ数 200
• 波源 (100 mm, 100 mm)

MATLABプログラム

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

参考文献

本稿を作成するにあたり，以下の文献を参考にした。

• 宇野 亨 編著，何 一偉，有馬 卓司 著，「数値電磁界解析のための FDTD 法――基礎と実践――」，”FDTD Method for Computational Electromagnetics -- Fundamentals and Practical Applications --”，コロナ社，2016年5月25日

数値電磁界解析のためのFDTD法- 基礎と実践 -

FDTD法による電磁界およびアンテナ解析 [ 宇野亨 ]

king-masashi.hatenablog.com

更新履歴

• 2008年2月15日　新規作成
• 2022年2月4日　加除修正
• 2022年9月4日　誤植を訂正
• 2022年10月12日　タイトルを「音響FDTD2次元」から「MATLAB で音波の伝搬を FDTD 法でシミュレーション」へ変更
• 2022年12月16日　吸収境界条件の説明を追加，テキストで参考文献を記載，その他全般的に加除修正
• 2023年11月24日　伝搬解析に用いたパラメータを追加
• 2024年1月4日　タイトルを「MATLAB で音波の伝搬を FDTD 法でシミュレーション」から「MATLAB を使って音波の伝搬を FDTD 法でシミュレーション」へ変更
• 2024年1月28日　検索エンジン向けタイトル「音波の伝搬をMATLABでシミュレーションする方法」，SNS 向けタイトル「MATLABを使って音波の伝搬をシミュレーション！FDTD法で挑戦！#MATLAB #音波 #シミュレーション」を追加