閉空間内における音響シミュレーションに対する数値計算法及び境界条件の検討
理学専攻・情報科学コース

1

はじめに2.3

音は, 我々にとって身近な流体現象である. コンサー
トホールなどにおいては, 建設後に音響障害が発生することを防ぐため, 設計前に数値シミュレーションを行い, 音の響き方を予測しておく必要がある. 自動車や工場内で使用される機械など, 工業分野での騒音対策は今日では当然であり, 設計時に騒音の正確な予測をしておくことは重要である.
本研究では, 閉空間内の音の伝播や反射を, 数値シミュレーションにより視覚的および定量的にとらえ,最終的に実際のホール内の音響シミュレーションを行うことが目的である.
先行研究 [1] において, 一部の境界面から反射波が発生していることを改善するため, 吸収境界条件と時間発展の計算手法を変更した. その結果, 音源から近い境界面における数値誤差を解消することができたが, 音源から遠い面において反射波が発生した.
これを改善するため, 本研究では境界条件の計算モデルとして PML[5] を, 空間方向の計算手法として結合コンパクト差分を用いる.

2計算方法

2.1

室内空間のモデル化

室内の簡単なモデルとして直方体領域を考え, 床面を x − z 平面とし, 床面の長辺を x 軸, 短辺を z 軸, 鉛直方向を y 軸とする, 3 次元直交座標を用いた. 音速と, 短辺の長さを基準に計算領域は (2, 1, 1) とした. 格
子数は x, y, z 方向それぞれに 200, 100, 100 ととり, 等間隔格子を用いる. 音源は（x, y, z ）= (0.5, 0.42, 0.5)を中心とした球と考え, 1 波長分の振動をする音波の観測を行った. 境界面は吸収境界とした.

雲井 みのり解法

2.3.1 空間の離散化
結合コンパクト差分とは, 複数の高階微分の差分を同時に含んだ式を作り, 差分の精度に対応して複数の独立な式を解くことにより, 各階微分の差分近似値を求める手法である. 結合コンパクト差分式は,( ′)( ′′)′′′fi′ + a1 fi+1+ fi−1+ b1 ∆x fi+1− fi−1)( )2 ( ′′′)
d1 (′′′fi+1 − fi−1(3)+ fi−1=+c1 ∆x fi+1∆x()()
a2′′′′′fi′′ +f ′ − fi−1+ b2 fi+1+ fi−1
∆x i+1( ′′′))
d2 (′′′+c2 ∆x fi+1− fi−1
= ( )2 fi+1 − 2fi + fi−1 (4)∆x))
a3 ( ′b3 ( ′′′′′++ fi−1− fi−1fi′′′ + ( )2 fi+1f∆x i+1∆x( ′′′))
d3 (′′′+c3 fi+1+ fi−1
= ( )3 fi+1 − fi−1∆x(5)であり, この式の係数を, Nihei,Ishii[2] に基づき,
16d3 + 255
4d3 + 45
8d3 + 195
, b1 = −, c1 =
,
240
1200
1800
8d3 + 315,d1 =
240
11d2 − 15
3d2 − 7d2 − 3a2 =, b2 = −, c2 =,
1616488d3 + 154d3 + 15
a3 = d3 , b3 = −, c3 =,2060()− 165d2 − 450
d2 = 9.12992, d3 =28d2 − 80
a1 =とし, 1 階微分を 8 次精度で近似した.

2.3.2 時間の離散化
時間発展については, 高解像度多段時間法で計算した. Tam[3] に従い,2.2 基礎方程式
基礎方程式として, 圧縮性非粘性流れに対する音速,静止密度, 代表長さを基準とした無次元線形オイラー方程式を用いた.

u(n+1) = u(n) + ∆t

3
∑
j=0

bj (∂u (n−j)
)
∂t

(6)

b0 = 2.305580888383, b1 = −2.4910075998482,
b2 = 1.5743409331815, b3 = −0.3858914221716

∂ρ∂t
∂V
∂t=

⃗ +Q+ρ−∇ • V
P ML(1)のようにとることで, 3 次精度となるが高い周波数で=⃗P M L−∇ρ + V(2)

もより正しい計算ができるようにした.

⃗ は流速,t はここで,ρ は空気の静止密度からのずれ, V時間である. 式 (1) の Q は音源項を表すが, 本研究では音源から発せられる波から数値誤差により余分な波が発生しないように, 空間的に数点を使ったガウス分布の振幅を時間的に振動させることによってなめらか⃗P M L は後述の PML,Vな音源を再現した. なお, ρP ML

境界を用いた場合に生じる項である. (PML を用いないときは 0 である).

2.4

境界条件

2.4.1 球面波近似
完全吸収境界条件を正確に課すのは困難である. 先行研究では波の進行方向に対する 1 次元波動方程式を近似した特性方程式 [4] を利用した. たとえば,x 軸に負方向に進む波に対しては, 音速を 1 とした時,ρ+u=0(7)あるいは∂ )−ρ|x=0 = 0∂x∂t()∂∂−+u|x=0 = 0∂x∂tv, w = 0( ∂(8)(9)(10)

となる.
しかし,3 次元点音源から発せられる音は, 音源からの距離を r として外に広がる球面波
ρ=f (r − t)r(11)として表される. このとき, 境界条件を 1 次元で近似して計算すると, 波面が 3 次元的にならず境界付近で数値誤差として反射波が発生した. そこで 3 次元性を考慮するため, 音源を中心とした極座標で近似した特性方程式を利用した. これを球面波近似と呼ぶ.

2.4.2

Perfectly Matched Layer(PML)⃗P M L は以下のもの
基礎方程式 (1) , (2) の ρP M L , Vを用いる.

+∫∫ ∫m4 + n4 ρdt + n5ρdtdt∫∫ ∫∂∂udt + n1udtdt
m1∂x∂x∫∫ ∫∂∂m2vdt + n2vdtdt
∂y∂y∫
∫ ∫∂
∂m3
wdt + n3wdtdt
∂z∂z
∫
∫ ∫
⃗ + n4 V⃗ dt + n5⃗ dtdtm4 VV∫∫ ∫m1 ∇ ρdt + n1 ∇ρdtdt

m1=

σy + σz , m2 = σz + σx , m3 = σx + σy ,m4=

σx + σy + σz

ρP M L=+++VP M L

=(12)

(13)

(a) the high resolution (b) the eighth-order comfourth-order central diﬀer- bined compact diﬀerence
withboundary
ence method with boundary scheme
conditionofspherical condition of PML(Case2)
approximation(Case1)

Fig. 1: Density distribution on the plane of z=0.5 at
t=0.8(a) the high resolution (b) the eighth-order comfourth-order central diﬀer- bined compact diﬀerence
withboundary
ence method with boundary scheme
conditionofspherical condition of PML(Case2)
approximation(Case1)

Fig. 2: Density distribution on the plane of z=0.5 at
t=1.2

n1

=σy σz , n2 = σz σx , n3 = σx σy ,n4=n 1 + n 2 + n 3 , n 5 = σx σy σzある. 比較手法 (2a) では反射波どうしが干渉し, 直接

ここで, σx , σy , σz は吸収率であり, 本研究では Tam[3]に従い 40 と設定した. また, 本研究では全境界に短辺音程度の大きさの誤差が発生している. 一方, 本研究(2b) では, 反射波どうしの干渉によって大きな数値誤の長さの 1/10 の吸収層を設定した.

差はないが, 反射波によって直接音に乱れが発生している.

2.5

比較するケース

比較する手法は, 基礎方程式に本研究と同じく線形オイラー方程式を用い, 時間発展についても同じ
く 3 次精度の高解像度多段時間法を用いた. 空間方向の解法については, 高解像度 4 次精度中心差分法を用い, 境界については, 吸収境界に対して球面波近似で計算を行った. 本研究との相違点を以下の表に示す

空間方向
境界条件

34

まとめ

基礎方程式として無次元線形オイラー方程式を用い,その解法として, 時間発展に高解像度多段時間法を用い, 空間方向には 8 次精度の結合コンパクト差分を用いて計算を行った. 吸収境界において, 球面波近似を利用して計算すると, 反射波が発生する境界面があったので, 計算モデルとして PML を採用することで, 吸収境界面で反射がおさえることができた.

Case1(比較手法)

Case2(本研究)参考文献

高解像度 4 次精度中8 次精度結合コンパ心差分

クト差分[1] 雲井, 石井, 河村, 永田, “高精度・高解像度音響シミュレーションにおける境界条件の検討”, 第 30 回数値流体シンポジウム
予稿集 (2016)

球面波近似

PML

結果と考察

[2] Tomonori Nihei, Katsuya Ishii, “fast solver of the shallow
water equations on a sphere using a combined compact
diﬀerence scheme”, Journal of Computational Physics
187, 639-659(2003)

Fig.1 は t = 0.8 における z = 0.5 の面の密度分布である. 赤いほど密度が高く, 青いほど低くなるように[3] Tam, C. K. W. T., ”Computational Aeroacoustics: A
Wave Number Approach(Cambridge Aerospace Series),”Cambridge University Press.(2012)表示している. 比較手法 (1a) では吸収境界であるにも[4] Engquist, B. and Majda, A., ”Absorbing boundary conditions for the numerical simulation of waves,” Math.
Comp. 31,629-651(1977)

関わらず, 上部境界から反射波が少し発生しているが,本研究の手法 (1b) では反射波が見られない.

Fig.2 は t = 1.2 における z = 0.5 の面の密度分布で[5] Berenger, J-P, “A Perfectly Matched Layer for the Absorption of Electomagnetic Waves” Journal of Computational Physics 114, 185-200(1994)