第 8 章 スペクトル法による膜鳴楽器の音響振動連成解析 145
8.2 数値計算例
が成り立つ。これより薄板とシェルへの加振圧力ベクトルfpとfsはTp とTsを用いて fp =−TTp
{fs
Ms
}
=−TTpλ (8.37a)
fs =TTs {fs
Ms
}
=TTs λ (8.37b)
と表される。ここで，λ={fs Ms}Tとおいた。
薄板振動場，シェル振動場，音場の連立方程式
ここまで立ててきた連成条件(8.27)，(8.29a)，(8.29b)，(8.31)，(8.33a)，(8.33b)，(8.37a)，(8.37b) を式(8.9)，(8.10)，(8.11)，(8.12)に代入し，薄板とシェルの束縛条件式(8.35)とともに書き並べると，
[Kp−ω2Mp
]ζ =f −TTpλ+Ginϕin−Goutϕout (8.38) [Ks−ω2Ms
]d=WsTTsλ+WsHinϕin−WsHoutϕout (8.39) [Kin−ω2Min
]ϕin=−WinPinζ−WinQind (8.40) [Kout−ω2Mout−B]
ϕout=WoutPoutζ+WoutQoutd (8.41)
Tpζ =Tsd (8.42)
となり，さらに未知数を左辺に移項し，すべての方程式をマトリクス方程式の形で書くと，
Kp−ω2Mp O TTp −Gin Gout
O Ks−ω2Ms −WsTTs −WinHin WoutHout
Tp −Ts O O O
WinPin WinQin O Kin−ω2Min O
−WoutPout −WoutQout O O Kout−ω2Mout−B
ζ d λ ϕin ϕout
=
f 0 0 0 0
(8.43) が得られる。これが薄板振動場，シェル振動場，内部および外部音場を連成した方程式である。薄板振動 場への加振圧力f を与えて上式を解けば，薄板振動場，シェル振動場の振動変位，そして内部，外部音場 の速度ポテンシャルが同時に得られる。
0.29m
0.46m
Rigid
Membrane
S(0.2, 0, 0) R1(1, 0, 1)
R2(0, 1, 1) R3(0, 0, 1)
図8.5: 計算したティンパニの加振点と受音点の配置
数，ルジャンドル級数の打ち切り次数はともに10次とした。Coupled NGMの自由度は310とし，この ときの格子点を図8.6に示しておく。Coupled BEMは膜振動場の500 Hzの波長に対して要素寸法比が 約0.2となるように576個の要素に分割した。BEMは要素寸法比が波長の1/5以下となるように設定す ることが推奨されており [47]，それより高い周波数では十分な精度は保てないものと予想される。なお，
参照解には膜振動場の1000 Hzの波長に対して要素寸法比が約0.2となるように2226個の要素に分割し たCoupled BEMを用いた。
二つの手法で計算した各受音点における周波数応答関数を図8.7から8.9に示す。Coupled BEMは波 長に対する要素寸法比が1/5より大きくなる500 Hzあたりから参照解との誤差を生じていることが分か る。それに対してCoupled NGMはそれよりも少ない自由度でありながら高い周波数まで参照解との一 致が見られる。
受音点ごとに見ていくと，Coupled BEMには見られるがCoupled NGMには見られないピークがい くつかあることが分かる。例えば，受音点R2における432 Hzに二つのCoupled BEMでは小さなピー クがあるが，Coupled NGMではそれがない。膜の円周方向のモードの次数をn，半径方向のモードの次 数をmとして，それらの組み合わせを(n, m)と書くとすると，この432 Hzは(3, 1)モードに相当し，
図8.5のy軸を節線の一つとして線対称な振動モードである。そのため，この節線上にある受音点R2で は振動変位や音圧が極小となる。これは(1, 1)モードについても同様であり，受音点R1の219 Hzにあ るピークが受音点R2ではいずれの計算結果においても見られない。Coupled BEMは非対称に分割され
0.0 0.1 0.2 0.3 0.4 0.5 0.6
x
−0.6
−0.4
−0.2 0.0 0.2 0.4 0.6
z
図8.6: Coupled NGMに用いた解析領域内の格子点
た要素を用いており，(1, 1)モードでは波長が十分に長く，影響しなかったが，(3, 1)モードはこの非対 称に分割された要素によって厳密に対称な振動にならず，音圧が極小になっていないと考えられる。それ に対して，Coupled NGMは円周方向をフーリエ級数展開しているため，この節線上では理論的に振動変 位や音圧が0となる。また，受音点R3は軸上の点であり，この点では軸対称モードしか観測されないは ずである。実際に受音点R1やR2で観測されているほとんどの非軸対称モードのピークが受音点R3で は見られない。しかし，R3においてCoupled BEMのみ(4, 1)モードに相当する533Hzにピークが見 られる。これについてもBEMの非対称な要素分割が影響しているものと考えられる。
これらのピークが非対称な要素分割による影響であるとすれば，より細かく要素分割して計算精度を改 善することでこれらのピークは各モードの節において観測されなくなると予想される。そこで，このこと についてより詳細に調べるため，Coupled BEMの要素数を変化させながら周波数応答関数を計算した。
図8.10はCoupled BEMの要素数を576，2226，4746と変化させたときの受音点R3における533 Hz 付近の周波数応答関数を示している。要素数を増やしていくとピークの大きさと幅が小さくなっていくこ とが分かる。つまり，このピークは要素分割の影響によって生じたものであると言える。
次に，周波数応答関数の計算に要した計算コストを表8.1にまとめる。計算時間は1000/512 Hzおき
に512 tapsの計算に要した時間を表している。計算環境と連立一次方程式の求解に用いた関数は次の通
りである。
0 200 400 600 800 1000
Frequency [Hz]
−30
−20
−10 0 10 20 30
R el at iv e SP L [d B ]
NGM BEM Reference
図8.7: Coupled NGMとCoupled BEMにより計算したティンパニの周波数応答関数(受音点R1)
0 200 400 600 800 1000
Frequency [Hz]
−30
−20
−10 0 10 20 30
R el at iv e SP L [d B ]
NGM BEM Reference
図8.8: Coupled NGMとCoupled BEMにより計算したティンパニの周波数応答関数(受音点R2)
0 200 400 600 800 1000
Frequency [Hz]
−30
−20
−10 0 10 20 30
R el at iv e SP L [d B ]
NGM BEM Reference
図8.9: Coupled NGMとCoupled BEMにより計算したティンパニの周波数応答関数(受音点R3)
500 510 520 530 540 550 560
Frequency [Hz]
−30
−20
−10 0 10 20 30
Relative SPL [dB]
BEM (576) BEM (2226) BEM (4746)
図8.10: Coupled BEMの要素数を変化させて計算したティンパニの周波数応答関数(受音点R3)
表8.1: 周波数応答関数の計算に要した各手法の計算コスト
自由度 必要メモリ容量 [MB] 計算時間[sec]
Coupled NGM 310 1.5 61
Coupled BEM 576 5.3 593
Coupled BEM with OpenMP 576 5.3 242
• CPU: Intel Xeon 3.30 GHz
• Memory: 64 GB
• Library: NAG C Library Mark 23
• Subroutine: nag complex gen lin solve
図8.7から8.9の周波数応答関数から分かるようにCoupled BEMは高い周波数での精度を犠牲にしてい るにもかかわらず，Coupled NGMはその約1/10の時間で計算できている。同等の精度で比較すれば，
より顕著な差があらわれる可能性もある。Coupled NGMは一度係数マトリクスを構築してしまえば，そ のほとんどの成分を異なる周波数の係数マトリクスにも流用することができる。一方，Coupled BEMは 周波数が変わるたびに要素内の積分の計算が必要となり，係数マトリクスの構築に計算時間がかかってい る。この計算は並列計算プラットフォームであるOpenMPを用いて簡単に並列化が可能であり，表8.1
にはOpenMPを用いて並列化したCoupled BEMの計算時間も記している。計算時間を並列化していな
い場合の半分以下に短縮できているが，Coupled NGMのほうが高速である。なお，Coupled NGMのマ トリクス方程式(8.43)の係数マトリクスのうち，連成の過程で構築したマトリクスの成分はほぼ0であ り，係数マトリクス全体も疎行列と見なすことができる。従って，連立一次方程式の解法に疎行列向けの 直接法や反復法を用いることでさらに高速化できる可能性がある。