JP5024613B2 - 音場解析装置 - Google Patents

Description

本発明は、任意散乱体まわりの音場解析を行う音場解析装置に関する。

任意の解析条件に対して解析的に音場を求めるのは不可能であって、ＦＥＭ（有限要素法）、ＢＥＭ（境界要素法）などの数値解析法が一般に用いられる。これらは領域または境界を離散化して音場を近似的に求める解析法である。ＢＥＭは自由境界条件（放射条件）を設定し易く、境界のみの離散化で済むため外部問題に多く適用されている。ただし、境界の要素分割は粗くても波長の１／５程度は必要とされる。仮に波長の１／５程度の要素分割であっても、高周波数を扱う音場解析では要素の総数は膨大となる。一般に要素数の増加は計算時間の増大を招く。計算機の記憶容量は限られるため、解析対象の周波数はおのずと中低域に限定される場合が少なくない。

ＢＥＭやＦＥＭなどの数値解析法に対し、離散化を要しない解析法の一つに固有関数展開法が挙げられる。これは音場の一般解を固有関数で展開し、その展開係数を境界条件に適合するように決定する方法である。２次元外部問題において円形散乱体まわりの音場は固有関数展開を用いて解析解が求められる。また境界形状が球で与えられる３次元問題についても解析解が既に得られている。このほかにも、楕円柱の振動によって生じる放射音場をＭａｔｈｉｅｕ関数で展開し、境界条件から展開係数を求めることによって解析解を導いたり、回転楕円体表面上にある矩形ピストンの放射インピーダンスを回転楕円座標系に適合する固有波動関数で展開して解く方法が知られている。また、two-variable expansions methodによって有限長円筒まわりの散乱音場を求める方法や無限長円柱、半無限長円柱などの単純な境界形状の場合も、境界条件に適合した座標系と固有関数を導入することによって解析解が得る方法も知られている。

解析解は境界形状が単純な場合など限られた条件下でしか得られないが、ＦＥＭやＢＥＭに比べて解の物理的な考察が容易であり、計算機の記憶容量が少なくて済むなどの利点がある。そのため、その解析結果は数値解析法の精度検証やベンチマークなどにしばしば用いられる。また、これらの解析法はスピーカの音響機器の設計などにも広く適用されている。

しかしながら、任意の境界条件に対して固有関数と展開係数の双方を陽な形式で記述できる解を導くのは非常に困難である。その場合、固有関数は陽な形式で与え、展開係数は何らかの数値解析的な方法によって求める方法が採用される。たとえば、有限円筒管の音響透過問題を境界条件式から導いた展開係数に関する連立１次方程式に帰着させて解く方法が知られている。

またＴ−ｍａｔｒｉｘ法（ＴＭＭ）をはじめとして音場の積分方程式に基づく解析例も多い。これは、音場の積分方程式に速度ポテンシャルの固有関数展開形を代入して各次の積分方程式を導き、これらを連立させて展開係数を求める解析法である。ＴＭＭは、音場の積分方程式に基づくことから任意の境界形状が扱えるという利点があり、ＴＭＭと従来のＨｅｌｍｈｏｌｔｚ境界積分方程式法（ＨＩＥＭ）を比較し、その特徴を整理した例もある。この解析法は音響分野、電磁分野などで適用例は多く、たとえば円盤による散乱音場解析、有限長円筒による散乱音場解析、車室内の音場解析などが挙げられる。

また、振動する有限長円筒から放射される音場の解析手法として、音場の固有関数展開係数を境界条件にできる限り適合するように最小二乗法を用いて求める手法を定式化した例もある。また、この手法とほぼ同様の手法を用いて、球面波が入射する場合の有限長円筒による散乱音場を解析した例もある。この解析法は境界積分方程式を用いないため、ＢＥＭで外部問題を扱う際に問題となる解の非一意性を回避できると考えられ、最小二乗法を用いて展開係数を求める方法、Ｋｉｒｃｈｈｏｆｆ−Ｈｕｙｇｅｎｓの積分公式に基づいて展開係数を求める解法、およびＢＥＭによる音場解析結果を比較し、解の非一意性回避について論じられている。その結果から、最小二乗法を用いる方法では解の一意性が確保されることが示唆されている。

上記の解析法の多くは任意形状の境界を扱えるように定式化されてはいるが、実際の解析例は円筒形や円形などの軸対称の境界を対象としたものが多い。散乱体の境界の軸対称性を考慮することによって計算コストを低減させることができ、物理的な考察がより容易になるからである。

ただし、散乱体の軸対称性を利用した解析は固有関数展開に基づくものばかりでなく、ＢＥＭを適用した例も見られる。この場合、任意の軸対称散乱体に関する境界積分は子午線（回転体と回転軸を合む平面との交線）上の線積分に変換され、積分計算が効率的になることが示されている。また２次元等角写像法（conformal mapping method）を用いて軸対称散乱体による遠方音場を解析した例も知られている。

一方、固有関数展開に基づいて任意断面形状の境界を扱った解析例として、２次元任意形状の連続体の固有値解析がある。これは、固有関数展開形で表した波動関数が境界上に離散的配置された節点上で近似的に境界条件を満たすものとして、固有値マトリックス方程式を解く問題に帰着させる方法である。これは選点法と呼ばれる。また同様の手法を３次元問題に拡張し、選点法を用いて音場を解析した例もある。この解析法は、領域形状が凸形または滑らかな凹形ならば良好な結果が得られる。ただし、境界上の節点の配置法に明確な根拠は示されていない。数値解析上の不安定さを回避するため不規則な揺らぎを与えるなど、節点配置は経験的に定められている。また節点の密度について「１波長あたり数個の点は必要であろう」と述べられていることから、必要な計算機の記憶容量はＢＥＭと同程度になると考えられる。

また、固有関数展開法を出発点とする境界展開法による波動解析も提案されている（例えば非特許文献１、２参照）。これは、水中における線対称または軸対称の物体まわりの散乱波解析を扱ったものである。解析周波数は低く、建築音響分野では用いない境界条件を設定していることなどから、この解析法をそのまま建築音響分野に適用することは難しい。また級数和の収束性などの解析的な検討も十分ではない。しかしながら、解析法の適用範囲を境界条件の複雑さ（たとえば、散乱体形状の複雑さ）に結びつけて直感的に理解し易く、物理的な考察も容易であるという利点がある。
清川哲志、小林浩、日野幹雄，軸対象構造物による波の散乱と波力，土木学会論文報告集，Ｎｏ．３２１，ｐｐ．１０３−１１２，１９８２ 清川哲志、小林浩，面対称柱体の水中振動による付加質量特性の研究，土木学会論文報告集，Ｎｏ．３２１，ｐｐ．７９−９０，１９８２

一般に音場解析には有限要素法（ＦＥＭ）や境界要素法（ＢＥＭ）などの数値解析法が用いられるが、これらの数値解析法は解析領域や境界を微小要素に分割（離散化）して数値計算するため、計算量は多い。外部騒音伝搬問題などを扱う音場解析は解析領域と周波数範囲が広く、計算機能力の上限から適用範囲はおのずと限定されたものになる。

一方、非特許文献１、２のように、ＦＥＭやＢＥＭと異なり、領域や境界の離散化を必要としない準解析的な解析法による音場解析例が見られる。それらの多くは音場の支配方程式（積分方程式）から導かれるものである。ただし、対象となる音場は円筒形などの軸対称構造物が多く、実用レベルにあるとは言い難い。また、非特許文献１、２は建築音響分野の研究ではなく、水中における線対称または軸対称の物体まわりの散乱波解析を扱ったものであるため、対象とする解析周波数は低く、建築音響分野では用いない境界条件を設定していることなどから、この解析法をそのまま建築音響分野に適用することはできないという問題がある。

本発明は、このような事情に鑑みてなされたもので、準解析手法に基づく任意散乱体まわりの音場解析を行う音場解析装置を提供することを目的とする。

本発明は、音波入射条件情報を入力する音波入射条件入力手段と、散乱体形状情報を入力する散乱体形状入力手段と、境界条件情報を入力する境界条件入力手段と、前記音波入射条件情報に基づいて、直接波のＦｏｕｒｉｅｒ展開を行う直接波Ｆｏｕｒｉｅｒ展開手段と、前記散乱体形状情報に基づいて、散乱波のＦｏｕｒｉｅｒ展開を行う散乱波Ｆｏｕｒｉｅｒ展開手段と、前記境界条件情報、前記直接波Ｆｏｕｒｉｅｒ展開手段が出力する直接波展開情報及び散乱波Ｆｏｕｒｉｅｒ展開手段が出力する散乱波展開情報を入力し、前記散乱体形状情報に基づいて境界条件を展開する境界条件展開手段と、前記境界条件展開手段が出力する展開済みの境界条件情報から連立１次方程式を導出する方程式導出手段と、前記方程式導出手段によって導出された前記連立１次方程式を解くことによって散乱波情報を算出する散乱波算出手段と、前記直接波Ｆｏｕｒｉｅｒ展開手段から出力される直接波展開情報と前記散乱波算出部から出力される散乱波情報を合成して出力する合成手段とを備えたことを特徴とする。

本発明によれば、建築物設計・企画の初期段階における音場の特徴の把握および現象の理解に利用でき、基本設計、騒音対策案の策定や対策箇所の絞込みなどへ活用することができる。またＢＥＭなどの数値解析法の精度検証用ツールとして利用することが可能となる。本発明による解析法は線対称構造物まわりの２次元音場を定式化したものであるが、本解析法と同様の定式化によって、以下の音場解析問題にも対応することができる。
（１）適切なＢｅｓｓｅｌ関数を導入すれば、室内（内部）問題への適用を容易に実現することができる。
（２）一般的な非対称な散乱体まわりの音場を扱うことができる。これは、Ｆｏｕｒｉｅｒ余弦展開ではなく、一般的なＦｏｕｒｉｅｒ展開を用いて実現すればよい。
（３）適切な境界条件を設定すれば、強制振動する構造体からの放射音場を扱うことが容易に実現することができる。
（４）３次元問題への拡張を容易に行うことができる。

以下、本発明の一実施形態による音場解析装置を図面を参照して説明する。まず、本発明による音場解析装置の解析原理を説明する。以下の説明において、「ｒ」、「ｑ」は、ベクトルを表す。ただし、「ｒ」、「ｑ」がスカラー量である場合は、「スカラーｒ」または「スカラーｑ」と表記する。また、式中においては、太字の「ｒ」、「ｑ」がベクトルを表し、細字の「ｒ」、「ｑ」がスカラー量を表す。

＜座標系と記号定義＞
解析対象は、線対称形の散乱体まわりの２次元音場である。座標系はどのようにとっても一般性は失わない。図１に示すように、解析領域Ω内に線対称形の散乱体Ｃがあり、その境界を∂Γとする。散乱体の線対称軸をｘ軸とし、散乱体Ｃの内部にあってｘ軸上に原点を設定した円座標系を設定する。ただし、境界∂Γ上の任意の点ｑ＝（スカラーｑ，φ）は原点から見通せる位置にあるものとする、ここに＋ｘ方向から入射波Φ^〈ｉ〉（ｒ）が到来する。このとき、入射波Φ^〈ｉ〉（ｒ）もまたｘ軸に対して線対称な分布を示すものとする。ここでは、後述する入射条件を設定した場合の解法を述べる。

音場の時間変化がｅｘｐ（−ｉωｔ）で表されるならば、音場はＨｅｌｍｈｏｌｔｚ方程式で記述される。以下では時間項ｅｘｐ（−ｉωｔ）を省略し、本問題をＨｅｌｍｈｏｌｔｚ問題として扱う。速度ポテンシャルのΦ（ｒ）は入射波Φ^〈ｉ〉（ｒ）および散乱波Φ^〈ｓ〉（ｒ）の１次結合として次式で表される。

速度ポテンシャルΦ（ｒ）はＳｏｍｍｅｒｆｅｌｄの放射条件式（５９）、および境界∂Γ上で次式が成り立つものとする。

ここにｎ＝（ｕ，υ）はｑにおける法線ベクトルを表す。またａ_１≠０とする。

＜固有関数展開による散乱波表現＞
入射波と散乱体形状はいずれもｘ軸について線対称であることから、散乱波Φ^〈ｓ〉（ｒ）もまたｘ軸について線対称に分布するのは明らかである。このとき、散乱波Φ^〈ｓ〉（ｒ）はθの偶関数として記述することができる。その周期は２πである。すなわち散乱波Φ^〈ｓ〉（ｒ）はθに関してＦｏｕｒｉｅｒ余弦展開が可能である。Ｓｏｍｍｅｒｆｅｌｄの放射条件（５９）を考慮すれば、散乱波Φ^〈ｓ〉（ｒ）はＨａｎｋｅｌ−Ｆｏｕｒｉｅｒ展開形として式（６９）と同様に次式で表される。

ここに｛ｃ_ｍ；ｍ＝０，±１，±２，．．．｝は未知の展開係数であり、一般に複素数である。本問題は展開係数ｃ_ｍを求める問題に帰着される。展開係数ｃ_ｍは境界条件（１）を満足するように定めてやればよい。

＜境界条件の展開＞
図２に示すように、受音点ｒ＝（スカラーｒ，θ）を境界∂Γ上の点ｑ＝（スカラーｑ，φ）に近づけた極限を考える。このとき、境界条件（１）から次式が成立する。

境界∂Γ上の任意点ｑ（スカラーｑ，φ）は原点から見通せるならば、動径スカラーｑは方位角φの１価関数として記述することが可能である。このとき、スカラーｑはφの偶関数であり、その周期は２πである。したがって、Φ^〈ｓ〉（ｑ）およびΦ^〈ｉ〉（ｑ）もφの偶関数として記述でき、次式のようにφに関するＦｏｕｒｉｅｒ余弦展開形として表現できる。

ここに｛Ψ^〈ｓ〉 _ｎ；ｎ＝０，±１，±２，．．．｝および｛Ψ^〈ｉ〉 _ｎ；ｎ＝０，±１，±２，．．．｝は展開係数である。

同様に考えて、法線ベクトルｎもφの関数として記述できることから、∂Φ^〈ｓ〉（ｑ）／∂ｎおよび∂Φ^〈ｉ〉（ｑ）／∂ｎも次のようにφの関数としてＦｏｕｒｉｅｒ余弦展開形で記述することが可能である。

ここに｛Ψ^〈ｓ〉 _ｎ’；ｎ＝０，±１，±２，．．．｝および｛Ψ^〈ｉ〉 _ｎ’；ｎ＝０，±１，±２，．．．｝は展開係数である。

式（５）、（７）、（９）、（１１）を式（４）に代入すると次式となる。

式（１３）が任意のφに関して成り立つ必要十分条件は、Ｆｏｕｒｉｅｒ級数展開の性質より次式で表される。

入射波は既知であるため、式（１４）の右辺は直接求めることができる。一方、左辺は未知の展開係数ｃ_ｍを含む。すなわち、式（１４）は展開係数ｃ_ｍに関する連立１次方程式である。以下では、式（１４）左辺のΨ^〈ｓ〉 _ｎとΨ^〈ｓ〉 _ｎ’を具体的に求める。いまスカラーｑをφの関数ｑ（φ）と表すとする。このとき、式（２）を式（６）に代入して次式を得る。

ここにΨ^〈ｓ〉 _ｍ，ｎは次式で表される。

一方、∂Φ^〈ｓ〉（ｑ）／∂ｎは便宜的に直交座標系ｘｙで表して次式のようになる。

法線ベクトルｎはφの関数として次のように表される。

ここにスカラーｑ’（φ）はスカラーｑ（θ）のθに関する１階偏微分においてθ＝φとした結果を表す。

式（１７）の∂Φ^〈ｓ〉（ｑ）／∂ｘ，∂Φ^〈ｓ〉（ｑ）／∂ｙは、φ＝ｔａｎ^−１（ｙ／ｘ），スカラーｑ（φ）＝√（ｘ^２＋ｙ^２）を考慮して、それぞれ次式となる。

式（１９）および式（２０）を式（１７）に代入すると次式となる。

ここにＮ_１（φ）、Ｎ_２（φ）は次式で表される。

式（２１）を式（１０）に代入して次式を得る。

ここにΨ^〈ｓ〉 _ｍ，ｎ’は次式で与えられる。

式（１５）、（２４）を式（１４）に代入して次式を得る。

ただし、Υ_ｍ，ｎおよびΘ_ｎは次式で表される。

式（２６）はｃ_ｍを未知数とする無限連立１次方程式である。これを解けば、散乱波の展開係数ｃ_ｍを求めることができる。しかしながら、この無限連立１次方程式は数値解析的に解かざるを得ない。本発明では、級数和を適当な次数で打ち切って有限連立１次方程式に近似する方法を採用する。次数ｍ，ｎの上限をそれぞれＭ，Ｎとするとき、式（２６）は次のようなマトリックス方程式で近似的に表される。

式（２９）を解いて得られる展開係数ｃ_ｍを用いて、式（２）から散乱波Φ^〈ｓ〉（ｒ）は近似的に次式で求められる。

式（３０）を式（６０）に代入すれば、速度ポテンシャルΦ（ｒ）が得られる。

＜入射条件に応じた解法＞
３種の入射条件（円筒波の散乱、二つの円筒波の散乱、平面波の散乱）における散乱音場を定式化する。この入射波はいずれもｘ軸について線対称に分布するものとする。線対称な入射波の多くは、これら３種の入射波を１次結合させて表現することができる。

（１）円筒波の散乱
図３に示すように、散乱体の外部かつｘ軸上に点音源ｒ’＝（スカラーｒ’，０）にあるとする。入射波の速度ポテンシャルΦ^〈ｉ〉（ｒ）は基本解を用いて式（７０）で表されるものとする。式（７０）を散乱波と同じようにＢｅｓｓｅｌ−Ｆｏｕｒｉｅｒ展開形で表せば式（７２）、すなわち次式となる。

受音点ｒを境界∂Γ上の点ｑに近づけた極限を考える。先に説明したように、境界上の点ｑの動径スカラーｑが方位角φの偶関数として記述できることに着目すれば、入射波およびその法線方向微分もまたφの関数として次のように記述できる。

入射波およびその法線方向微分はスカラーｒの関数と、スカラーｒ’の関数の積に分離された形で表現される。このとき、式（８），（１２）からΨ^〈ｉ〉 _ｎおよびΨ^〈ｉ〉 _ｎ’は次式となる。

ここにΨ^〈ｉ〉 _ｍ，ｎおよびΨ^〈ｉ〉 _ｍ，ｎ’はそれぞれ次式で表される。

式（３５）、（３６）から明らかなように、Ψ^〈ｉ〉 _ｍ，ｎおよびΨ^〈ｉ〉 _ｍ，ｎ’は音源の位置や大きさによらない。すなわち、この積分は入射条件に無関係なものとして評価することができる。したがって、連立１次方程式（２９）の右辺Θｎは次のようになる。

式（３７）を式（２９）に代入し、これを解いて得られたら用いて式（３０）から散乱波Φ^〈ｓ〉（ｒ）を求められる。入射波をＢｅｓｓｅｌ−Ｆｏｕｒｉｅｒ展開形で表した場合、Ψ^〈ｓ〉 _ｍ，ｎおよびΨ^〈ｉ〉 _ｍ，ｎ被積分関数に合まれるＢｅｓｓｅｌ関数が異なる点を除けば同じ形になる。またΨ^〈ｓ〉 _ｍ，ｎ’とΨ^〈ｉ〉 _ｍ，ｎ’についてもＢｅｓｓｅｌ関数以外は同じ形となる。Ｈａｎｋｅｌ関数Ｈ^（１） _ｍ（ｚ）が第１種Ｂｅｓｓｅｌ関数Ｊ_ｍ（ｚ）と第２種Ｂｅｓｓｅｌ関数Ｙ_ｍ（ｚ）の１次結合、すなわち

で表されることに着目すれば、次の関係が成り立つのは明らかである。

ここにＲは複素数の実部を表す。

したがって、Ψ^〈ｓ〉 _ｍ，ｎ、Ψ^〈ｓ〉 _ｍ，ｎ’を求めると、直ちにΨ^〈ｉ〉 _ｍ，ｎ、Ψ^〈ｉ〉 _ｍ，ｎ’が得られることになる。実際の音場解析において、音源の位置や数などの入射条件を変化させて解析し、散乱音場の特徴を考察する場合は多い。後述するが、二つの円筒波が入射する場合および平面波入射する場合も積分Ψ^〈ｓ〉 _ｍ，ｎ，Ψ^〈ｓ〉 _ｍ，ｎ’，Ψ^〈ｉ〉 _ｍ，ｎ，Ψ^〈ｉ〉 _ｍ，ｎ’はそれぞれ式（１６）、（２５）、（３５）、（３６）と全く同じ形で求められる。すなわち散乱体の形状が決まれば、これらの積分は評価できる。言い換えれば、散乱体の音響的な性質（大きさ、形状、境界の吸音条件など）が変化しない限り積分を再評価する必要はない。そのため、数多くの入射条件を扱う場合などでは、解析全体にかかるコストの低減が期待できる。

（２）二つの円筒波の散乱
図４に示すように、散乱体外部にあってｘ軸について線対称関係にある二つの点音源ｒ’＝（スカラーｒ’，θ’）、ｒ''＝（スカラーｒ’，−θ’）があるとする。このとき、入射波Φ^〈ｉ〉（ｒ）は基本解（７０）を用いて次式で表されるものとする。

なおθ’＝０の場合、式（４１）は式（７０）に一致する。式（４１）をＢｅｓｓｅｌ−Ｆｏｕｒｉｅｒ展開形で表すと次式となる。

式（４２）と式（７２）の違いはｃｏｓｍθ’の有無であり、これはφに無関係である。したがって式（４２）から式（２９）の右辺Θ_ｎを求める手順は前述した手順と全く同じであり、次式で与えられる。

当然ながら、Ψ^〈ｉ〉 _ｍ，ｎおよびΨ^〈ｉ〉 _ｍ，ｎ’はそれぞれ式（３５）、（３６）で与えられる。式（４３）を式（２９）に代入し、これを解いて得られたｃ_ｍを用いて、式（３０）から散乱波Φ^〈ｓ〉（ｒ）を求めることができる。この解析条件では、図５に示すようにｘ軸上に完全反射境界があるとし、音源ｒ''を完全反射境界に対する音源ｒ’の虚像と見なすことができる。すなわち、本解析条件は、線形散乱体Ｃの半分を含む半無限空間における音場を表すともいえる。

（３）平面波の散乱
図６に示すように、平面波Φ^〈ｗ〉（ｒ）が入射角πで到来するとする。このとき、入射波の速度ポテンシャルΦ^〈ｉ〉（ｒ）は式（８２）、すなわち次式で表されるものとする。

式（８２）をＢｅｉｓｓｅｌ−Ｆｏｕｒｉｅｒ展開形で表すと式（８３）、すなわち次式となる。

式（８３）から式（２９）の右辺Θ_ｎを求める手順は前述の手順と全く同じであり、次式で与えられる。

ここにΨ^〈ｉ〉 _ｍ，ｎおよびΨ^〈ｉ〉 _ｍ，ｎ’はそれぞれ式（３５）、（３６）で与えられる。式（４４）を式（２９）に代入し、これを解いて得られたｃ_ｍを用いて式（３０）から散乱波Φ^〈ｓ〉（ｒ）を求めることができる。

＜解析解との対応＞
次に、円筒波が入射する条件（図３）を例にとり、散乱体Ｃの形状が円形であるときに境界展開法による解は解析解（８１）に一致することを示す。散乱体Ｃが半径ｑ_０の円ならば、動径ｑ（φ）およびその微分は次式となる。

このとき、式（１６）、（２５）、（３５）、（３６）からΨ^〈ｓ〉 _ｍ，ｎ，Ψ^〈ｓ〉 _ｍ，ｎ’，Ψ^〈ｉ〉 _ｍ，ｎ，Ψ^〈ｉ〉 _ｍ，ｎ’はそれぞれ次式となる。

三角関数の直交性から次式が成り立つ。

式（４６）〜式（５０）を式（２７）、式（３７）に代入して、式（２６）を整理すると次のようになる。

式（２６）と式（５１）を比較すれば明らかなように、境界形状が円の場合は級数和（Σ^∞ _ｍ＝０）は消失する。すなわち連立１次方程式（２９）を解く必要がなくなり、展開係数ｃ_ｍは陽に次式で求められる。

式（５２）を式（２）に代入して散乱波Φ^〈ｓ〉（スカラーｒ，θ）を得る。最終的に速度ポテンシャルΦ（ｒ）は次式で表される。

式（５３）は解析解（８１）に一致する。平面波が入射する場合も同様にして、本解析法による解は解析解（９０）に一致することが確認される。以上から、境界展開法は解析解を包含しており、適用範囲の広い解析法であるといえる。

＜音場の支配方程式＞
音場の時問変化がｅｘｐ（−ｉωｔ）で表されるとき、受音点ｒにおける速度ポテンシャルΦ（ｒ）は次の非同次Ｈｅｌｍｈｏｌｔｚ方程式の解として求められる。

ただしρ（ｒ）は音源の寄与を表す既知関数とする。

＜基本解＞
Ｇｒｅｅｎ関数（Green's function）は振幅の大きさが１の点音源からの寄与を表す。受音点をｒ、音源をｒ’とするとき、Ｈｅｌｍｈｏｌｔｚ方程式に対するＧｒｅｅｎ関数Ｇ（ｒ；ｒ’）は、次の非同次方程式（inhomogeneous equation）の解として与えられる。

ここにδはＤｉｒａｃのデルタ関数を表す。

自由空間中、すなわち境界が無限遠にある場合のＧｒｅｅｎ関数は基本解（fundamental solution）に等しい。３次元の基本解Ｇ_３ ^±（ｒ；ｒ’）はスカラーｒ＝｜ｒ−ｒ’｜として次式で表される。

２次元の基本解Ｇ_２ ^±（ｒ；ｒ’）は次式で表される。

ここにＨ_０ ^（１）は第１種０次Ｈａｎｋｅｌ関数を表す。

＜放射条件＞
領域Ωが開領域である場合、無限遠で速度ポテンシャルは有界であり、かつ無限遠からは反射波が帰ってこないと見なせる。これを満たす境界条件はＳｏｍｍｅｒｆｅｌｄの放射条件と呼ばれる。３次元問題におけるＳｏｍｍｅｒｆｅｌｄの放射条件は、スカラーｒ＝｜ｒ｜として次式で表される。

一方、２次元問題では次式となる。

＜固有関数展開法による音場解析＞
固有関数展開法（eigenfunction expansion method）は、速度ポテンシャルを基礎方程式の固有関数で展開し、境界条件を満足するように展開係数を決定する方法である。ここでは、２次元外部問題として円形散乱体による音波の散乱を扱う。これは解析解（厳密解）が得られる問題である。

（１）円筒波入射の場合
図７に示すように、解析領域Ω内に半径ｑ_０の円形散乱体Ｃ（境界：∂Γ）があるとする。座標系はどのようにとっても一般性は失わない。ここでは散乱体Ｃの中心を原点に設定した円座標系を用いる。受音点ｒ＝（スカラーｒ，θ）および音源ｒ’＝（スカラーｒ’，０）はいずれも散乱体Ｃの外部にあるとする。また散乱体Ｃ境界∂Γ上の任意の点をｑ（ｑ０，φ）とする。速度ポテンシャルの時間変化をｅｘｐ（−ｉωｔ）とすると、受音点ｒにおける速度ポテンシャルの空間依存項Φ（ｒ）は、Ｈｅｌｍｈｏｌｔｚ方程式（５４）の解として与えられる。また音場は線形であることから、Φ（ｒ）は音源からの直接波Φ^〈ｉ〉（ｒ）と散乱体による散乱波Φ^〈ｓ〉（ｒ）の１次結合として次式で表される。

当然ながら、Φ^〈ｉ〉（ｒ）およびΦ^〈ｓ〉（ｒ）もまたＨｅｌｍｈｏｌｔｚ方程式を満足する。速度ポテンシャルΦ（ｒ）はＳｏｍｍｅｒｆｅｌｄの放射条件（５９）を満足し、散乱体Ｃの境界∂Γ上では完全反射条件、すなわち次式を満足する。

以下では、まず散乱波Φ^〈ｓ〉（ｒ）に着目して説明する。前述したようにΦ^〈ｓ〉（ｒ）は、次の同次Ｈｅｌｍｈｏｌｔｚ方程式の解である。

ここに∇^２はＬａｐｌａｃｉａｎであり、２次元円座標系では∇^２＝∂^２／∂ｒ^２＋（１／ｒ）∂／∂ｒ＋∂^２／∂θ^２である。ここで、ｒ、θはスカラー量である。

散乱音場は、音源と散乱体の中心を結ぶ直線に対して対称になるのは明らかである。すなわち、Φ^〈ｓ〉（ｒ）はθの偶関数として記述できよその周期は２πである。このときΦ^〈ｓ〉（ｒ）は次に示すようにＦｏｕｒｉｅｒ余弦展開が可能である。

ここに｛Ψ_ｍ（ｒ）；ｍ＝０，１，２．．．｝はＦｏｕｒｉｅｒ展開係数である。またε_ｍ，ｎは次式で与えられるＮｅｕｍａｎｎ因子である。

式（６３）を式（６２）に代入して整理して次式を得る。

式（６６）が任意のθに関して成り立つ必要十分条件は次式で与えられる。

式（６７）は、ｍ次Ｂｅｓｓｅｌ方程式である。式（６７）の解のうち、放射条件（５９）を満たすものは第１種ｍ次Ｈａｎｋｅｌ関数である。
したがって、式（６７）の解Ψ_ｍ（ｒ）は適当な係数｛Ｃ_ｍ；ｍ＝０，１，２，．．．｝を用いて次のように表されるものとする。

ここにＨ_ｍ ^（１）はｍ次１種Ｈａｎｋｅｌ関数である。

式（６８）を式（６３）に代入して次式を得る。

一方、直接波Φ^〈ｉ〉（ｒ）は点音源からの放射波の寄与を表す。このときΦ^〈ｉ〉（ｒ）は基本解（５７）を用いて次式で表されるものとする。

式（７０）に含まれるＨａｎｋｅｌ関数は、Ｂｅｓｓｅｌ関数の加法定理を用いて次式で表される。

ここに｛スカラーｒ＞スカラーｒ’；スカラーｒ⇔スカラーｒ’｝はこの式がスカラーｒ＞スカラーｒ’に対して成り立ち、スカラーｒ’＞スカラーｒに対してはｒとｒ’を入れかえた式が成り立つことを表すものとする。

式（７１）を式（７０）に代入すれば、Φ^〈ｉ〉（ｒ）をスカラーｒ、スカラーｒ’およびθを変数分離した形で次式のように表すことができる。

式（７２）はその形からθに関するＦｏｕｒｉｅｒ余弦展開形と考えることができる。速度ポテンシャルΦ（ｒ）を求めるには、式（７２）、（６９）から境界条件（６１）を満足するように展開係数ｃ_ｍを決定すればよい。いまスカラーｒ’＞スカラーｒとすると、直接波Φ^〈ｉ〉（ｒ）および散乱波Φ^〈ｓ〉（ｒ）のスカラーｒに関する１階偏微分はそれぞれ次式のようになる。

ここにＪ'_ｍ（ｋスカラーｒ）およびＨ_ｍ ^（１）’（ｋスカラーｒ）はそれぞれＪ_ｍ（ｚ）およびＨ_ｍ ^（１）（ｚ）のｚに関する１階偏微分においてｚ＝ｋスカラーｒとした結果を表す。

式（７３）および式（７４）を式（６１）に代入して整理すると次のようになる。

式（７５）が任意のφに関して成り立つ必要十分条件は次式である。

したがって展開係数ｃ_ｎは次式で与えられる。

式（７７）を式（６９）に代入して次式を得る。

式（７２）および式（７８）を式（６０）に代入して次式を得る。

なお境界∂Γ上で次に示す同次境界条件

が成り立つならば、Φ（ｒ）は次のようになる。

（２）平面波入射の場合
図７における点音源ｒ’に代えて、次式に示す平面波Φ^〈ｗ〉（ｒ）が入射角πでｘ負方向へ進行するとする。

式（８２）は、指数関数部を展開して次式のように表すことができる。

また∂Φ^〈ｗ〉（ｒ）／∂スカラーｒは次のようになる。

当然ながら、境界条件（６１）においてΦ^〈ｉ〉をΦ^〈ｗ〉に置き換えた式は成立し、次のようになる。

式（８５）が任意のθに関して成り立つ必要十分条件は次式である。

したがって展開係数ｃｍは次式で与えられる。

式（８７）を式（６９）に代入して次式を得る。

式（８３）および式（８８）を式（６０）に代入して次式を得る。

境界∂Γ上で式（８０）が成り立つならば、Φ（ｒ）は次のようになる。

次に、図８、９を参照して、本実施形態の構成を説明する。図１は同実施形態の構成を示すブロック図である。この図において、符号１は、前述した解析原理を用いて音場解析処理を実行するコンピュータである。符号１１は、音波入射条件情報を入力して、直接波のＦｏｕｒｉｅｒ展開を実行する直接波Ｆｏｕｒｉｅｒ展開部である。符号１２は、散乱体（建造物）の形状情報を入力して、散乱波のＦｏｕｒｉｅｒ展開を実行する散乱波Ｆｏｕｒｉｅｒ展開部である。符号１３は、境界条件情報、直接波Ｆｏｕｒｉｅｒ展開部１１が出力する直接波情報及び散乱波Ｆｏｕｒｉｅｒ展開部１２が出力する散乱波情報を入力して、散乱体形状に基づいて境界条件を展開する境界条件展開部である。符号１４は、境界条件展開部１３が出力する展開済みの境界条件情報から連立１次方程式を導出する方程式導出部である。符号１５は、方程式導出部１４によって導出された連立１次方程式を解くことによって散乱波情報を算出する散乱波算出部である。符号１６は、直接波Ｆｏｕｒｉｅｒ展開部１１から出力される直接波情報と散乱波算出部１５から出力される散乱波情報を合成して出力する合成部である。

次に、図９を参照して、現実の解析条件から特徴を抽出し、解析モデルを再構成する方法について説明する。たとえば、図９（ａ）のような交通量の多い幹線道路を騒音源として、対象建物２への騒音伝搬を解析する場合、騒音対策検討の初期段階において、詳細なモデル構築はコストがかかる。そのため、図９（ｂ）のように解析モデルを簡略化することは一般に行われる。具体的には、建物形状２を左右対称形とした対象建物３と見なし、音源を線音源などの単純な条件に置き換えるなどである。このとき、音場の線対称性が成り立つと見なせれば、境界展開法が適用でき、ＢＥＭよりも効率的な解析が可能になる。この方法によって再構成された解析モデルがコンピュータ１への入力情報となる。

次に、図８を参照して、図８に示すコンピュータ１の解析動作を説明する。まず、解析者は、現実の解析条件から特徴を抽出し、解析モデルを再構成する。これによって、図９（ｂ）に示す解析モデルが構成されることになる。そして、解析者は、音波入射条件情報、散乱体形状情報及び境界条件情報を入力する。これを受けてコンピュータ１は、これらの入力情報を読み取る。そして、直接波Ｆｏｕｒｉｅｒ展開部１１は、音波入射条件を入力して、直接波をＦｏｕｒｉｅｒ展開する。また、散乱波Ｆｏｕｒｉｅｒ展開部１２は、散乱体形状情報を入力し、散乱波をＦｏｕｒｉｅｒ展開する。

次に、境界条件展開部１３は、境界条件情報、直接波Ｆｏｕｒｉｅｒ展開部１１が出力する直接波情報及び散乱波Ｆｏｕｒｉｅｒ展開部１２が出力する散乱波情報を入力して、散乱体形状に基づいて境界条件を展開する。そして、方程式導出部１４は、境界条件展開部１３が出力する展開済みの境界条件情報から連立１次方程式を導出し、散乱波算出部１５は、この連立１次方程式を解くことによって散乱波情報を算出する。続いて、合成部１６は、直接波Ｆｏｕｒｉｅｒ展開部１１から出力される直接波情報と散乱波算出部１５から出力される散乱波情報を合成して出力する。解析者は、合成部１６が出力した結果を参照して、音場の評価を行う。

次に、級数和の打ち切り次数の決定について説明する。前述したように、散乱波Φ^〈ｓ〉（ｒ）の展開係数｛ｃ_ｍ；ｍ＝０，１，２，．．．，Ｍ｝を求めるためには、連立１次方程式の具体的な形を定める必要がある。この連立１次方程式の元数（未知数の個数）は、散乱波に関する級数和の打ち切り次数Ｍに等しい。有意な解をもつためには二つの級数和の打ち切り次数Ｍ、Ｎに関してＭ≦Ｎとする必要がある。仮にＭ＜Ｎとした場合、連立１次方程式は過剰な連立１次方程式となるので、最小二乗法により最適解を求めることになる。最適なＭとＮの組み合わせは、散乱体の形状や周波数などにより決まると考えられるが、その値は現状では明らかではない。ここでは単純に考えてＭ＝Ｎとして、連立１次方程式の元数の決定法について述べる。このとき式（２９）左辺の行列はＭ×Ｍの正方行列となる。

いま、式（２９）右辺Θｎに含まれる積分項Ψ^〈ｉ〉 _ｍ，ｎおよびΨ^〈ｉ〉 _ｍ，ｎ’に着目する。前述したようにΨ^〈ｉ〉 _ｍ，ｎおよびΨ^〈ｉ〉 _ｍ，ｎ’は入射波によらず、散乱体の形状が決まれば求めることができる。またＦｏｕｒｉｅｒ級数展開の性質からΨ^〈ｉ〉 _ｍ，ｎおよびΨ^〈ｉ〉 _ｍ，ｎ’はｍ，ｎが大きくなるにつれて０に近づくことが予想される。図１０に示すように、解析領域Ω内に矩形の散乱体Ｃ（寸法が２Ｗ×２Ｈ）を考える。以下では、基準長をＬとしてＷ＝Ｈ＝ＬであるときのΨ^〈ｉ〉 _ｍ，ｎおよびΨ^〈ｉ〉 _ｍ，ｎ’の収束性について検討する。なお散乱体Ｃの境界上では完全反射条件を満足するものとする。

矩形散乱体Ｃに対して、次数ｍ、ｎをそれぞれ０〜３０次まで変化させた場合の｜Ψ^〈ｉ〉 _ｍ，ｎ｜および｜Ψ^〈ｉ〉 _ｍ，ｎ’｜を図１１に示す。上段からそれぞれｋＬ＝２．５、５．０、１０．０の場合である。各段の左側は｜Ψ^〈ｉ〉 _ｍ，ｎ｜、右側は｜Ψ^〈ｉ〉 _ｍ，ｎ’｜を表す。いずれも解析結果の最大値で正規化した値である。図１１から、｜Ψ^〈ｉ〉 _ｍ，ｎ｜および｜Ψ〈^ｉ〉 _ｍ，ｎ’｜はｍ，ｎが小さくかつｍ〜ｎのときに大きな値をもつことがわかる。この傾向はｋＬによらない。またｍ、ｎのいずれか一方が大きくなるにつれて、いずれの積分項も０に近づく。いいかえれば、ｍ，ｎのいずれか一方が十分大きければ、Ψ^〈ｉ〉 _ｍ，ｎ、Ψ^〈ｉ〉 _ｍ，ｎ’が音場におよぼす影響は小さいと考えてよい。以上のことからｍ＝ｎ＝ｍ_０の場合を考え、十分小さな正数をε_１、ε_２として

の関係が成立するような最大のｍ_０を用いれば、散乱波を表現するのに十分な次数が確保されると見なす。このときのｍ_０を打ち切り次数、すなわち連立１次方程式の元数Ｍとする。ここではε_１＝ε_２＝１．０×１０^−１２とする。

次に、本発明による解析結果と、従来の境界要素法（ＢＥＭ）による解析結果との比較結果について説明する。図１２に、解析に用いた６種（円形、矩形Ａ、矩形Ｂ、八角形、楕円Ａ、楕円Ｂ）の散乱体形状を示す。図１３は、図１２に示す散乱体形状のそれぞれが、音場に及ぼす影響を解析した結果を示す図である。図１３において、実線は、本発明による解析結果を示し、破線は、従来の境界要素法を用いて解析した結果を示している。また、図１４は、本発明による解析解の安定性を示す図である。図１４において、実線は、本発明による解析結果を示し、破線は、従来の境界要素法を用いて解析した結果を示している。また、図１４の右側（（ａ）、（ｃ）、（ｅ））は、散乱波の相対音圧レベルを示す、左側（（ｂ）、（ｄ）、（ｆ））は散乱波音圧レベルの差を示している。図１３、図１４から明らかなように、本発明による解析解と従来の境界要素法による解析解とは良く一致していることが分かる。

このように、本発明による解析方法は、比較的単純な形状の散乱体まわりの音場を対象とした準解析的な波動場解析法である。これはＦＥＭ（有限要素法）、ＢＥＭ（境界要素法）などの領域や境界の離散化を伴う数値解析法とは異なり、固有関数展開法に基づくものである。具体的に言えば、波動場を固有関数で展開し、境界条件式をさらにＦｏｕｒｉｅｒ級数展関することによって音場の展開係数に関する連立１次方程式を解く問題に帰着させる方法である。この解析法は以下に示す特徴を有している。
（１）領域や境界を離散化する必要がない準解析解法であり、解析解（厳密解）を包含する。
（２）連立１次方程式の大きさ（元数）は同一解析条件における境界要素法に比べて小さく（約１／１０）、計算機の記憶容量の点からいえば有利である。
（３）矩形など角部を有する散乱体まわりの音場を解析する際、角部をわずかに丸めることによって実用上十分な精度を確保しつつ計算時間を短縮することが可能である。
（４）散乱体形状が円形、あるいは縦横比が１に近い場合は十分な精度を確保できる。
（５）境界要素法における見かけの固有周波数問題のように、特定の周波数における精度悪化は生じない。

また、本発明による解析法と同様の定式化によって、以下の音場も扱うことが可能である。
（ａ）線対称音場ばかりでなく、任意形状の散乱体まわりの非対称な音場。このとき境界展開には一般的なＦｏｕｒｉｅｒ級数展開を用いることになる。
（ｂ）線対称な構造物がその圏外方向に一定周期で振動する場合の放射音場。
（ｃ）線対称な室内の音場を扱うような内部問題。このとき、解析領域内で有界となる適当なＢｅｓｓｅｌ関数を用いることによって定式化される。

なお、図１における処理部の機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより音場解析処理を行ってもよい。なお、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムが送信された場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリ（ＲＡＭ）のように、一定時間プログラムを保持しているものも含むものとする。

また、上記プログラムは、このプログラムを記憶装置等に格納したコンピュータシステムから、伝送媒体を介して、あるいは、伝送媒体中の伝送波により他のコンピュータシステムに伝送されてもよい。ここで、プログラムを伝送する「伝送媒体」は、インターネット等のネットワーク（通信網）や電話回線等の通信回線（通信線）のように情報を伝送する機能を有する媒体のことをいう。また、上記プログラムは、前述した機能の一部を実現するためのものであってもよい。さらに、前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるもの、いわゆる差分ファイル（差分プログラム）であってもよい。

線対称散乱体まわりの音場解析のための記号定義を示す説明図である。 境界上の点と法線ベクトルの関係を示す説明図である。 線対称散乱体と円筒波の関係を示す説明図である。 線対称散乱体と二つの円筒波の関係を示す説明図である。 完全反射境界を有する半無限領域を示す説明図である。 線対称散乱体と平面波の関係を示す説明図である。 固有関数展開法による音場解析のための記号定義を示す説明図である。 本発明の一実施形態の構成を示すブロック図である。 図８に示す解析モデルの再構成例を示す説明図である。 矩形散乱体を示す説明図である。 境界展開に基づく積分項の収束性を示す説明図である。 計算モデルの散乱体形状を示す説明図である。 解析結果の一例を示す図である。 解析結果の一例を示す図である。

符号の説明

１・・・コンピュータ、１１・・・直接波Ｆｏｕｒｉｅｒ展開部、１２・・・散乱波Ｆｏｕｒｉｅｒ展開部、１３・・・境界条件展開部、１４・・・方程式導出部、１５・・・散乱波算出部、１６・・・合成部

Claims (1)

1. 音波入射条件情報を入力する音波入射条件入力手段と、
   散乱体形状情報を入力する散乱体形状入力手段と、
   境界条件情報を入力する境界条件入力手段と、
   前記音波入射条件情報に基づいて、直接波のＦｏｕｒｉｅｒ展開を行う直接波Ｆｏｕｒｉｅｒ展開手段と、
   前記散乱体形状情報に基づいて、散乱波のＦｏｕｒｉｅｒ展開を行う散乱波Ｆｏｕｒｉｅｒ展開手段と、
   前記境界条件情報、前記直接波Ｆｏｕｒｉｅｒ展開手段が出力する直接波展開情報及び散乱波Ｆｏｕｒｉｅｒ展開手段が出力する散乱波展開情報を入力し、前記散乱体形状情報に基づいて境界条件を展開する境界条件展開手段と、
   前記境界条件展開手段が出力する展開済みの境界条件情報から連立１次方程式を導出する方程式導出手段と、
   前記方程式導出手段によって導出された前記連立１次方程式を解くことによって散乱波情報を算出する散乱波算出手段と、
   前記直接波Ｆｏｕｒｉｅｒ展開手段から出力される直接波展開情報と前記散乱波算出部から出力される散乱波情報を合成して出力する合成手段と
   を備えたことを特徴とする音場解析装置。

Images