FDTD法による音響シミュレーション

  • by

pythonでFDTDシミュレーション

音響シミュレーションそのものについて、実装することに関する資料が多くない
と思ったのでノートブック形式で解説するような記事にしました。

少しばかり長くなりますがお付き合いください。

なお書籍「はじめての音響数値シミュレーションプログラミングガイド(日本建築学会編、コロナ社 2011)」を
だいぶ参考にして 本記事を書いています。詳しく知りたい方はこちらを参照ください。

Contents

  1. 音響シミュレーションについて
  2. FDTD法について
  3. 実装
  4. まとめ

1. 音響シミュレーション

音響の話に入る前に、そもそもシミュレーションとはなんなのか、
自分なりにまとめてみたいと思います。

シミュレーションとは、
モデルにもとづいて、コンピュータ上で現象を再現すること
だと考えてます。

もちろんsimulationという英語が元になっていますが、
動詞の形だとsimulate(〜のふりをする(他動詞))です。

たとえば中学、高校あたりで習う自由落下の式は、

\[y = \frac{1}{2}gt^2\] \[y: 変位[\rm{m}], g: 重力加速度[\rm{m/s^2}], t: 時間[\rm{s}]\]

で表すことができます。

これは「自由落下」という物理現象を、時間(もしくは重力加速度)と
変位の関係だけに着目して抽出したモデルといえます。

この変位yの値が時間変化によってどのように動くか、
tに代入する時間を変えてその結果を表すことは
シミュレーションにあたります。

このシミュレーション結果の正しさはモデルの正しさに影響を受けます。
自由落下の例では空気抵抗など一切考慮していないので、
実際の動きを完全にまねることはできません。

ちなみに実際の動きを100%再現するモデル化は存在しません。
また再現精度をあげるほどモデルは複雑化するので、
計算資源を過剰に浪費しかねません。
そこでどこまでの精度が必要なのか、どの程度の処理の軽さが必要なのか、
目的に合わせて議論する必要があります。

話休題。
音響シミュレーションとは読んで字のごとく、
音響の現象をコンピュータ上で再現することです。

音響の現象も様々で、音の伝わる様子や反射の様子、
もしくは吸音される仕組みや音響流などの特殊な現象もあります。

音の伝搬(伝わること)のシミュレーションには、主に2種類あり、
その中でもいくつか手法が存在します(Fig. 1)。
今回扱うFDTD法は中でも、音の伝わる様子を表すのが得意です。

Fig. 1 音響シミュレーションの分類と例

2. FDTD法

FDTDはFinite-Difference Time-Domainの頭文字をとっています。
日本語では「時間領域有限差分法」と呼びます。

内容をひとことで言うと、
音の振る舞いを表す式を差分によって表現したものを、
時間変化させていくことで音の伝搬をシミュレートする方法
です。

元は電磁気学の分野で使われていた手法ですが、
同じ波の現象である音にも応用されています。

またFDTD法の特徴のひとつとして、スタガード格子(FIg. 2)があげられます。
これは差分を利用するFDTDの特性を活かしたもので、
音圧と音速を定義する座標が互い違いになるようにしたものです。

Fig. 2 x-tのスタガード格子

元となる式(運動方程式、連続の式)

FDTDのもととなるモデルの式は、

\[ \frac{\partial p}{\partial t}
+ \kappa
(\frac{\partial u_{x}}{\partial x}
+\frac{\partial u_{y}}{\partial y}
+\frac{\partial u_{z}}{\partial z}) = 0
\] \[ p: 音圧, \kappa: 気体弾性率, u_{x(y, z)}: x(もしくはy, z)軸方向の粒子速度
\]

で表される連続の式と、

\[ \frac{\partial p}{\partial t} + \rho \frac{\partial u_{x}}{\partial t}=0 \\
\frac{\partial p}{\partial t} + \rho \frac{\partial u_{y}}{\partial t}=0 \\
\frac{\partial p}{\partial t} + \rho \frac{\partial u_{z}}{\partial t}=0
\] \[ \rho: 媒質密度
\]

で表される各軸方向の運動方程式から成り立ちます。

式の導出には尾本先生による解説 ( 「波動方程式から理解する音響学」
ございますのでこちらを参考にされるといいかと思います。

いきなり式が出てきましたが大切なのは個々の意味です。
微分項( ∂ / ∂ x)や( ∂ / ∂ t)は変化量を意味しています。
( ∂ / ∂ x) はx軸の変位による変化量、( ∂ / ∂ t) は時間変化量で、
具体的には( ∂p / ∂ t)は「音圧の時間変化量」 になります。

この微分(つまり変化量)を差分、つまり引き算の形で表します。
Fig. 3 に微分を差分で近似するイメージを示します。

Fig. 3 微分を差分で近似するイメージ

式で表すと、

\[ \frac{df(x)}{dx} \approx \frac{f(x+h)-f(x-h)}{2h}
\]

と言った感じで表せます。

この差分近似をつかってx軸の運動方程式を変形させると、

\[ u^{n+1}_{x}(i+1)
=u^{n}_{x}(i+1)
– \frac{\Delta t}{\rho \Delta x} [p^{n}(i+1)-p^{n}(i)] \]

となります。(詳しくは上で紹介した書籍で!)

数学が苦手な方だと、ウッとなるかもしれませんが、大したことはありません。
Fig. 4をもとに順に意味を追っていきましょう。

Fig. 4 運動方程式の変形式における音速と音圧の関係

連続の式についても簡単のためにx軸のみにして、
微分を差分近似して、今度は音圧に関してまとめると、

\[ p^{n+1}(i)
=p^{n}(i)-\frac{\kappa \Delta t}{\Delta x}[u^{n+1}_x(i+1)-u^{n+1}_{x}(i)] \]

と、なります。

同じように図にしてみます。(Fig. 5)

Fig. 5 連続の式の変形式における音圧と音速の関係

これらの関係をまとめたのがFig. 6です。

Fig. 6 スタガード格子上では差分が同じ位置で定義可能

ここまでの内容に音源情報と境界線(面)上での振る舞いを追加することで、
ひとつ前の時刻の値から音速、音圧を計算し進めることができます。

境界での振る舞いは音響インピーダンスによって決定します。
しかし境界の扱いは不安定になる要素であり、完全吸音面を設置したり、
高次の近似をすることなど様々な試みがされています。

3. 実装

FDTD法シミュレーションプログラム
Google ColabGithubにて公開してあります。
notebook環境+α(numpy, matplotlib, scipy)があれば動きます。

一通り実行するとFIg. 7, 8のような結果が得られます。

Fig. 7 シミュレーション結果:音圧分布
Fig. 8 シミュレーション結果:受音点での時間音圧変化

OpenAcousticsにて (GPL v2下で) 公開されているコードをもとに、
Python3でも動くように、かつnotebook形式で実行できるようにしました。
また各所にFDTD法の説明や根拠となる式なども追加しました。

プログラムの流れとしては、

  1. パッケージ等のインポート
  2. シミュレーション条件の設定
  3. 室形状の設定
  4. 音源の設定
  5. シミュレーション
  6. 可視化

のようになっています。
ちなみに室形状および音源はひとつ選んで実行していただく形になっています。

4. まとめ

シミュレーションそのものの話からFDTD法シミュレーションの中身の説明、
そしてプログラムの公開までしました。

初めてのオープンソースの改変・公開なので何か問題等あれば、
ご連絡ほどよろしくお願いします。


三笠茶屋くんぺい

旧三笠ホテルの向かい側にある喫茶店

小さなフラワーアレンジメントが素敵

茶屋と言いつつもカレーがおすすめ

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です