目次

• はじめに

• 時間調和解析

• 固有振動数解析

• ヘルムホルツ方程式

• ヘルムホルツ方程式入門

• 時間領域モデルから周波数音響モデルを導く

• モデルパラメータの設定

• 音源タイプ

• 単極音源

• 二重極音源

• 損失媒体における音の伝搬

• 複数材料媒体における音の伝搬

• 時間領域モデリングと周波数領域モデリングの比較

• 波動方程式：時間領域モデリング

• ヘルムホルツ方程式：周波数領域モデリング

• 正確性の比較

• 音響境界条件

• インピーダンス境界条件

• 定式化

• 導出

• 時間調和解析におけるインピーダンス境界条件

• 吸収境界条件

• 定式化

• 導出

• 時間調和解析における吸収境界条件

• Sound Hard境界条件（壁）

• 定式化

• 導出

• 時間調和解析におけるSound Hard境界条件

• 固有振動数解析におけるSound Hard境界条件

• 法線速度境界条件

• 定式化

• 導出

• 時間調和解析における法線速度境界条件

• Sound Soft境界条件

• 定式化

• 導出

• 時間調和解析におけるSound Soft境界条件

• 固有振動数解析におけるSound Soft境界条件

• 圧力源境界条件

• 定式化

• 導出

• 時間調和解析における圧力源境界条件

• 放射境界条件

• 定式化

• 導出

• 時間調和解析における放射境界条件

• フロケ周期境界条件

• 定式化

• 導出

• 時間調和解析におけるフロケ周期境界条件

• 完全整合層(PML)

• 時間調和解析におけるPML

• 用語集

• 参考文献

はじめに

音響学とは圧力の変化によって音波をモデル化する物理の一領域である．音響系をモデル化するためには2つのアプローチが広く使われる．一つは時間領域で音響をモデル化するというもので，もう一つは周波数領域でモデル化するというものである．このチュートリアルは周波数領域での音響モデリングに焦点を当て，モデルとしてヘルムホルツ(Helmholtz)偏微分方程式(PDE)を利用する．ここで紹介する周波数領域の音響モデリングは，音響モデリングの第一歩として読むべきチュートリアル「時間領域における音響学」で導入されている概念の上に構築する．2つのアプローチは関連しているので，この関係はこのチュートリアル全体で提示する．

周波数領域の音響学のメインモデルはヘルムホルツ方程式である．ヘルムホルツ方程式は独立方程式である．ヘルムホルツ方程式は時間独立の方程式なので，時間領域内の音響学のモデリングに使われる時間依存の波動方程式と比較して，より効率的に解くことができる．しかしヘルムホルツ方程式は調和時間依存性を持つ音響系をモデル化するときにのみ適用される．つまり，非調和音声信号は時間領域でモデル化しなければならず，ヘルムホルツ方程式を使うことによる利点は使えない．

周波数領域における2つのタイプの解析は「時間調和解析」および「固有振動解析」で紹介する．時間調和解析も固有振動解析もさまざまなタイプの境界条件を使ったヘルムホルツ偏微分方程式に基づいている．このチュートリアルでは境界条件についても取り上げる．時間調和解析の目的は，周波数範囲上で音響系の周波数応答を計算することである．その一方，固有振動解析は音響系の固有モードと固有振動を解くために適用される．時間調和モデルの実際の解析はParametricNDSolveを使って，固有振動解析はNDEigensystemを使って行われる．

音響系モデリングについての拡張例はモデルコレクションで見ることができる．

このチュートリアルで使われる記号とその単位は，用語のセクションに要約されている．

時間調和解析

時間調和解析は音響系の周波数応答に関係する音響学の一領域である．音声信号が特定の周波数の正弦関数として表せる場合，それは時間調和といわれる．時間調和解析では，音響系はある範囲の周波数上でいくつかの調和音響信号にさらされ，関心のある周波数でのデバイスの性能が解析される．この種の解析は，高周波数で音を減衰させるよう設計されたローパスフィルタ等の周波数依存の音響系を構築する場合に重要である．

調和刺激に応えて，結果の音圧場 も同様に時間調和であることを示すことができる[1]．任意の空間位置 が角周波数 の正弦波時間依存性を持つ場合，音圧場は時間依存であるといわれる．

調和音圧場 の一般式は以下のように書くことができる：

ここでは与えられた位置における振幅を表し， は における初期の位相シフトである．2つの音声信号に位相の差がない場合，は「同相」といわれる．位相の差が の場合，2つの音声信号は「逆相」といわれる．

解析的な簡便性のため，通常，時間調和関係(2)は複素指数関係(CER)として知られる複素形式で表される：

慣習的に，CER式は次のように簡単に表されることが多い：

ここで，複素数式の実部は実関数 を表しているものと解釈される．

CER式(3)は複素平面の回転ベクトル として解釈することができる．次の図は，この動作を表している．

回転ベクトル は複素振幅関数として知られる．振幅関数 は角周波数 で反時計回りに回転する．指定された時間 において，実軸への の投影は一瞬の音圧を表し，ベクトル長は局所的振幅に当たる．

振幅関数とその複素共役をそれぞれ および で表すと，局所的振幅は以下で計算することができる：

時間調和関係(4)が波動方程式に挿入されると，方程式は時間独立のヘルムホルツ(Helmholtz)方程式に簡約される．波動方程式からのヘルムホルツ方程式の導関数は時間領域モデルから周波数音響モデルを導くというセクションで説明する．ここでは，ヘルムホルツPDEモデル(5)の角周波数 を使うと，未知の音声場が周波数領域について解けることを理解することが大切である：

項 および はそれぞれ単極音源と二重極音源である．

しかし，計算された解 は，時間調和関係(6)を使って簡単に時間領域に変換することができる．

固有振動数解析

単極音源 と二重極音源 がヘルムホルツPDEモデル(7)から削除されると，この方程式はsource-freeのヘルムホルツ方程式に簡約される．

方程式(8)は となるような固有値問題として扱うことができ，NDEigensystemで解くことができる．ここで微分演算子 は(9)の左辺に当たり，は固有系の固有値を表す．

source-freeヘルムホルツ方程式を満足する固有値 の集合は，以下のいずれかによって，対応する固有振動数 を与える：

あるいは

それぞれの固有値 と対になる振幅関数 は，固有モードと呼ばれる．

固有振動数 は自然振動数ともいわれ，音響系の共鳴を決定する．音響共鳴を示すために，開口の管と閉口の管を考えよう．両方とも空気で満たされており，長さは である．

慣習により，閉口の管とは一方の端だけが閉じている管を指す．

中圧は，開口部における周囲の基準圧に等しくなければならないので，音圧は となるようなゼロで固定される．しかし，閉口部では前進運動が不可能なため，音圧は累積されその最大値に到達する．このような境界条件のため，管は一定の振動数，つまり固有振動数でのみ音波を維持することができる．

固有振動数 のそれぞれにおいて，管の中で定常波が形成される．固有モード1に対応する最低の固有振動数は，基本周波数と呼ばれる．「時間領域における音響学」で示すように，同じ振幅の進行波と比較すると，定常波は弱い音源が必要である．つまり，音源はそれぞれの固有振動数で最も音響系を活性化するということである．

したがって，楽器，音響フィルタ，コンサートホール等の音響を利用（または妨害）する音響系を設計する場合，固有振動数解析を考慮することが大切である．

ヘルムホルツ方程式

ヘルムホルツ方程式入門

周波数領域における音響系の挙動は，関心のある周波数領域からの特定の周波数に対するヘルムホルツ方程式を繰り返し解くことによって調べられる．ヘルムホルツ方程式(10)は，特定の角周波数 における調和音圧場をモデル化するために使われる：

ヘルムホルツ方程式の従属変数は音圧 である．音圧波は密度 の媒体内で音速 で伝搬する．音圧場は，による角周波数 に関連した周波数 における調和音刺激に応答してモデル化される．

音圧は周囲の基準圧からの局所圧の逸脱 （ は位置ベクトル）として理解することができる．項 と はそれぞれ単極音源と二重極音源を表す．これらの音源については音源タイプのセクションで説明する．

インタラクティブなPDE演算子の使用については，ドキュメントのいろいろなセクションで説明されている．「偏微分方程式の数値解法」をご参照いただきたい．

チュートリアル「時間領域における音響学」で示されているように，過渡の音響モデルも同様に設定することができる．

周波数領域の音響モデルでは，時間微分項は時間調和関係(11)を使うことによって周波数依存項 に変換される．この導出方法は「時間領域モデルから周波数音響モデルを導く」で見ることができる．

このチュートリアルの例では次のモデルパラメータを使う．これらのパラメータがシミュレーション領域を定義する．

次の1Dの例は，周波数領域の音響モデルのシミュレーションを示している．最初のステップで時間調和解析が実行され，後続するステップで同じ音響モデルの固有系解析が行われる．この2つの解析タイプの間の関係が明らかになる．

次のステップでは，固有値解析を行う．

それぞれの固有振動数 において，振幅応答はその最大値に達する．

時間領域モデルから周波数音響モデルを導く

ヘルムホルツ方程式は，調和時間依存性のある波動方程式(12)から導かれる．一般的な波動方程式は以下で与えられる：

ここで項 と はそれぞれ単極音源と二重極音源である．

もし音圧場が時間調和ならば，特定の周波数 に対する時間の圧力変動は振幅関数 によって複素平面上で表すことができる：

同様に，単極音源と二重極音源はそれぞれ振幅関数 と で表すことができる：

(13)の二階時間微分を取ると次のようになる：

(14)の傾きを取る：

(15)，(16)，(17)を(18)に代入すると，波動方程式は以下のようになる：

共通項 を因数分解すると，方程式は時間非依存の非斉次ヘルムホルツ方程式(19)に簡約される：

モデルパラメータの設定

音響シミュレーションでは，支配偏微分方程式の正確な数値解を得るために，音波の波長 は十分細かいメッシュで解像しなければならない．ここでは時間調和波にMaxCellMeasureを設定する関数を作成する．

時間調和波では波長は である．最大辺長のデフォルトの解像度では. あたり12ノードに設定される．つまり，音波が伝搬する各方向において波長 に対して少なくとも12個の要素があるということである．音波を正確に解くには，通常これで十分である[20]．しかし，正確性の異なる要求に合うよう別の解像度を割り当てることができる．

多くの例で，調和音波を生成するために放射境界条件を，波反射を避けるために吸収境界条件を使う．

音源タイプ

ヘルムホルツ方程式モデル(21)は，単極音源 と二重極音源 の2つのタイプの時間調和音源を含んでいる．次のセクションでは，周波数領域でのモデリングのためにこれらの音源をどのように設定するかを示す．音源の物理的意味は「時間領域における音響学」で詳しく説明する．

単極音源

ヘルムホルツ方程式(22)の単極音源 は，音声を等方的に放射する点光源をモデル化する．音響単極音源の例として，半径が拡張と収縮を繰り返す小さい球がある[23]．

単極音源を利用するためには，単極音源の強さ と音源の位置 を指定しなければならない．単極音源項 は以下のように書ける：

ここで は音源の位置における，正規化されたディラック(Dirac)のデルタ関数である．

デルタ関数を正規化する方法は以下のようにいろいろある[24,25]：

ここで は、正規化されたデルタ関数 のサポートを制御する正規化パラメータである．通常 はメッシュ間隔 に相当するサイズを持たなければならない．

時間調和解析というセクションで述べているように，ヘルムホルツ方程式の解は複素値の振幅関数 である．これには位相と振幅の両方の情報が含まれている．音圧の振幅は絶対値あるいは複素弾性率 に対応する．

1D領域では，音圧振幅は各特定の周波数における単極音源 について一定である．音源領域 内では，使用された正規化デルタ関数 の離散的性質により，にわずかな偏差がある．より細かいメッシュを使うことで，数値誤差を減少させることができる．

与えられた単極音源の強さ では，音圧振幅は周波数が高くなるにつれて減少する．一次元の振幅に対する解析解[26]は次で与えられる：

音圧 を時間の関数として考えるとより役に立つ．そのためには，振幅関数 は調和波関係(27)を介して変換し とする必要がある．

単極音源は で等方的に音声を放射する点光源である．ここで青い線は過渡音圧 ，灰色の線は視覚的確認のために挿入された解析音圧である．

比較として，同様に2Dの単極音源を構築する．

2Dと3Dの単極の場合，音源から放射された音波が広がるにつれ，波面は幅広くなる．したがって，指定された単極音源の強さ では，音圧振幅は音源の位置への距離とともに減衰する．

二重極音源

二重極音源 は，強さが等しいが反対の位相であり距離 だけ離れた2つの単極音源からなる．この分離距離は，音源の波長と比較すると小さい．音響の二重極の例として，正弦曲線状に振動する小さく硬い球が挙げられる[28]．単極音源 とは異なり，二重極音源は等方的に音声を放射しない．

「時間領域における音響学」で，二重極音源 は2つの等しい単極音源 でモデル化できるということを示した．ここではその中間ステップは示さず，音源項 を持つ二重極音源を直接モデル化する．

二重極音源を利用するためには，二重極音源の強さ と音源の位置 を指定しなければならない．二重極音源項 は次のように書ける：

ここで は音源の位置における正規化デルタ関数である．

「時間調和解析」というセクションで示したように，ヘルムホルツ方程式の解は複素数値の振幅関数 である．これには位相と振幅両方の情報が含まれる．音圧の振幅は絶対値あるいは複素弾性率に対応する．

1D領域では，音圧振幅は二重極音源領域の外側で一定である．埋め込まれた単極音源は反対の位相なので，音源の位置における振幅 は合計するとゼロになる．1Dの音圧振幅の解析解[29]は以下で与えられる：

しかし，上の数値解のプロットは，使用された正規化デルタ関数 の離散性のため，解析解とわずかに異なる．数値誤差を小さくするためにはより微細なメッシュを使うことができる．

音圧 を時間の関数として考えるとより役に立つ．そのためには，振幅関数 を調和波関係(30)で に変換しなければならない．

二重極音源は音声を等方的に放射しない．結果の音波は正弦波であるが位相は逆である．

比較として，同様に2Dの二重極音源を構築する．

2Dと3Dの二重極では，放射された音波はすべての方向に均等に広がる訳ではない[31]．それゆえ指定された二重極音源の強さ では，音圧振幅は空間方向と音源位置への距離の両方に依存する．

損失媒体における音の伝搬

「時間領域における音響学」にあるように，指定された多孔性 ，流れ抵抗 ，有効体積弾性率 の損失媒体で音の伝搬をモデル化することができる．変更された波動方程式は次のようになる：

調和波関係 を挿入し，共通項 を因数分解すると，方程式(32)は次のようになる：

方程式(33)は，周波数領域における音の減衰をモデル化するのに使われる，修正ヘルムホルツ方程式である．

有効体積弾性率 は周波数依存であり[34]，実験式で近似することができる．

損失媒体モデルを例示するために，多孔質吸音材料で満たされた管を伝搬する音波を考える．この材料の多孔性と流れ抵抗はそれぞれ と で与えられる．

音波の反射を防ぐため，右側に吸収境界条件を加える．修正ヘルムホルツ方程式(35)に対応するためにNeumannValueを [36]に設定する．

(37)の有効体積弾性率 は周波数依存の変数なので，振幅の減衰率も周波数によって異なる．右端がわずかに波打っているのは，吸収境界への数値反射によるものであり，より細かいメッシュを使うと軽減することができる．

音圧波の振幅は，右に伝搬するにつれて減衰する．減衰率は周波数によってわずかに異なる．

時間領域モデリングと周波数領域モデリングの比較

音響系は，時間領域でも周波数領域でもモデル化することができる．「時間領域における音響学」で説明しているように，波動方程式は時間領域の音波の過渡解を求めるのに使われる．音圧場の励起が時間調和の場合は常に，周波数領域の音圧場の定常状態解を直接解くためにヘルムホルツ方程式を使うことができる．どちらのアプローチにも，以下の表にまとめたように長所と短所がある．

    

	時間領域モデリング	周波数領域モデリング
支配PDE	波動方程式	ヘルムホルツ方程式
従属変数	過渡 p(t,X)	定常 p(X)
調和励振	あり	あり
非調和励振	あり	なし
計算コスト	高	低
精度	低	高

換言すると，音響系が時間非調和依存性を持つ場合，それは時間領域でモデル化しなければならず，ヘルムホルツ方程式を使う利点は使えない．しかし他のすべての場合では，周波数領域でモデルを使うことに利点がある．

この動作を表すために，次の例では右に移動する波の音響系を考える．時間領域で1回，周波数領域で1回シミュレーションが行われる．両方のモデルで，シミュレーション領域は波長 の4倍に設定され、周波数は任意に で選ばれる．

以下の例では左側に放射境界条件を加えて調和音波を生成し，右側に吸収境界条件を置いて波の反射を避ける．

まずモデルを時間領域で解析し，周波数領域解析を実行する．互いに結果を比較する．

波動方程式：時間領域モデリング

波動方程式モデルでは，シミュレーションの終了時間 は，動的定常状態に達するための音圧場に必要な時間として定義される．

結果は右に移動する過渡の音圧波を示す．

次にヘルムホルツ方程式で同じモデルを構築する．

ヘルムホルツ方程式：周波数領域モデリング

方程式の設定は変わらない．変更はすべて変数とパラメータを介して行われる．

ヘルムホルツ方程式は音圧場の定常解を直接計算するので，解く過程で時間積分を行う必要はない．このため，計算コストが大幅に削減される．

ヘルムホルツ方程式の結果は定常状態の音圧場 である．解を時間領域に変換するために，調和波関係を使う．Esc

しかし，過渡解を周波数領域に変換するためには，フーリエ解析[38]が必要である．このプロセスは過渡音の信号を調和信号の和に分解するものである．それぞれの調和信号には特定の周波数と相対振幅があるため，過渡信号を周波数領域Escにマップすることが可能なのである．

ヘルムホルツ方程式を解くための計算的コストはずっと低いため，ヘルムホルツ偏微分方程式を異なる周波数で反復的に解くことが可能であり，これによって周波数領域モデリングに適したものにすることができる．

次に波動方程式モデルとヘルムホルツ方程式モデルの正確性を比較する．

正確性の比較

波動方程式モデルは時間依存のPDEであるがヘルムホルツ方程式モデルはそうではないため，前者は数値的時間積分による余分な誤差の影響を受ける．したがって，ヘルムホルツモデルは波動方程式よりも正確である．

音響境界条件

「時間領域における音響学」では一般的な音響境界条件の詳細と，それをどのようにして時間領域でモデル化するかについて説明した．反復を避けるため，次のセクションではこれらの境界条件を周波数領域で構築する方法だけを示す．導出や物理的説明に関心のある読者は，「時間領域における音響学」をご一読いただきたい．

音響学で最も一般的な境界条件はDirichletCondition，NeumannValue，PeriodicBoundaryConditionでモデル化でき，次の4つのタイプに分類することができる：

• ロビンタイプ

• インピーダンス境界条件

• 吸収境界条件

• 放射境界条件

• ノイマンタイプ

• Sound Hard境界条件 - 壁

• 法線速度境界条件

• 圧力源境界条件（ノイマンモデル）

• ディリクレタイプ

• Sound Soft境界条件

• 圧力源境界条件（ディリクレモデル）

• 周期タイプ

• フロケ周期境界条件

一般に，ソルバアルゴリズムではNeumannValue[g-q p(X) ,X∈Γ_b]を使って， が成り立つように境界 上の流束を指定する．

しかし音響モデルでは二重極音源 は領域内でしか指定できず，境界のどの部分上でも常にゼロに等しいため となる．

それぞれの境界条件について，特定の境界条件が時間調和解析，固有周波数解析，その両方のどれに適用されるのかを以下の方法で述べる．

インピーダンス境界条件

定式化

境界に指定されたインピーダンスがある場合，インピーダンス境界条件は以下で与えられる：

導出

音波が別の媒体に伝達したり，多孔質表面等の部分的に反射する境界に遭遇した場合，音波の一部がその接触面で反射し接触面全体に伝わる．

この入射波，反射波，伝達波の間の関係を定式化するために使われる特性の一つに，固有音響インピーダンスがある．固有音響インピーダンス は音波が媒体を移動するときの音の粒子速度に対する音圧の比であり，以下で定義される：

「時間領域における音響学」から，インピーダンス境界条件は指定された境界インピーダンス を使って以下のように定式化できることが分かっている：

調和波関係 を挿入し共通項 を因数分解すると，インピーダンス境界条件は下のようにして周波数領域で適用することができる：

ここで波動数 は音速に対する角周波数の比を示す．

インピーダンス境界条件は以下で使える．

時間調和解析におけるインピーダンス境界条件

インピーダンス境界条件を説明するために，次の例では右側が多孔質表面である円筒を考える．多孔質表面は部分的に反射する境界なので，インピーダンス境界として扱う．

音波は時間領域で調べた方が直感的である．解は調和波関係(39)()を使って時間領域に変換することができる．

インピーダンス境界は部分的に反射する境界をモデル化しているので，反射されるよりも多くのエネルギーが右側に移動している．これにより結果の波動は右に移動しているように見える．圧力振幅の最大値と最小値は空間で固定されている（破線）点に注意のこと．この種の波動は部分重複波と呼ばれる．

最大振幅と最小振幅の間の比は定在波比(SWR)と言われる．

定在波比は周波数に独立であることが示される．未知の境界の場合，領域の固有音響インピーダンス が与えられると，定在波比を使ってインピーダンス と反射係数 を測定することができる：

ここで と はそれぞれ反射波と入射波の振幅である．

吸収境界条件

定式化

特定のタイプの入力波および波の起点から境界,までの距離 がある場合，吸収境界条件は以下で与えられる：

• 平面波：

• 円筒波：

• 球面波：

導出

通常，無限にまで拡張されるシミュレーション領域は，シミュレーションにおいて実現可能なオプションとは言えない．吸収境界条件は無限領域をモデル化するのに使われる方法である．吸収境界条件(ABC)は到来波を吸収することで動作するため，モデルが無限範囲を持つかのように動作させることができる．シミュレーションを無限範囲でモデル化する方法は吸収境界条件だけではない．吸収境界条件の代りに完全整合層 (PML)を使うこともできる．

時間領域のチュートリアルから，吸収境界条件は以下で与えられることが分かる：

調和波関係 を挿入すると，周波数領域の吸収境界条件が与えられる：

平面波は ，円筒波は ，球面波 が使われる．

吸収境界条件は以下で使える．

時間調和解析における吸収境界条件

例として，計算領域が から に設定された無限長の円筒を見る．領域の計算をモデル化するために，右端に平面波吸収境界条件を加える．

結果の音圧場は単純に右移動の調和波なので，圧力振幅は全ての周波数領域でで固定される．

到来波は，シミュレーション領域が無限範囲を持つかのように，右側の境界で吸収される．

Sound Hard境界条件（壁）

定式化

壁境界 では，sound hard境界条件は以下で与えられる：

導出

sound hard境界では，前進運動ができないため，粒子速度の垂直成分 はゼロである：

(40)を運動量保存方程式 に代入し，調和波関係(41) を適用すると，sound hard境界条件は周波数領域で以下のように定式化することができる：

境界のどの部分にも境界条件が指定されていない場合は，デフォルトでノイマン境界条件が使われる．これは，与えられた境界において境界条件が指定されていない場合はsound hard境界がデフォルトの境界条件であることを意味する．

sound hard境界条件は以下のように使える．

時間調和解析におけるSound Hard境界条件

時間調和解析の例として，一端が閉じた円筒を見てみる．

振幅場の形はいくつかの最小点と最大点を示しており，それぞれ節，腹と呼ばれる．この時間領域では，節は定常波が変位を持たない位置であり，腹は最大の変位を持つ位置である．周波数領域では，変位は最小値，最大値として表される．sound hard境界では前進運動は不可能なので，音圧 は右端の最大値，つまり腹の一つで固定される．最大値は放射境界によって設定された振幅の2倍に相当する．

この結果は時間領域で調べるとより直感的である．

音波は右にも左にも動かず，時間とともに単に振動する．この種の波は定常波と言う．これは2つの波が反対方向に移動することによる重なりによって形成される．

この場合，右に向かう波は放射境界によって，左に向かう波はsound hard境界から反射された波によって生成される．進行波が定常波を与えるようにどのように重なり合うかに興味のある方は，時間領域における音響学をご参照いただきたい．

固有振動解析におけるSound Hard境界条件

時間調和解析とは異なり，固有振動解析は固有振動 および音響系の対応する固有モード（固有関数）を求めることを目的とする．次の例では，両端が閉じた円筒を考える．

最初の固有振動数は音のない解に対応し，ゼロである．したがって，固有振動数と固有モードの最初のペアは自明な解であり，固有モード0で表される．

すでに説明した通り，sound hard境界では前進運動は不可能なので，音圧はその最大値に固定される．

法線速度境界条件

定式化

境界上の指定された音の粒子速度が である場合，法線速度境界条件は以下で与えられる：

導出

非零の時間調和粒子速度が境界で指定される場合，この種の境界は法線速度境界と呼ばれる：

(42)を運動量保存方程式 に代入し調和波関係(43) を適用すると，法線速度境界条件は，周波数領域で次のように定式化できる．

法線速度境界条件は以下で使うことができる．

時間調和解析における法線速度境界条件

次の例では，右側の境界に調和振動が導入され，既知の速度振幅 で振動する．音場は次に示すようにNeumannValueを使って計算することができる．

結果の音波は，左に移動する調和波に過ぎないため，全ての周波数に対して領域全体で振幅分布が固定される．

法線速度境界は右端に調和波を生成し，それが左に伝搬する．

「時間領域における音響学」で説明してあるように，法線速度境界は，固有音響インピーダンス が既知のとき圧力源境界で置換できる．

Sound Soft境界条件

定式化

sound soft境界条件は以下で与えられる：

導出

sound soft境界では，媒体圧力は周囲の基準圧力に等しく設定される．つまり，であり境界の音圧はゼロで固定されているということである：

sound soft境界条件は次で使うことができる．

時間調和解析におけるSound Soft境界条件

時間調和解析の例として，一方だけが開いた円筒を考える．開いている端は，音波の運動に制約がないため，sound soft境界として扱われる．

のとき音圧は で最小になる．ここにはsound soft境界が位置しており，節と呼ばれる．音圧場の最大値は放射境界によって設定された振幅の2倍である．

sound hard境界条件と同様に，結果の音波は反対向きに移動する2つの波の重なりによって形成される定常波である．

この場合，右向きの波は放射境界によって生成され，左向きの波はsound soft境界から反射された波である．移動する波が定常波を与えるように重なる現象に興味のある方は，時間領域における音響学をご参照いただきたい．

固有振動解析におけるSound Soft境界条件

固有振動解析におけるsound soft境界の挙動を示すために，両端が開いた円筒を考える．

予想通り，sound soft境界の音圧はゼロで固定されている．

圧力源境界条件

定式化

境界T上の指定された圧力振幅が のとき，圧力源境界条件は以下で与えられる：

導出

境界で非零の時間調和音圧が指定されると，圧力源境界条件を考える．圧力源境界条件はDirichletConditionおよびNeumannValueを使って指定することができる：

音圧源境界条件は以下で使うことができる．

時間調和解析における圧力源境界条件

圧力源境界条件のディリクレ条件を定式化するためには，調和波関係(44) で方程式(45)を書き直すだけである：

ここで は圧力源の振幅を表す．

結果の波は左に移動する調和波に過ぎないので，振幅分布は領域を通してに固定される．

法線速度境界と同様に．圧力源は右端に調和波を生成し，それが左に伝搬する．

圧力源はNeumannValueでモデル化することもできる．「時間領域における音響学」で示す通り，圧力源に対するNeumannValueの設定は以下で与えられる．

調和波関係(46) を挿入すると，圧力源境界条件は周波数領域において以下のように定式化することができる．

上のアニメーションは，ディリクレモデルと同じ効果を示している．圧力源境界のノイマンモデルとディリクレモデルの間のトレードオフに興味のある方は，「時間領域における音響学」の中の対応するセクションをご参照いただきたい．

放射境界条件

定式化

境界 上の指定された入射音圧が ，波の方向ベクトルが のとき，放射境界条件は以下で与えられる：

境界の法線ベクトル ，波の方向ベクトル ，波の入射角 の関係は以下に示す通りである：

従って，方程式(47)は波の入射角 で以下のように表すことができる：

1D領域で放射境界を適用する場合，入射角 は常にゼロ（直角入射）であり，方程式(48) は次のように簡約される：

1Dの放射境界の例は次のセクションに挙げる．斜入射の反射音波を示す2Dの例はこちらに挙げる．

導出

放射境界は，圧力源境界と吸収境界の特性を組み合せたハイブリッドの境界条件である．入力波を生成するために時間調和の音圧 が境界で指定されるが，単純な圧力源とは異なり，放射境界では出力波はほとんど反射なしで計算領域を出ることができる．

「時間領域における音響学」で示しているように，放射境界条件は以下で与えられる：

調和波関係(49) を挿入すると，放射境界条件は周波数領域で以下のように定式化することができる．

放射境界条件は以下で使える．

時間調和解析における放射境界条件

次の例では，一端が閉じた半無限の円筒を考える．右端は調和音波が計算領域に入る放射境界条件として扱われ，左端はsound hard境界として設定される．

結果の音圧場は，入力波と出力波の両方の重なりによって形成される定常波を示す．この場合，における放射境界は左に移動する入力波を生成し，における暗示的なsound hard境界はその波を右側に出力波として反射する．これによりこの領域は、制約条件のない放射境界から離れる．進行波が重なって定常波になる過程に興味のある方は，「時間領域における音響学」をご参照いただきたい．

フロケ周期境界条件

定式化

音圧を音源境界から周期境界へマップする関数 がある場合，フロケ周期境界条件は次のように書くことができる：

ここで は周期境界から音源境界までのオフセット距離，は音波の方向を示す波数ベクトルである．

導出

フロケ周期境界条件は，空間的に周期的な領域内で時間調和音波をモデル化するために使われる．つまり，音源境界として参照される1つの境界から得られる情報は，マッピング関数 を使って，周期境界として参照される別の境界にマップすることができるのである．フロケの周期境界条件は，音響PEDモデルでPeriodicBoundaryConditionで設定される．

フロケの周期境界条件は以下で使える：

時間調和解析におけるフロケ周期境界条件

例として， 平面における以下の周期構造を考える．この構造は および 方向に沿って無限に拡張するものと想定される．周期境界条件を使うと，構造を切り取って単位セルをシミュレーション領域として取ることができる．

このモデルでは，音波は底辺から入射角 で入ってくる．これはとして表される2Dの放射境界でモデル化される．左端に周期境界，右端に音源境界を設定して，音波が領域の右側を通過するときに同じ大きさで左側に現れるようにする．

さらに, 境界を通過する出力波のシミュレーションを行う必要がある．この場合，出力角が非正規であるため吸収境界条件は使えない．その代りに，「完全整合層」という別の方法を使う．完全整合層についての詳細は完全整合層(PML)にある．

完全整合層を実装するために，完全整合層領域を含むようにシミュレーション領域を拡張する：

通常の領域と完全整合層領域の間の内部境界を維持する2D領域を設定するために，手動でメッシュを作成する必要がある．メッシュ生成についての詳細は「要素メッシュの生成」チュートリアルをご一読いただきたい．

実装した完全整合層でよい結果を得るために，メッシュ生成にデフォルトよりも細かい格子を使う．ここでは最大の辺長を に設定する．これは領域の （縦）方向と （横）方向に少なくとも10個の要素があることを意味する．

次に左端にフロケ周期境界条件を設定する．周期境界（左端）と音源境界（右端）の間のオフセット距離はである．

斜入射の入力音波のシミュレーションを行うために，領域の底辺に2Dの放射境界条件を適用する．

アニメーションの質を向上させる方法はこちらに記載してある．

このアニメーションは領域上の非正規な音波の伝播を示している．波が右の境界を通過するとき，空間周期性によって左側に再び現れる．音の媒体は音響的に無損失と想定されるため，波の大きさは常に一定レベルに保たれる．

完全整合層(PML)

完全整合層(PML)は無限範囲のシミュレーション領域をモデル化する方法である．したがって，PMLは吸収境界条件の代替方法と言える．以下のセクションでは，周波数領域のモデリングで使われるヘルムホルツの偏微分方程式に対するPMLの実装方法を説明する．完全な導出と説明は，「時間領域における音響学」の付録セクションに記載されている．

PMLを実装するためには，2つのことが必要になる．まずシミュレーション領域が拡大されなければならない．この拡張はPMLがアクティブとなる領域である．次に，偏微分方程式の座標変換が行われる．

PML座標変換[51]後の三次元ヘルムホルツ方程式は以下によって与えられる：

ここで ，， はPMLの吸収係数である．各次元のPML減衰を制御するために3つの補助変数 ，， を導入する．

である一次元の場合，方程式(52)は次のように簡約される：

吸収係数 は選択しなければならないチューニングパラメータであり，PML内でから に線形に増加するよう設定される．薄いPMLにはより大きいの値が要求されるが，大きいは数値的な反射を増加させる傾向がある．

PMLは以下で使える：

時間調和解析におけるPML

ここでは例として，右にからまでのPMLを持つ，からまでの範囲の計算領域を使って，から までの1D領域のシミュレーションを考える．

グラフから分かるように，もとの非PML領域では である．

非PML領域の音圧は入射波振幅 において固定されており，PML領域内でゼロに減衰する．PML界面における数値的反射を最小化するために，吸収係数 とPMLの幅は注意深く選ばなけらばならない．これら2つのパラメータの間のトレードオフは，「時間領域における音響学」に述べてある．

波はPML内で減衰するが，もとの領域内では変化しない．減衰率は周波数とは関係ない．

用語集

    

記号	解説	単位
ρ	媒体の密度	[kg/m3]
c	媒体内の音速	[m/s]
p	音圧	[Pa]
p	局所的な音の振幅	[Pa]
	指定された境界圧力	[Pa]
	音圧の共役	[Pa]
t	時間	[s]
tend	シミュレーション終了時間	[s]
X	位置ベクトル	[m]
s	方向変換	N/A
F	オプショナルの二重極源	[N/m3]
	二重極源の強さ	[N/m3]
θ	二重極指向角	[rad]
Q	オプショナルの単極源	[1/s2]
	単極源の強さ	[1/s2]
d	二重極源の分離距離	[m]
λ	音の波長	[m]
Ω	シミュレーション領域	[m]
k	波数	[rad/m]
f	音波の周波数	[Hz]
ω	音波の角振動	[rad/s]
δ	ディラックのデルタ関数	N/A
	正規化デルタ関数	N/A
γ	正規化パラメータ	[m]
h	メッシュ間隔	[m]
Xs	音源の位置	[m]
κ	体積弾性率	[Pa]
α	減衰因数	[m2/(s·N)]
ϕ	多孔性	N/A
Vv	空隙容量	[m3]
VT	合計容量	[m3]
Rf	流れ抵抗	[kg/(m3·s)]
β	ガウスパルスの標準偏差	[m]
ζ	音の粒子変位	[m]
v	音の粒子速度	[m/s]
	指定された境界速度	[m/s]
T	音波周期	[t]
Z	特性インピーダンス	[Pa·s/m]
Zb	境界インピーダンス	[Pa·s/m]
Ar	反射波の振幅	[Pa]
Ai	入射波の振幅	[Pa]
σ	PMLの吸収係数	[rad/(s·m)]
σmax	吸収係数の最大値	[rad/(s·m)]

参考文献

1.  Ihlenburg, Frank. The Medium-Frequency Range in Computational Acoustics: Practical and Numerical Aspects. Journal of Computational Acoustics, Vol.11, No. 2 175-193, 2003.

2.  Heutschi, Kurt. Lecture Notes on Acoustics I. Swiss Federal Institute of Technology Zurich, 2016.

3.  Johnson, Steven. Notes on Perfectly Matched Layers (PMLs). MIT, 2010.

4.  Bilbao, Stefan and Hamilton, Brian. Directional Source Modeling In Wave-Based Room Acoustics Simulation. IEEE, 2017.

5.  Peskin, Charles. The Immersed Boundary Method. Cambridge University, 2002.

6.  Russell, Daniel, Titlow, Joseph and Bemmen, Ya-Juan. Acoustic monopoles, dipoles and quadropoles: An experiment revisited. American Journal of Physics 67, 660, 1999.

7.  Vita, Micro. The Wave Equation with a Source. Oklahoma State University.

8.  J. De Moerloose and M. A. Stuchly, Behavior of Berenger's ABC for evanescent waves, IEEE Microwave and Guided Wave Letters, vol. 5, no. 10, pp. 344-346, Oct. 1995.

9.  J. Berenger, Evanescent waves in PML's: origin of the numerical reflection in wave-structure interaction problems, IEEE Transactions on Antennas and Propagation, vol. 47, no. 10, pp. 1497-1503, Oct. 1999.

10.  J. De Moerloose, Jan & Stuchly, Maria. Reflection analysis of PML ABCs for low-frequency applications, IEEE Microwave and Guided Wave Letters, vol. 6., no. 4, pp. 177-179, Apr. 1996.

11.  E. Turkel and A. Yefet. Absorbing PML boundary layers for wave-like equations, Applied Numerical Mathematics, vol. 27, pp. 533-557, 1998.

12.  G. Pan, A. Abubakar and T. Habashy. An effective perfectly matched layer design for acoustic fourth-order frequency-domain finite-difference scheme, Geophysical Journal International, vol. 188, pp. 211-222, 2012.

13.  T. J. Cox, P. D'Antonio. Acoustic absorbers and diffusers: Theory, design, and application, London: Spon Press, 2004.

14.  E. W. Weisstein. Fast Fourier Transform, MathWorld--A Wolfram Web Resource. Retrieve from: http://mathworld.wolfram.com/FastFourierTransform.html.

関連するテクニカルノート

▪
• 時間領域における音響学
▪
• 周波数領域における音響学
▪
• 音響学モデルコレクション
▪
• PDEModelsの概要