音響やオーディオ技術に関する専門用語、音響やオーディオの研究開発において必須な物理、数学等の学問で用いられる用語、音響やオーディオ技術のトレンドとして注目している技術に関する用語などを紹介します。
キーワード：音響解析、有限差分時間領域法、FDTD
近年、音響の数値解析に用いられる手法は主に3種類があります。有限要素法（Finite Element Method：FEM）、境界要素法（Boundary Element Method：BEM）と有限差分時間領域法（Finite Difference Time Domain Method：FDTD）です。前二者は定常状態（注1）の予測によく用いられるのに対し、FDTD法は時間軸に沿って、初期状態から指定の時間あるいは定常状態になるまですべての時刻を計算する手法という特徴があります。
まずスピーカー開発にFEM法とBEM法がどう応用されているかを簡単に紹介します。FEM法は固体と流体の連成解析を効率的に行えるため、振動板の形状と材質の組み合わせによって発生する音がどのように変化するかを検証することができます（図1(a)参照）。BEM法は境界上の音響特性を算出することができるため、スピーカーを中心とした半径1mや2mの音響特性を計算することによりスピーカーの指向周波数特性を検証することができます（図1(b)参照）。
図1 スピーカー開発に利用するFEM法とBEM法
FDTD法はスピーカー開発に利用する例はまだ少なく、音波伝搬を重視する建築音響分野で広く用いられています[1、2]。
FDTD法は本来電磁波の解析に用いられる手法です。マクスウェル方程式に基づき、磁界と電界相互の関係式を用いて、観測空間における電磁波分布の時間的変化を算出することができます。同様なアプローチにより、同じ波動現象である音波においても、FDTD法による音響解析（注２）が可能となります。
注１：定常状態とは、時間によって変化が生じない状態です。
注２：音響解析とは、空間における任意点での音波形を算出し、さらに各音響指標（周波数特性、残響時間、遮音性能や明瞭度など）を求めます。
FDTD法による音場解析の面白さは、音波が障害物に到来した際の反射、回折、干渉などの現象を求めることができる点です。この時間変化を動画として可視化すれば、音波伝搬の様子を明示することができます。
ここで、音波伝搬解析の例を示します。図2と3は1平方メートルの空間領域を仮定し、その周りに対しては、自由音場のような吸収境界条件を設定しています。図2では、空間中央にインパルス音源を与えると、音波が音源を中心とした円として伝搬することを確認できます。さらに図3では、右下に正方形の障害物を配置することにより、障害物による反射と回折の現象を観測できます。
図2 FDTD法による音響解析動画
図3 FDTD法による音響解析動画：右下に障害物がある場合
次に、FDTD法による音響解析において必要となる基本的な要素を紹介します。それは定式化、安定条件、境界条件および音源の設定です。これらの要素については、より複雑な手法がありますが、いずれも微積分や計算機の理論の応用ですので、本記事では割愛します。
（１）音波伝搬の数理モデルとFDTD法の定式化 まず、音波伝搬のメカニズムから説明します。声帯やスピーカーの振動板などが音の発生源として振動し始めたとき、周辺の空気粒子はこの振動を受け、動き始めます。この動きに伴い、空気密度が変わり、圧力（音圧）変化も生じます（状態１）。そして圧力変化によって近傍の空気粒子もまた動き始めます（状態２）。この空気粒子の繰り返しの周期運動によって音波が形成され、音波が音源から周りに伝搬していきます。
以上の各状態をシンプルな一次元の数理モデルで表現します。
状態１：図4.aのように、観測したい空気格子に空気粒子が密度$ρ_0$と速度$u$で流入し、$u+du$の速度で流出する場合、この一瞬における格子内の質量変化率は以下のように表現できます。 \begin{gather} \frac{d \rho^\prime}{dt}dx = \rho_0 u - \rho_0 (u + du) = - \rho_0 du \nonumber \\[10pt] \frac{d \rho^\prime}{dt} = - \rho_0 \frac{du}{dx} \end{gather} ここで$ρ'$は格子内の空気密度です。また空気の状態方程式から$P=ρ'c^2$の関係があるので、式（１）は（２）となります。 \begin{gather} \frac{dP}{dt} = - \rho_0 c^2 \frac{du}{dx} \end{gather} ここで$c$は空気中の音速です。
状態２：図4.bのような空気格子の両側に圧力差$dP$が生じた場合、運動の第2法則により格子の中の空気粒子は速度$u$で動きます。その圧力と速度の関係は以下になります。 \begin{gather} \frac{dP}{dx} = -\rho_0 \frac{du}{dt} \end{gather}
図4 連続方程式と運動方程式の一次元モデル
式（２）と（３）は音波伝搬の連続方程式と運動方程式といいます。FDTD法はこの２つの微分方程式を差分の式に変換し、四則演算のみで音波伝搬を計算することができます。 式（２）の差分式は以下になります。 \begin{gather} \frac{P_i^{n+1} - P_i^n}{\Delta t} = - \rho_0 c^2 \frac{u_{i+1/2}^{n+1/2}-u_{i-1/2}^{n+1/2}}{\Delta x} \nonumber \\[12pt] P_i^{n+1} = P_i^n - \Delta t \left( \rho_0 c^2 \frac{u_{i+1/2}^{n+1/2}-u_{i-1/2}^{n+1/2}}{\Delta x} \right) \end{gather} 式（３）の差分式は以下になります。 \begin{gather} \frac{P_i^n - P_{i-1}^n}{\Delta x} = -\rho_0 \frac{u_{i-1/2}^{n+1/2}-u_{i-1/2}^{n-1/2}}{\Delta t} \nonumber \\[12pt] u_{i-1/2}^{n+1/2} = u_{i-1/2}^{n-1/2} - \frac{\Delta t}{\rho_0}\left( \frac{P_i^n-P_{i-1}^n}{\Delta x} \right) \end{gather} 上記の式に用いている記号は下記の通りです。
$P$：音圧
$u$：$x$方向の粒子速度
$∆t$：時間の刻み幅
$∆x$：空間の刻み幅
$i$：空間領域における格子の位置
$n$：時間ステップ
上記数式における1/2とは、空間グリッドと時間ステップにおける音圧と空気粒子速度が半分ずれることを意味しています。
この数理モデルを分かりやすく説明すれば、前の時刻の音圧と粒子速度を利用し、次の時刻の音圧と粒子速度を予測するという手法です。
以上の式を２次元に拡張すると音圧と粒子速度は図5のような分布になります。次の時刻の音圧を現時刻の音圧と上下左右の粒子速度から推測し、次の時刻の粒子速度を現時刻の粒子速度と同方向にある音圧（$x$方向なら左右、$y$方向なら上下の音圧）から推測します。
図5 音圧と空気粒子速度の空間配置
数式は以下となります。 \begin{gather} P_{i,j}^{n+1}=P_{i,j}^n - \Delta t \rho_0 c^2\left( \frac{u_{i+1/2,j}^{n+1/2}-u_{i-1/2,j}^{n+1/2}}{\Delta x} + \frac{v_{i,j+1/2}^{n+1/2}-v_{i,j-1/2}^{n+1/2}}{\Delta y} \right)\end{gather} \begin{gather} u_{i-1/2,j}^{n+1/2} = u_{i-1/2,j}^{n-1/2} - \frac{\Delta t}{\rho_0}\left( \frac{P_{i,j}^n-P_{i-1,j}^n}{\Delta x} \right)\end{gather} \begin{gather} v_{i,j-1/2}^{n+1/2} = v_{i,j-1/2}^{n-1/2} - \frac{\Delta t}{\rho_0}\left( \frac{P_{i,j}^n-P_{i,j-1}^n}{\Delta y} \right)\end{gather} 上記の式に用いている記号は下記の通りです。
$P$：音圧
$u$、$v$：$x$と$y$方向の粒子速度
$∆t$：時間の刻み幅
$∆x$、$∆y$：空間の刻み幅
$i$、$j$：空間領域における格子の位置
$n$：時間ステップ
以上説明してきたのは、FDTD法の最もシンプルな差分化手法です。ただし、微分方程式を差分化することによって、誤差が必ず生じます。より高精度（誤差の低減）の計算手法が多く存在していますが、一例として、空間的にもっと広い範囲の音圧と空気粒子速度を参照し、次の時刻の結果を推測するFDTD（２，４）法という手法があります[3]。
（２）安定条件
FDTD法による音響解析を定式化した後は空間幅と時間幅の設定を行います。一般には、音響解析の結果を求めるために、長時間をかけて定常状態まで計算します。FDTD法による計算が発散しないように適切な空間幅と時間幅の設定値があり、これを安定条件またはCFL条件と呼びます。一次元の場合、安定条件は下記の式で示します。 \begin{gather} c\Delta t \leq \Delta x \end{gather} 二次元の場合は下記の式です。
\begin{gather} c\Delta t \leq \frac{1}{\sqrt{ \Bigl( \frac{1}{\Delta x} \Bigr) ^2 + \Bigl( \frac{1}{\Delta y} \Bigr) ^2 }} \end{gather} 簡潔に説明すると、時間ステップ一回の計算による情報伝搬の距離は空間格子の刻み幅より短くしなければなりません。二、三次元の場合はさらに短くする必要があります。FDTD法は現在の空間における情報を利用し、少し先の未来を推測方法です。時間幅を大きく設定してしまうと、一気に遠く未来を推測することになってしまうため、現実との乖離が生じますます。一方、空間幅の設定は比較的自由ですが、解析したい周波数が高いほど空間幅を小さく設定しなければなりません。
（３）境界条件
計算を実行する上で、無限大の空間を前提とすることは現実的に不可能ですので、境界面を設定することが必要となります。
FDTD法による音響解析によく用いられる境界条件には二種類あります。
（３－１）剛壁条件
境界面において音波が完全に反射されるという条件です。一般的な部屋、残響室、拡散材の設計によく使われています。境界面上の空気の粒子速度を0に固定することにより設定できます。
（３－２）完全吸収条件
境界面において音波が完全に吸収されるという条件です。仮想の無響室や自由音場における解析を行いたい場合はこの条件が利用されます。剛壁条件に比べて完全吸収条件はより重要ですが、その設定も複雑です。完全吸収条件を使用することで、計算領域を小さくしても計算結果が境界面による反射に影響されないほか、メモリと計算時間も節約できます。ただし、完全吸収条件を正しく設定していない場合、不要な反射が生じ、計算結果に誤差が生じます。 ここでよく使われている３種類の完全吸収条件を紹介します：
（３－２－１）Murの一次吸収境界条件[4] \begin{gather} P_{0,j}^{n+1} = P_{1,j}^{n} + \left( \frac{c-\Delta x}{c+\Delta x} \right)\left( P_{1,j}^{n+1} - P_{0,j}^{n} \right) \end{gather} \begin{gather} P_{i,0}^{n+1} = P_{i,1}^{n} + \left( \frac{c-\Delta y}{c+\Delta y} \right)\left( P_{i,1}^{n+1} - P_{i,0}^{n} \right) \end{gather} 式11と12は、それぞれ左$（x=0）$と下$（y=0）$の境界の例です。これはよく使用されている吸収境界条件の中で、最も簡単な式です。式11で、左の境界面の音圧を計算するには、その横隣の音圧のみ利用します。同じく、式12で、下の境界面の音圧を計算するには、一個上の音圧のみ利用します。このため、境界面に垂直入射する音波は完全に吸音されますが、境界面に対し斜め入射する音波はその垂直入射成分のみ消えて、水平成分は反射音として残ってしまいます。
（３－２－２）Higdonの二次吸収境界条件[5]
水平成分による反射音をより抑えたい場合には、下記の式で示した二次吸収境界条件を使用します。
\begin{equation}\begin{split} P_{0,j}^{n+1}=2\left(\frac{c\Delta t-\Delta x}{c\Delta t+\Delta x}\right)\left(P_{1,j}^{n+1}-P_{0,j}^{n}\right)+\left(\frac{c\Delta t-\Delta x}{c\Delta t+\Delta x}\right)^2\left(P_{2,j}^{n+1}-2P_{1,j}^{n}+P_{0,j}^{n-1}\right)-\\ 2τ\left(\frac{c\Delta t-\Delta x}{c\Delta t+\Delta x}\right)\left(P_{0,j}^{n}-P_{1,j}^{n-1}\right)+2τP_{1,j}^{n-1}-τ^2P_{0,j}^{n-1} \end{split}\end{equation} \begin{equation}\begin{split} P_{i,0}^{n+1}=2\left(\frac{c\Delta t-\Delta y}{c\Delta t+\Delta y}\right)\left(P_{i,1}^{n+1}-P_{i,0}^{n}\right)+\left(\frac{c\Delta t-\Delta y}{c\Delta t+\Delta y}\right)^2\left(P_{i,2}^{n+1}-2P_{i,1}^{n}+P_{i,0}^{n-1}\right)-\\ 2τ\left(\frac{c\Delta t-\Delta y}{c\Delta t+\Delta y}\right)\left(P_{i,0}^{n}-P_{i,1}^{n-1}\right)+2τP_{i,1}^{n-1}-τ^2P_{i,0}^{n-1} \end{split}\end{equation} ここで$τ=0.995$です。
上記二つの境界条件は数理的に反射音を消す方法です。
（３－２－３）完全適合層（Perfect Matched Layer：PML） [6]
反射音をさらに抑えるための条件です。
この境界条件は無響室のくさび状吸音材に近似した概念です。入射音と境界面の音響インピーダンスを完全にマッチングさせることによって反射音を消します。ただし、無響室の吸音層と同じく、PML条件による吸音効果は境界層の数によって左右されます。境界層を厚く設定すればよい吸音効果が得られますが、その分メモリと計算時間を増やす必要があります。実例として、20層の吸音層を設ければ十分な吸音効果を得られるという報告があります[3]。
（４）音源
FDTD法による音源の設定について特に計算不能となる制限はありません。ただし、現実に存在しない音源を設定すれば、意味がない結果が得られてしまいます。サイン波とガウシアンパルス音源が音源として用いられています。
サイン波は単一周波数を観測する時に用います。当該周波数の計算過程を観測することで、安定条件と境条件の設定が適切かどうかを確認できます。
ガウシアンパルスは広い周波数帯域を含む信号です。空間の音圧分布を初期条件として与えるだけでよいので、時間領域における継続した計算が不要となります。音源の設定方法は以下になります。
\begin{gather} P(x,y)=e^{-\left(\frac{r}{R}\right)^2} \end{gather} \begin{gather} r=\sqrt{\left(x-x_s\right)^2+\left(y-y_s\right)^2}{、}{R=N\Delta h}\nonumber\end{gather} 上記の式に用いている記号は下記の通りです。
$P$：音圧
$e$：ネイピア数
$x$、$y$：空間領域における位置
$x_s$、$y_s$：空間領域における音源の中心
$N$：任意定数、パルスに含む周波数帯域を制御する
$∆h$：空間幅（空間格子は正方形の場合）
ほかにも多くの音源の設定手法がありますが、たとえ音源の方程式が連続だとしても、空間と時間幅の設定によって不連続とみなされ（幅が狭いパルスなど）、解析に誤差が生じることに注意が必要です。
以上、FDTD法における４つの基本要素を紹介しました。多少のプログラミングスキルがあればシンプルな音響解析が実現できます。前述のように空間領域での計算結果を動画にできるほか、時間領域での計算結果をwavファイルとして出力すれば、可聴化することもできます。実に便利な数値解析手法です。
それぞれの要素について、現在も研究が継続されています。その目的は、高精度化と誤差低減です。なお、計算機の進歩にともない、高速化も研究テーマの一つになっています。
このFDTD法による音響解析は、室内空間における拡散音場の検討と、開放空間における音波伝搬の予測が主要な応用分野となっています。FDTD法の数理モデルを拡張すれば、多孔質材、固体における音波伝搬の解析も可能になり、空気との連成解析もできます。つまり、室内に吸音材がある場合の音響解析と、壁の材料による隣室間の遮音効果の予測も可能になります。
FDTD法による音響解析のデメリットは計算時間と使用メモリが膨大になることです、そのため、まだ広く使われている音響解析の手法とはなっていません。しかしながら、近年の研究により、多くのケースで高速計算が実現できるようになってきました。したがって、FDTD法による音響解析は、今後もっと多くの場面に利用できるようになるのではないかと期待しています。
参考文献
[1] 牛山歩，長友宏，坂本慎一，橘秀樹, "3次元FDTD解析による駒場コンベンションホールの音場予測," 日本音響学会講演論文集，pp.983-984, 2005.
[2] 横田 考俊, 坂本 慎一, 橘 秀樹, 石井 聖光, “FDTD法による鳴き竜現象の数値解析と可聴化,” 日本建築学会環境系論文集, vol. 73, no. 629, pp. 849-856, 2008.
https://doi.org/10.3130/aije.73.849
[3] Tetsuya Sakuma, Shinichi Sakamoto, Toru Otsuru, Computational Simulation in Architectural and Environmental Acoustics, Springer, 2014.
https://doi.org/10.1007/978-4-431-54454-8
[4] G. Mur, "Absorbing Boundary Conditions for the Finite-Difference Approximation of the Time-Domain Electromagnetic-Field Equations," in IEEE Transactions on Electromagnetic Compatibility, vol. EMC-23, no. 4, pp. 377-382, 1981.
https://doi.org/10.1109/TEMC.1981.303970
[5] Robert L. Higdon, "Absorbing boundary conditions for difference approximations to the multi-dimensional wave equation," Mathematics of computation, vol. 47, no. 176, pp. 437-459, 1986.
https://doi.org/10.2307/2008166
[6] J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 114, no. 2, pp. 185–200, 1994.
https://doi.org/10.1006/jcph.1994.1159
Seki
まずスピーカー開発にFEM法とBEM法がどう応用されているかを簡単に紹介します。FEM法は固体と流体の連成解析を効率的に行えるため、振動板の形状と材質の組み合わせによって発生する音がどのように変化するかを検証することができます（図1(a)参照）。BEM法は境界上の音響特性を算出することができるため、スピーカーを中心とした半径1mや2mの音響特性を計算することによりスピーカーの指向周波数特性を検証することができます（図1(b)参照）。
FDTD法はスピーカー開発に利用する例はまだ少なく、音波伝搬を重視する建築音響分野で広く用いられています[1、2]。
FDTD法は本来電磁波の解析に用いられる手法です。マクスウェル方程式に基づき、磁界と電界相互の関係式を用いて、観測空間における電磁波分布の時間的変化を算出することができます。同様なアプローチにより、同じ波動現象である音波においても、FDTD法による音響解析（注２）が可能となります。
注１：定常状態とは、時間によって変化が生じない状態です。
注２：音響解析とは、空間における任意点での音波形を算出し、さらに各音響指標（周波数特性、残響時間、遮音性能や明瞭度など）を求めます。
FDTD法による音場解析の面白さは、音波が障害物に到来した際の反射、回折、干渉などの現象を求めることができる点です。この時間変化を動画として可視化すれば、音波伝搬の様子を明示することができます。
ここで、音波伝搬解析の例を示します。図2と3は1平方メートルの空間領域を仮定し、その周りに対しては、自由音場のような吸収境界条件を設定しています。図2では、空間中央にインパルス音源を与えると、音波が音源を中心とした円として伝搬することを確認できます。さらに図3では、右下に正方形の障害物を配置することにより、障害物による反射と回折の現象を観測できます。
次に、FDTD法による音響解析において必要となる基本的な要素を紹介します。それは定式化、安定条件、境界条件および音源の設定です。これらの要素については、より複雑な手法がありますが、いずれも微積分や計算機の理論の応用ですので、本記事では割愛します。
（１）音波伝搬の数理モデルとFDTD法の定式化 まず、音波伝搬のメカニズムから説明します。声帯やスピーカーの振動板などが音の発生源として振動し始めたとき、周辺の空気粒子はこの振動を受け、動き始めます。この動きに伴い、空気密度が変わり、圧力（音圧）変化も生じます（状態１）。そして圧力変化によって近傍の空気粒子もまた動き始めます（状態２）。この空気粒子の繰り返しの周期運動によって音波が形成され、音波が音源から周りに伝搬していきます。
以上の各状態をシンプルな一次元の数理モデルで表現します。
状態１：図4.aのように、観測したい空気格子に空気粒子が密度$ρ_0$と速度$u$で流入し、$u+du$の速度で流出する場合、この一瞬における格子内の質量変化率は以下のように表現できます。 \begin{gather} \frac{d \rho^\prime}{dt}dx = \rho_0 u - \rho_0 (u + du) = - \rho_0 du \nonumber \\[10pt] \frac{d \rho^\prime}{dt} = - \rho_0 \frac{du}{dx} \end{gather} ここで$ρ'$は格子内の空気密度です。また空気の状態方程式から$P=ρ'c^2$の関係があるので、式（１）は（２）となります。 \begin{gather} \frac{dP}{dt} = - \rho_0 c^2 \frac{du}{dx} \end{gather} ここで$c$は空気中の音速です。
状態２：図4.bのような空気格子の両側に圧力差$dP$が生じた場合、運動の第2法則により格子の中の空気粒子は速度$u$で動きます。その圧力と速度の関係は以下になります。 \begin{gather} \frac{dP}{dx} = -\rho_0 \frac{du}{dt} \end{gather}
式（２）と（３）は音波伝搬の連続方程式と運動方程式といいます。FDTD法はこの２つの微分方程式を差分の式に変換し、四則演算のみで音波伝搬を計算することができます。 式（２）の差分式は以下になります。 \begin{gather} \frac{P_i^{n+1} - P_i^n}{\Delta t} = - \rho_0 c^2 \frac{u_{i+1/2}^{n+1/2}-u_{i-1/2}^{n+1/2}}{\Delta x} \nonumber \\[12pt] P_i^{n+1} = P_i^n - \Delta t \left( \rho_0 c^2 \frac{u_{i+1/2}^{n+1/2}-u_{i-1/2}^{n+1/2}}{\Delta x} \right) \end{gather} 式（３）の差分式は以下になります。 \begin{gather} \frac{P_i^n - P_{i-1}^n}{\Delta x} = -\rho_0 \frac{u_{i-1/2}^{n+1/2}-u_{i-1/2}^{n-1/2}}{\Delta t} \nonumber \\[12pt] u_{i-1/2}^{n+1/2} = u_{i-1/2}^{n-1/2} - \frac{\Delta t}{\rho_0}\left( \frac{P_i^n-P_{i-1}^n}{\Delta x} \right) \end{gather} 上記の式に用いている記号は下記の通りです。
$P$：音圧
$u$：$x$方向の粒子速度
$∆t$：時間の刻み幅
$∆x$：空間の刻み幅
$i$：空間領域における格子の位置
$n$：時間ステップ
上記数式における1/2とは、空間グリッドと時間ステップにおける音圧と空気粒子速度が半分ずれることを意味しています。
この数理モデルを分かりやすく説明すれば、前の時刻の音圧と粒子速度を利用し、次の時刻の音圧と粒子速度を予測するという手法です。
以上の式を２次元に拡張すると音圧と粒子速度は図5のような分布になります。次の時刻の音圧を現時刻の音圧と上下左右の粒子速度から推測し、次の時刻の粒子速度を現時刻の粒子速度と同方向にある音圧（$x$方向なら左右、$y$方向なら上下の音圧）から推測します。
数式は以下となります。 \begin{gather} P_{i,j}^{n+1}=P_{i,j}^n - \Delta t \rho_0 c^2\left( \frac{u_{i+1/2,j}^{n+1/2}-u_{i-1/2,j}^{n+1/2}}{\Delta x} + \frac{v_{i,j+1/2}^{n+1/2}-v_{i,j-1/2}^{n+1/2}}{\Delta y} \right)\end{gather} \begin{gather} u_{i-1/2,j}^{n+1/2} = u_{i-1/2,j}^{n-1/2} - \frac{\Delta t}{\rho_0}\left( \frac{P_{i,j}^n-P_{i-1,j}^n}{\Delta x} \right)\end{gather} \begin{gather} v_{i,j-1/2}^{n+1/2} = v_{i,j-1/2}^{n-1/2} - \frac{\Delta t}{\rho_0}\left( \frac{P_{i,j}^n-P_{i,j-1}^n}{\Delta y} \right)\end{gather} 上記の式に用いている記号は下記の通りです。
$P$：音圧
$u$、$v$：$x$と$y$方向の粒子速度
$∆t$：時間の刻み幅
$∆x$、$∆y$：空間の刻み幅
$i$、$j$：空間領域における格子の位置
$n$：時間ステップ
以上説明してきたのは、FDTD法の最もシンプルな差分化手法です。ただし、微分方程式を差分化することによって、誤差が必ず生じます。より高精度（誤差の低減）の計算手法が多く存在していますが、一例として、空間的にもっと広い範囲の音圧と空気粒子速度を参照し、次の時刻の結果を推測するFDTD（２，４）法という手法があります[3]。
（２）安定条件
FDTD法による音響解析を定式化した後は空間幅と時間幅の設定を行います。一般には、音響解析の結果を求めるために、長時間をかけて定常状態まで計算します。FDTD法による計算が発散しないように適切な空間幅と時間幅の設定値があり、これを安定条件またはCFL条件と呼びます。一次元の場合、安定条件は下記の式で示します。 \begin{gather} c\Delta t \leq \Delta x \end{gather} 二次元の場合は下記の式です。
\begin{gather} c\Delta t \leq \frac{1}{\sqrt{ \Bigl( \frac{1}{\Delta x} \Bigr) ^2 + \Bigl( \frac{1}{\Delta y} \Bigr) ^2 }} \end{gather} 簡潔に説明すると、時間ステップ一回の計算による情報伝搬の距離は空間格子の刻み幅より短くしなければなりません。二、三次元の場合はさらに短くする必要があります。FDTD法は現在の空間における情報を利用し、少し先の未来を推測方法です。時間幅を大きく設定してしまうと、一気に遠く未来を推測することになってしまうため、現実との乖離が生じますます。一方、空間幅の設定は比較的自由ですが、解析したい周波数が高いほど空間幅を小さく設定しなければなりません。
（３）境界条件
計算を実行する上で、無限大の空間を前提とすることは現実的に不可能ですので、境界面を設定することが必要となります。
FDTD法による音響解析によく用いられる境界条件には二種類あります。
（３－１）剛壁条件
境界面において音波が完全に反射されるという条件です。一般的な部屋、残響室、拡散材の設計によく使われています。境界面上の空気の粒子速度を0に固定することにより設定できます。
（３－２）完全吸収条件
境界面において音波が完全に吸収されるという条件です。仮想の無響室や自由音場における解析を行いたい場合はこの条件が利用されます。剛壁条件に比べて完全吸収条件はより重要ですが、その設定も複雑です。完全吸収条件を使用することで、計算領域を小さくしても計算結果が境界面による反射に影響されないほか、メモリと計算時間も節約できます。ただし、完全吸収条件を正しく設定していない場合、不要な反射が生じ、計算結果に誤差が生じます。 ここでよく使われている３種類の完全吸収条件を紹介します：
（３－２－１）Murの一次吸収境界条件[4] \begin{gather} P_{0,j}^{n+1} = P_{1,j}^{n} + \left( \frac{c-\Delta x}{c+\Delta x} \right)\left( P_{1,j}^{n+1} - P_{0,j}^{n} \right) \end{gather} \begin{gather} P_{i,0}^{n+1} = P_{i,1}^{n} + \left( \frac{c-\Delta y}{c+\Delta y} \right)\left( P_{i,1}^{n+1} - P_{i,0}^{n} \right) \end{gather} 式11と12は、それぞれ左$（x=0）$と下$（y=0）$の境界の例です。これはよく使用されている吸収境界条件の中で、最も簡単な式です。式11で、左の境界面の音圧を計算するには、その横隣の音圧のみ利用します。同じく、式12で、下の境界面の音圧を計算するには、一個上の音圧のみ利用します。このため、境界面に垂直入射する音波は完全に吸音されますが、境界面に対し斜め入射する音波はその垂直入射成分のみ消えて、水平成分は反射音として残ってしまいます。
（３－２－２）Higdonの二次吸収境界条件[5]
水平成分による反射音をより抑えたい場合には、下記の式で示した二次吸収境界条件を使用します。
\begin{equation}\begin{split} P_{0,j}^{n+1}=2\left(\frac{c\Delta t-\Delta x}{c\Delta t+\Delta x}\right)\left(P_{1,j}^{n+1}-P_{0,j}^{n}\right)+\left(\frac{c\Delta t-\Delta x}{c\Delta t+\Delta x}\right)^2\left(P_{2,j}^{n+1}-2P_{1,j}^{n}+P_{0,j}^{n-1}\right)-\\ 2τ\left(\frac{c\Delta t-\Delta x}{c\Delta t+\Delta x}\right)\left(P_{0,j}^{n}-P_{1,j}^{n-1}\right)+2τP_{1,j}^{n-1}-τ^2P_{0,j}^{n-1} \end{split}\end{equation} \begin{equation}\begin{split} P_{i,0}^{n+1}=2\left(\frac{c\Delta t-\Delta y}{c\Delta t+\Delta y}\right)\left(P_{i,1}^{n+1}-P_{i,0}^{n}\right)+\left(\frac{c\Delta t-\Delta y}{c\Delta t+\Delta y}\right)^2\left(P_{i,2}^{n+1}-2P_{i,1}^{n}+P_{i,0}^{n-1}\right)-\\ 2τ\left(\frac{c\Delta t-\Delta y}{c\Delta t+\Delta y}\right)\left(P_{i,0}^{n}-P_{i,1}^{n-1}\right)+2τP_{i,1}^{n-1}-τ^2P_{i,0}^{n-1} \end{split}\end{equation} ここで$τ=0.995$です。
上記二つの境界条件は数理的に反射音を消す方法です。
（３－２－３）完全適合層（Perfect Matched Layer：PML） [6]
反射音をさらに抑えるための条件です。
この境界条件は無響室のくさび状吸音材に近似した概念です。入射音と境界面の音響インピーダンスを完全にマッチングさせることによって反射音を消します。ただし、無響室の吸音層と同じく、PML条件による吸音効果は境界層の数によって左右されます。境界層を厚く設定すればよい吸音効果が得られますが、その分メモリと計算時間を増やす必要があります。実例として、20層の吸音層を設ければ十分な吸音効果を得られるという報告があります[3]。
（４）音源
FDTD法による音源の設定について特に計算不能となる制限はありません。ただし、現実に存在しない音源を設定すれば、意味がない結果が得られてしまいます。サイン波とガウシアンパルス音源が音源として用いられています。
サイン波は単一周波数を観測する時に用います。当該周波数の計算過程を観測することで、安定条件と境条件の設定が適切かどうかを確認できます。
ガウシアンパルスは広い周波数帯域を含む信号です。空間の音圧分布を初期条件として与えるだけでよいので、時間領域における継続した計算が不要となります。音源の設定方法は以下になります。
\begin{gather} P(x,y)=e^{-\left(\frac{r}{R}\right)^2} \end{gather} \begin{gather} r=\sqrt{\left(x-x_s\right)^2+\left(y-y_s\right)^2}{、}{R=N\Delta h}\nonumber\end{gather} 上記の式に用いている記号は下記の通りです。
$P$：音圧
$e$：ネイピア数
$x$、$y$：空間領域における位置
$x_s$、$y_s$：空間領域における音源の中心
$N$：任意定数、パルスに含む周波数帯域を制御する
$∆h$：空間幅（空間格子は正方形の場合）
ほかにも多くの音源の設定手法がありますが、たとえ音源の方程式が連続だとしても、空間と時間幅の設定によって不連続とみなされ（幅が狭いパルスなど）、解析に誤差が生じることに注意が必要です。
以上、FDTD法における４つの基本要素を紹介しました。多少のプログラミングスキルがあればシンプルな音響解析が実現できます。前述のように空間領域での計算結果を動画にできるほか、時間領域での計算結果をwavファイルとして出力すれば、可聴化することもできます。実に便利な数値解析手法です。
それぞれの要素について、現在も研究が継続されています。その目的は、高精度化と誤差低減です。なお、計算機の進歩にともない、高速化も研究テーマの一つになっています。
このFDTD法による音響解析は、室内空間における拡散音場の検討と、開放空間における音波伝搬の予測が主要な応用分野となっています。FDTD法の数理モデルを拡張すれば、多孔質材、固体における音波伝搬の解析も可能になり、空気との連成解析もできます。つまり、室内に吸音材がある場合の音響解析と、壁の材料による隣室間の遮音効果の予測も可能になります。
FDTD法による音響解析のデメリットは計算時間と使用メモリが膨大になることです、そのため、まだ広く使われている音響解析の手法とはなっていません。しかしながら、近年の研究により、多くのケースで高速計算が実現できるようになってきました。したがって、FDTD法による音響解析は、今後もっと多くの場面に利用できるようになるのではないかと期待しています。
参考文献
[1] 牛山歩，長友宏，坂本慎一，橘秀樹, "3次元FDTD解析による駒場コンベンションホールの音場予測," 日本音響学会講演論文集，pp.983-984, 2005.
[2] 横田 考俊, 坂本 慎一, 橘 秀樹, 石井 聖光, “FDTD法による鳴き竜現象の数値解析と可聴化,” 日本建築学会環境系論文集, vol. 73, no. 629, pp. 849-856, 2008.
https://doi.org/10.3130/aije.73.849
[3] Tetsuya Sakuma, Shinichi Sakamoto, Toru Otsuru, Computational Simulation in Architectural and Environmental Acoustics, Springer, 2014.
https://doi.org/10.1007/978-4-431-54454-8
[4] G. Mur, "Absorbing Boundary Conditions for the Finite-Difference Approximation of the Time-Domain Electromagnetic-Field Equations," in IEEE Transactions on Electromagnetic Compatibility, vol. EMC-23, no. 4, pp. 377-382, 1981.
https://doi.org/10.1109/TEMC.1981.303970
[5] Robert L. Higdon, "Absorbing boundary conditions for difference approximations to the multi-dimensional wave equation," Mathematics of computation, vol. 47, no. 176, pp. 437-459, 1986.
https://doi.org/10.2307/2008166
[6] J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 114, no. 2, pp. 185–200, 1994.
https://doi.org/10.1006/jcph.1994.1159
Seki
