数値計算例

とおいたベクトルである。式(6.44a)-(6.44c)^{を使って式}(6.67)^のΦ¹_n^，Φ²_n^，Φ³_n^をϕ_n ^{で表せば，最終}

的に [

K1+K2+K3+jωC− ω²

c²diag(J)W ]

ϕ_n= diag(J)W q_n+un (6.71) というマトリクス方程式が得られる。ここで，

K1= [Iη⊗Dξ]^TW [

diag

(R²(Z_η²+R²_η) J

)

Iη⊗Dξ−diag

(R²(ZηZξ+RηRξ) J

)

Dη⊗Iξ

]

(6.72a) K2=Wdiag

(n²(RξZη−RηZξ)² J

)

(6.72b) K3= [Dη⊗Iξ]^TW

[

−diag

(R²(ZξZη+RξRη) J

)

Iη⊗Dξ+ diag

(R²(Z_ξ²+R²_ξ) J

)

Dη⊗Iξ

]

(6.72c) とおいた。NGMにおいては粒子速度境界条件はベクトルun，インピーダンス境界条件はマトリクスC として考慮されるため，SCMで行ったような方程式の置き換えなどは必要なく，式(6.71)^{をそのまま} ϕ_n^{について解けばよい。}

2m 1m

S

R1

R2 R3

図6.3: ^{計算した円筒形空洞}

6.2.2 円筒形空洞の強制振動問題

次に，境界条件としてインピーダンス境界をもつ円筒形空洞の音圧分布，および周波数応答関数を提案 手法と軸対称1次三角形アイソパラメトリック要素を用いたFEMによって計算し，精度を評価する。円 筒形空洞の寸法は図6.3の通りであるが，すべての壁面の比音響インピーダンスを15000 Pa·s/m(^垂直 入射吸音率で約0.1)^{とした。音源の位置}S^は，(rs, θs, zs) = (0.9,0,1.9)^とし，θ= 0^の(r, z)^平面上でr 方向51^点，z^方向101^点，計5151点の格子上で音圧を計算し，参照解とのL2相対誤差を算出した。な お，壁面がインピーダンス境界である場合の円筒形空洞内の音場については理論解析解はないため，最大 要素長0.02 m^，1 kHzの波長に対する要素寸法比0.06^，節点数10585の十分細かく要素分割したFEM による解を参照解とした。各手法におけるθ方向のフーリエ級数展開の打ち切り数は10^{次とした。}

自由度の比較

図6.5に横軸を未知数の数としたときの各周波数における音圧分布のL2相対誤差を示す。第4^章や第 5章の結果とは異なり，提案手法とFEMの差が小さく，提案手法は誤差が頭打ちになっているようにも 見える。この要因の一つとして，参照解に厳密解ではなく，FEMによる数値解を用いているためである と考えられる。500 Hz^や1000 Hzでは未知数の数に対する誤差がほぼ変化していないが，提案手法が参 照解の精度を上回っている可能性もある。

必要メモリ容量の比較

次に横軸を連立方程式の求解に必要なメモリ容量としたときの L2相対誤差の変化を図6.6に示す。

FEMの係数マトリクスは対称疎行列であるので，非零成分のみを圧縮行格納形式 [44]^{で格納したときの} メモリ容量を示している。一方，提案手法の係数マトリクスは非対称密行列である。直接法で連立方程式 を解く場合，方程式Ax=b^{のように未知ベクトル}x^に対して1回だけ作用する係数マトリクスA^を構 築する必要があり，特に二次元や三次元の問題のような大規模なマトリクスはそれを格納するメモリ容 量も大きくなる。もう一つの連立方程式の解法である反復法の場合には，近似解x^{(i)に対して係数マト} リクスAを作用させ，その結果が既知のベクトルb^{に近づくように近似解}x⁽ⁱ⁾を更新していく。このと き，係数マトリクスA^が

A=A1(A2+A3) (6.74)

のように分解できれば，マトリクスA^{とベクトル}x^(i)の積は

Ax⁽ⁱ⁾=A1(A2x⁽ⁱ⁾+A3x⁽ⁱ⁾) =A1y⁽ⁱ⁾, y⁽ⁱ⁾=A2x⁽ⁱ⁾+A3x⁽ⁱ⁾ (6.75) でも演算することができ，分解したマトリクスA1，A2，A3のそれぞれが疎行列であれば，メモリ容量 と演算量を大きく削減することができる。スペクトル法における係数マトリクスもそのような例の一つで あり，式(6.72a)^{のマトリクス}K1は密行列であるが，それを構成するIη⊗DξやDη⊗Iξ は疎行列と なっている。そのため反復法を用いる場合には，一つの密行列ではなく，複数の疎行列の形で格納するこ とで，メモリ容量を大きく削減することが可能である。本研究においては連列方程式の解法に直接法を用 い，その結果を図6.5^にはNGM (direct)として示しているが，反復法で直接法と同等の精度が得られる と仮定した場合の結果についてもNGM (iterative)として同図に示している。

提案手法において連立方程式の解法に直接法を用いる場合には，FEMよりも多くのメモリ容量が必要 になり，63 Hz^，125 Hz^，250 Hz^でL2相対誤差を0.01以下に抑えようとするとFEM^の約100^倍のメ モリ容量が必要である。一方，連立方程式の解法に反復法を用いる場合には，FEM^{と同等，もしくはそ} れ以下のメモリ容量で同等の精度が得られる。特に500 Hz^や1000 Hzにおいては誤差によってはFEM の1/50程度に削減することも可能である。ただし，図6.5^のNGM (iterative)の結果は直接法によって 得られたものであり，反復法は得られる解の精度がアルゴリズムや反復回数などのパラメータに依存する ため，今後それらを考慮したより詳細な検証が必要である。

音圧分布の比較

L2相対誤差は音圧分布の誤差を一つの値として評価できるため，誤差の収束の傾向を調べる際に便利 である。しかし，直観的に実際の音圧分布における違いと対応付けることは難しい。そこで，同程度の未 知数，および必要メモリ容量を用いる提案手法とFEMによって計算した音圧分布についても示してお く。特に差が顕著であった500 Hzと1000 Hzのθ= 0の(r, z)平面上における音圧分布を図6.7と6.8 に示す。提案手法における自由度は図6.5(e)において参照解との誤差がほぼ収束している288^とし，こ のとき必要メモリ容量は連立方程式の解法に反復法を用いるとすると約55 kB^{であった。}FEM^の自由度 はそれよりも大きい561とし，必要メモリ容量は約65 kBであった。なお，参照解としたFEM^の自由 度は10585^{，必要メモリ容量は約}1.2 MB^である。

提案手法は参照解とほぼ差がない音圧分布が得られている。それに対してFEM^{は提案手法よりも多い} 自由度，メモリ容量を用いているにもかかわらず，計算された音圧分布は参照解と大きく異なっている。

このように実際の音圧分布で比較しても，提案手法がFEMより精度がよいことが分かる。

周波数応答関数の比較

手案手法とFEMにより計算した円筒形空洞の周波数応答関数についても示しておく。受音点は図6.3 に示したR1(0.9,0,0.1)^，R2(0.9, π/2,0.1)^，R3(0.9, π,0.1)^の3点とした。計算に用いた提案手法の自由 度は128^{，必要メモリ容量は約}25 kB^，FEM^{の自由度は}231^{，必要メモリ容量は約}27 kB^{である。各手} 法におけるθ方向のフーリエ級数展開の打ち切り数はこれまでと同様に10^{次とした。}

図6.9^から6.11に各受音点における周波数応答関数を示す。提案手法により計算した周波数応答関数 は参照解と一致している。それに対して必要メモリ容量が提案手法と同程度のFEM^{は，低い周波数では} 参照解と合っているが，周波数が高くなるにつれて誤差が生じていることが分かる。

6.2.3 一般的な軸対称空洞の強制振動問題

最後に円筒形以外の軸対称空洞としてティンパニのケトル内部の周波数応答関数と音圧分布を提案手法 で計算する。ケトルの形状と音源および周波数応答関数を求めた受音点の位置を図6.13^{に示す。音源の} 位置S^{と受音点の位置}R1^，R2^，R3^は

• S: (r(ξ= 0.9, η= 0.9), θ, z(ξ = 0.9, η= 0.9)) = (0.27,0,−0.024)

• R1: (r(ξ= 0.1, η=−0.9), θ, z(ξ= 0.1, η =−0.9)) = (0.017,0,−0.43)

• R2: (r(ξ= 0.9, η= 0.9), θ, z(ξ = 0.9, η= 0.9)) = (0.27, π/2,−0.024)

• R3: (r(ξ= 0.9, η= 0.9), θ, z(ξ = 0.9, η= 0.9)) = (0.27, π,−0.024)

とした。壁面の境界条件は音響的に剛とした。提案手法の自由度は50^{，比較する}FEM^{の自由度は}92 とし，このときどちらもメモリ容量は約10kBであった。フーリエ級数展開の打ち切り次数は10^とし た。また，参照解として1kHzの波長に対して要素寸法比0.05^，自由度1470^のFEM^{でも同様の計算を} 行った。

提案手法に用いた物理座標系における格子点の分布を図6.13^{に示す。実際には第}4^{章の要領で図の}x 座標が正の格子点のみを用いて計算している。

図6.14^から6.16^{は提案手法と}FEMにより計算した周波数応答関数を示している。これまでの結果と 同様に，提案手法にって計算した周波数応答関数は参照解と一致しているのに対して，それより多い自由 度，同程度のメモリ容量を用いたFEMは周波数が高くなるにつれて参照解からずれている。

最後に提案手法とFEM^{で計算した}1000 Hz^{の音圧分布を図}6.17に示す。提案手法は参照解にかなり 近い音圧分布が得られている。FEMによって計算した音圧分布も概形は参照解に近い結果が得られてい るが，音圧レベルが極小となる節の位置に着目してみると，わずかではあるが参照解からずれていること が分かる。このように一般曲線座標系を用いてティンパニのケトルのような一般的な軸対称音場を計算す る場合においても，提案手法が有効であるということが分かる。