第 6 章 スペクトル法による軸対称空洞内の音場解析 97
7.2 数値計算例
7.2.1 剛体球の散乱音場
となる。Γ3以外の境界条件は粒子速度が0であるとすると，式(6.66)は，
−
∑M l=0
Dil(ξ)TwlwjΦ1m(ξl, ηj) +mwiwjΦ2m(ξi, ηj)−
∑M l=0
Djl(η)TwiwlΦ3m(ξi, ηl) +k2Ji,jwiwjϕm(ξi, ηj)
+ℓj(1) Ji,M
|Ji,M|
∑∞ n=0
anm(k, Rd)ψ(ξi,1)wi
∑M l=0
ψ(ξl,1)ϕm(ξl,1)wl =−Ji,jwiwjqm(ξi, ηj) (7.16) と書き換えることができる。上式をi= 0,1, . . . , M, j = 0,1, . . . , M について書き並べると，最終的に次 のようなマトリクス方程式を立てることができる。
[
K1+K2+K3−B−ω2
c2diag (J)W ]
ϕn= diag (J)W qn (7.17) 新しく導入したマトリクスBがDtN写像を表しており，その成分は
BM(M+1)+i,M(M+1)+j = Ji,M
|Ji,M|
∑∞ n=0
anm(k, Rd)ψ(ξi,1)wiψ(ξj,1)wj, i, j= 0,1, . . . , M
BM(M+1)+i,M(M+1)+j = 0, otherwise (7.18)
と表される。なお，上式はLegendre陪関数の無限級数であるが，実際の数値計算においては有限の次数 で打ち切る必要がある。その他のマトリクスについては前章で導出したものである。
となる。この微分方程式は特異なStrum-Liouville問題と呼ばれ，このとき解ϕは特定の境界条件を満た す必要がない [6]。つまり，Neumann境界条件を与えてしまうと図7.4(b)のように正しい解が得られな いものと考えられる。
そこで，ξ方向の節点として境界に節点をもたないLegendre-Gaussノード [7]を用いる。これにより，
境界条件を与えないということが可能になる。η方向は境界条件が必要であるため，通常のLGLノード を用いる。図7.5に仮想境界の半径Rd= 3.5，ξとη方向の節点数をともに20としたときの節点の配置 を示す。
Legendre-Gaussノードを用いたNGMによる音圧分布を図7.4(c)に示す。理論解析解と近い音圧分
布が得られていることが分かる。
仮想境界の半径とLegendre陪関数の打ち切り次数
ここではDtN写像の二つのパラメータ，仮想境界の半径とLegendre陪関数の打ち切り次数を変化さ せながら計算を行い，それらと計算精度の関係について調べる。仮想境界の半径を3，3.5，4 mとして，
それぞれ500 Hzの波長に少なくとも格子が二つ以上含まれるように離散化し，横軸にLegendre陪関数
の打ち切り次数，縦軸に理論解析解(打ち切り次数60)とのL2相対誤差をとったものを図7.6に示す。
仮想境界の半径が大きくなるほど，誤差が大きくなっていることが分かる。仮想境界の半径を大きく とると領域が広くなり，当然必要な節点数も大きくなるため，計算負荷という観点からも好ましくない。
よって仮想境界の半径はできるだけ小さくしたほうがよいと言える。
次にLegendre陪関数の打ち切り次数と誤差の関係を見てみると，打ち切り次数を大きくするにつれて
いったん誤差は小さくなるが，さらに大きくしていくと高い周波数から順に誤差が大きくなっている。こ
3.0m 2.0m
0.5m
Source
図7.3: 剛体球と音源の配置
0 1 2 3 4 5 6
x
[m]
−6
−4
−2 0 2 4 6
z
[m ]
0 20 40 60 80
Re lat ive SP L [ dB ]
(a)理論解析解
0 1 2 3 4 5 6
x
[m]
−6
−4
−2 0 2 4 6
z
[m ]
0 20 40 60 80
Re lat ive SP L [ dB ]
(b) NGM(Legendre-Gauss-Lobattoノード)
0 1 2 3 4 5 6
x
[m]
−6
−4
−2 0 2 4 6
z
[m ]
0 20 40 60 80
Re lat ive SP L [ dB ]
(c) NGM(Legendre-Gaussノード)
図7.4: θ方向に用いる節点の種類による音圧分布の違い
れは式(7.14)の積分を数値積分で処理した式(7.15)の誤差によるものだと考えられる。今，図7.3のθ
方向はLegendre-Gaussノードで離散化しており，用いる数値積分はLegendre-Gauss求積である。θ方 向の節点数をNθ とすると，この数値積分は被積分関数が2Nθ −1次の多項式までは厳密に計算できる が，それより高次の多項式の積分は近似となる。式(7.14)の被積分関数において，ϕmはNθ個の節点値 から内挿されるNθ−1次の多項式，また，ψ(ξ′,1)のうち，一般曲線座標系で表されたr(ξ′,1)もまた同 様にNθ−1次の多項式である。Legendre陪関数Pn|m|はn次の多項式であるから，被積分関数全体では 2Nθ+n−2次の多項式となっている。つまり，nが2以上の積分は近似となり，nが大きくなるにつれ て誤差も大きくなる。たとえ被積分関数の次数が2Nθ−1より大きくても，その空間的な変化が緩やか
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
x
−4
−3
−2
−1 0 1 2 3 4
z
図7.5: 剛体球の周りの格子点(Rd= 3.5，ξ，η方向の節点数がともに20のとき)
で，2Nθ−1次の多項式で近似できれば誤差は小さい。しかし，波数kが大きくなると解ϕmの空間的な 変化が大きくなり，数値積分の誤差も大きくなると予想される。
このことを確認するために，式(7.14)の積分を数値積分で計算する際，Legendre陪関数の次数によっ て被積分関数がどの程度近似されているのか調べた。図7.6(a)において，周波数を125 Hzと250 Hzと したとき，誤差の減少が横ばいになる前の次数の被積分関数をプロットしたものが図7.7である。青の実 線は理論解析解から求まる本来積分したい被積分関数，赤の破線はそれを43点のLegendre-Gaussノー ド上で離散化してLegendre多項式で内挿したもの，つまり数値積分によって積分される被積分関数を表 している。Legendre陪関数の次数を大きくしていくと，近似する被積分関数の誤差が大きくなっている。
また，Legendre陪関数の次数が同じ20次であっても，125 Hzより250 Hzの被積分関数のほうが誤差 が大きい。これらのことから，Legendre陪関数の次数，もしくは波数が大きくなると数値積分の誤差に よってDtN写像の精度が低下してしまう。つまり，スペクトル法における節点数との関係により，適切 な打ち切り次数があるということである。
数値積分の精度を上げる単純な方法は積分点を増やすことである。NGMの場合には積分点を節点とす るため，すなわち節点を増やすことが積分点を増やすことに相当する。そこで，r方向とθ方向の節点数 を倍にして図7.6(a)と同様の計算を行った。図7.8にその結果を示す。図7.6(a)では打ち切り次数を36 次以上にすると誤差が大きくなる傾向にあったが，図7.8ではそれが見られず，改善されている。以上の
0 10 20 30 40 50 60
Truncation number
10-5 10-4 10-3 10-2 10-1 100 101 102
L2
re lat ive er ror
63 Hz 125 Hz 250 Hz 500 Hz
(a)仮想境界の半径3 m (1849節点)
0 10 20 30 40 50 60
Truncation number
10-5 10-4 10-3 10-2 10-1 100 101 102
L2
re lat ive er ror
63 Hz 125 Hz 250 Hz 500 Hz
(b)仮想境界の半径3.5 m (2500節点)
0 10 20 30 40 50 60
Truncation number
10-5 10-4 10-3 10-2 10-1 100 101 102
L2
re lat ive er ror
63 Hz 125 Hz 250 Hz 500 Hz
(c) 仮想境界の半径4 m (3249節点)
図7.6: Legendre陪関数の打ち切り次数と誤差
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 θ
−100
−80
−60
−40
−20 0 20 40 60 80
Integrand
Exact Intepolated
(a) 125 Hz，20次
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 θ
−100
−80
−60
−40
−20 0 20 40 60 80
Integrand
Exact Intepolated
(b) 125 Hz，22次
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 θ
−100
−80
−60
−40
−20 0 20 40 60 80
Integrand
Exact Intepolated
(c) 125 Hz，24次
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 θ
−200
−150
−100
−50 0 50 100 150 200 250
Integrand
Exact Intepolated
(d) 250 Hz，16次
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 θ
−200
−150
−100
−50 0 50 100 150 200 250
Integrand
Exact Intepolated
(e) 250 Hz，18次
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 θ
−200
−150
−100
−50 0 50 100 150 200 250
Integrand
Exact Intepolated
(f) 250 Hz，20次
図7.7: DtN写像における被積分関数の実部
0 10 20 30 40 50 60
Truncation number
10-6 10-5 10-4 10-3 10-2 10-1 100 101 102
L2
re lat ive er ror
63 Hz 125 Hz 250 Hz 500 Hz
図7.8: Legendre陪関数の打ち切り次数と誤差(仮想境界の半径3 m，節点数7396)
結果から，提案手法におけるLegendre陪関数の打ち切り次数は，単純に大きくすれば精度が上がるわけ ではなく，境界上の節点数についても考慮する必要があるということが言える。
アドミッタンス終端との比較
領域型の解法で開領域音場を解く最も単純な方法はアドミッタンス終端 [50]と呼ばれる方法である。
この方法は無限領域を有限領域で打ち切り，境界条件として媒質の特性インピーダンスρcを与えるもの であり，無限遠におけるSommerfeldの放射条件を近似的に与えることに相当する。ここでは境界条件と してアドミッタンス終端を適用した場合と提案手法により計算した散乱音場の音圧分布を比較する。
図7.9と7.10に各手法で求めた250 Hzと500 Hzの音圧分布を示す。NGMの節点数は7396(500 Hz で1波長あたりの格子数4)，DtN写像の打ち切り次数は30とした。アドミッタンス終端では疑似反射の 影響によると思われる誤差が生じている。アドミッタンス終端は波面が仮想境界近傍で平面波に近くなる 場合に有効な方法である。そのため，仮想境界を音源や散乱体から遠くに設定すれば精度を改善できると 予想されるが，NGMによる有限領域を広くすることになり，それに応じて節点数も多くする必要がある。
それに対して提案手法では理論解析解に近い音圧分布が得られており，アドミッタンス終端に対する優位 性を確認できる。提案手法は境界で厳密解である理論解析解と連成しているため，原理的には仮想境界の 半径に制約はない。そのため，有限領域を小さく設定することにより，節点数を削減することができる。