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

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

1. 音波とは何か

音波（sound wave）とは、空気やその他の物質、すなわち媒質（medium）中に存在する波動である。媒質粒子が振動することで、その振動が媒質中を伝搬していく現象である。
したがって、音波を解析する際には、媒質粒子の振動そのものと、その振動が媒質中を伝搬する様子の両方を考慮する必要がある。

以前に作成した音響 FDTD 1 次元バージョンを拡張し、MATLAB によって音響 FDTD 2 次元バージョンのプログラムを作成した。

king-masashi.hatenablog.com

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

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

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

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

2-1. 伝搬解析に用いたパラメータ

音波の伝搬解析に用いたパラメータ(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)

2-2. MATLABプログラム

以下に、音響 FDTD 2 次元バージョンの 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

3. 参考文献

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

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

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

[rakuten:book:10685814:detail]

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 法でシミュレーション」へ変更
• 2025年11月3日　Copilot を用いてリファクタリング
• 2024年1月28日　検索エンジン向けタイトル「音波の伝搬をMATLABでシミュレーションする方法」，SNS 向けタイトル「MATLABを使って音波の伝搬をシミュレーション！FDTD法で挑戦！#MATLAB #音波 #シミュレーション」を追加
• 2025年1月11日　記事の概要を追加
• 2026年1月11日　見出しに番号を付与