(57)【要約】【課題】本発明は、境界要素法による波動音響解析方法を適用して、広域な３次元空間の音場解析を実現できる３次元音響波動解析方法を提供する。【解決手段】３次元音場の解析領域Ωを３次元空間に設定し、速度ポテンシャルΦを解とする支配方程式、及び該支配方程式に対する境界条件を算定する。この後、支配方程式の地盤面に関する境界条件を満足する境界条件適合型応答関数を算定し、該境界条件適合型応答関数を導入して支配方程式をこれと等価な境界積分方程式に変換する。次に、解析領域Ω内の境界面を、対象とする周波数に応じた大きさの面要素に分割して、離散化手法を用いて前記境界積分方程式を解くことにより境界上の速度ポテンシャル分布を算定する。最後に、境界積分方程式により算定された境界上の速度ポテンシャル分布を、前述した境界積分方程式に代入することにより、解析領域Ωの速度ポテンシャルΦを算定する。【選択図】  図２
2004333224.tif 000002
【特許請求の範囲】
【請求項１】
境界要素法による波動音響解析法を介して、少なくとも地盤面及び建物群の外壁面に境界面を有する広域な３次元空間を対象領域とした３次元音場を、速度ポテンシャルを用いて定式化する３次元音響波動解析方法であって、
速度ポテンシャルを解とする支配方程式に、波動場の波源に対する応答を表す関数であり、前記支配方程式の地盤面に関するインピーダンス境界条件と同一の条件を満足する関数を導入して、
前記支配方程式を、前記対象領域内で一様な境界積分方程式に変換することを特徴とする３次元音響波動解析方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、広域な３次元空間を対象とした３次元音響波動解析方法に関する。
【０００２】
【従来の技術】
従来より、音場解析方法には、幾何学音響解析法、及び波動音響解析法が用いられている。幾何学音響解析法は、音の波動性を無視してエネルギーの伝搬工程を幾何学的に解析する手法であり、波動性を考慮しないことから、回折や多重反射を近似的に取り扱う。このため、特に波動性が問題となる低周波数域では、音場の予測精度に課題が生じる。
また、波動音響解析方法は、波動性を考慮するが、解析解を得ることに困難な場合が多く、差分法、有限要素法、境界要素法等の数値解法に基づいて波動現象を解析する。これら差分法や有限要素法は、解析領域全体を格子や要素に分割するために計算量が膨大となり、広域の音場解析への適用は困難である。一方で、境界要素法は、計算量が少ないといった利点を有する。
【０００３】
【非特許文献１】
日本建築学会著、「実務的騒音対策指針（第二版）」、技報堂出版株式会社、ｐｐ１８−２１
【０００４】
しかし、吸音性地盤面を考慮した一般的な条件に対する解析は、計算に時間を要するため実用的ではない。このため、現状では、建物内部を対象とした解析や２次元問題として扱うにとどまる場合が多い。
【０００５】
【発明が解決しようとする課題】
上記事情に鑑み、本発明は、境界要素法による波動音響解析方法を適用して、広域な３次元空間の音場解析を実現できる３次元音響波動解析方法を提供することを目的としている。
【０００６】
【課題を解決するための手段】
請求項１記載の３次元音響波動解析方法は、境界要素法による波動音響解析法を介して、少なくとも地盤面及び建物群の外壁面に境界面を有する広域な３次元空間を対象領域とした３次元音場を、速度ポテンシャルを用いて定式化する３次元音響波動解析方法であって、速度ポテンシャルを解とする支配方程式に、波動場の波源に対する応答を表す関数であり、前記支配方程式の地盤面に関するインピーダンス境界条件と同一の条件を満足する関数を導入して、前記支配方程式を、前記対象領域内で一様な境界積分方程式に変換することを特徴としている。
【０００７】
【発明の実施の形態】
以下、本発明に係る３次元音響波動解析方法について、図１から図６を参照して詳細を以下に示す。
なお、本発明の３次元音響波動解析方法で用いる境界要素法による波動音響解析法は、境界要素法の支配方程式に波動方程式あるいはそのＦｏｕｒｉｅｒ変換によって得られるｈｅｌｍｈｏｌｔｚ方程式を適用して、波動音場を定式化するものとして一般に知られており、音場解析の分野で一般的に用いられている速度ポテンシャルを支配方程式の解とし、該速度ポテンシャルを波動音場の定式化に用いている。
【０００８】
このような境界要素法による波動音響解析法を適用し、広域な３次元空間の三次元音場を定式化する３次元音響波動解析方法を、図１に示すフロー図を用いて簡略に述べる。
まず、３次元音場の解析領域Ωを３次元空間に設定し、地盤面や建物群の外壁面に形成される境界面г、観測点ｐ、点音源ｑｎを設定する（ステップｓ１）。
次に、解析領域Ωにおいて、観測点ｐにおける速度ポテンシャルΦを解とする非同次３次元ｈｅｌｍｈｏｌｔｚ方程式を境界要素法の支配方程式に適用し、該支配方程式に対する境界条件を算定する（ステップｓ２）。
この後、支配方程式の地盤面に関する境界条件と同じ条件を満足する関数（以降、境界条件適合型応答関数と称す）を算定し、該境界条件適合型応答関数を導入して支配方程式をこれと等価な境界積分方程式に変換する。（ステップｓ３）。
次に、解析領域Ω内の境界面を、対象とする周波数に応じた大きさの面要素に分割して、離散化手法を用いて前記境界積分方程式を解くことにより境界上の速度ポテンシャル分布を算定する（ステップｓ４）。
最後に、境界積分方程式により算定された境界面г上の速度ポテンシャル分布を、前述した境界積分方程式に代入することにより、解析領域Ωの速度ポテンシャルΦを算定する（ステップｓ５）。
【０００９】
以下に、上述する行程よりなる３次元音響波動解析方法を、ステップｓ１からステップｓ５の流れに沿って詳述する。
【００１０】
（ステップｓ１）
図２に、３次元音場を定式化する解析領域Ωを示す。本実施の形態では、地盤面１上に広域な３次元空間を設定しており、該３次元空間には、鉛直上向き方向をｚ正方向、水平方向で直交する２方向をｘ、ｙ方向とする直交座標系ｘｙｚが設定されている。ここで、前記地盤面１は、ｚ＝０に設定されており、解析領域Ωは、ｚ≧０の半無限空間としている。また、該解析領域Ωには、複数の点音源｛ｑｎ（ξｎ、ηｎ、ζｎ）｜ｎ＝１、２、・・・、Ｎ｝、観測点ｐ（ｘ、ｙ、ｚ）、任意個数の建物群Ｃが存在しており、該建物群Ｃ各々の外壁面は、境界面гｉ｛гｉ｜ｉ＝１、２、・・、ｍ、・｝として表すこととする。なお、本実施の形態では、該建物群Ｃ各々の境界面гｉをまとめて境界面гと表すこととし、建物群Ｃの境界面гは、（１）式で表される。
また、各境界面гに対して、外向き法線ベクトルをｎとし、地盤面１と境界面гの音響アドミッタンス比をそれぞれＡｇ、Ａと設定している。
【００１１】
【数１】
2004333224.tif 000003
【００１２】
なお、本実施の形態では、点音源ｑｎ、観測点ｐに、一般に音源解析に適用されているスピーカ及びマイクロホンが備えられており、これらはパワーアンプ３及びマイクアンプ４を介して計測装置２に連動されている。３次元音響波動解析方法は、点音源ｑｎ、観測点ｐに係る入出力データ、及び前記計測装置２のプログラム領域に備えられており以下に詳述する３次元音響波動解析方法に係る演算式を用いて、前記計測装置２の演算領域で演算することにより、先に設定した解析領域Ωの３次元音場を定式化するものである。なお、スピーカやマイクロホン等用いる装置は上述するものに必ずしもこだわるものではない。
【００１３】
（ステップｓ２）
解析領域Ωの３次元音場は、速度ポテンシャルΦを用いて定式化するものであり、該速度ポテンシャルΦは、（２）式で表すことができる。本発明では、速度ポテンシャルΦの依存項であるφ（ｐ）を算定すればよいため、以下では時間依存項ｅｘｐ（−ｉωｔ）を省略する。
【００１４】
【数２】
2004333224.tif 000004
【００１５】
観測点ｐにおける速度ポテンシャルφ（ｐ）は、（３）式に示すように、境界要素法の支配方程式として適用した非同次３次元ｈｅｌｍｈｏｌｔｚ方程式の解として表される。
【００１６】
【数３】
2004333224.tif 000005
【００１７】
これら（３）式に示す支配方程式に対する境界条件は、Ｓｏｍｍｅｒｆｅｌｄの放射条件を（４）式、地盤面１に関するインピーダンス境界条件を（５）式、及び境界面гに関するインピーダンス境界条件を（６）式に表すことができる。なお、音響アドミッタンス比が０である境界面は、完全反射面となる。
【００１８】
【数４】
2004333224.tif 000006
【数５】
2004333224.tif 000007
【数６】
2004333224.tif 000008
【００１９】
ところで、（３）式に示す支配方程式は、Ｇｒｅｅｎの定理とＧｒｅｅｎ関数Ｇを用いて（７）式に示す積分方程式に帰着できる。
【００２０】
【数７】
2004333224.tif 000009
【００２１】
このとき、ある音源ｑと観測点ｐに関して、Ｇｒｅｅｎ関数Ｇは（８）式に示す方程式および（９）式に示す放射条件を満足する関数として与えられる。なお、個々の音源ｑｎを区分する必要のない場合には、音源をｑとして表す。
【００２２】
【数８】
2004333224.tif 000010
【数９】
2004333224.tif 000011
【００２３】
（ステップＳ３）
（７）式で用いたＧｒｅｅｎ関数Ｇは、（８）式及び（９）式の他に、速度ポテンシャルφ（ｐ）と同じ境界条件を満足すれば、（７）式の右辺の第２項における境界面гに関する積分を０とし、また、速度ポテンシャルφ（ｐ）と同じ境界条件を満たさない場合には、その境界面г上の積分を数値解析することとなる。
このとき、Ｇｒｅｅｎ関数Ｇが、速度ポテンシャルφ（ｐ）と地盤面１に関する同じ境界条件を満足する場合には、前記地盤面１に関するインピーダンス境界条件（５）式より、Ｇｒｅｅｎ関数Ｇは（１０）式を満足する。
【００２４】
【数１０】
2004333224.tif 000012
【００２５】
従来では、（８）（９）（１０）式を満足するＧｒｅｅｎ関数Ｇを算定していたが、本発明では、Ｇｒｅｅｎ関数Ｇに代わり、同じ条件を満足する境界条件適合型応答関数を算定する。なお、ここで用いる境界条件適合型応答関数は、波動場の波源に対する応答を表す関数である。
（１０）式を満足する境界条件適合型応答関数は、基本解ｅｘｐ（ｉｋｒ）／４πｒと、図３に示すｒ０の関数β（ｒ０）を用いて（１１）式で表されるものと仮定する。
【００２６】
【数１１】
2004333224.tif 000013
【００２７】
【数１２】
2004333224.tif 000014
【数１３】
2004333224.tif 000015
【００２８】
観測点ｐ  が地盤面に近い場合、ｒ０  は、（１４）式と見なせる。このとき、形式的に、β（ｒ０）は、観測点ｐ  のｚ  座標には依存しないものの、すなわち∂β（ｒ０）／∂ｚを０と見なして、（１０）式及び（１１）式から（１５）式が算定される。
【００２９】
【数１４】
2004333224.tif 000016
【００３０】
【数１５】
2004333224.tif 000017
【００３１】
【数１６】
2004333224.tif 000018
【００３２】
さらに、これら（１５）式に示したβ（ｒ０）を（１１）式に代入することにより、境界条件適合型応答関数が算定される。これら境界条件適合型応答関数は、一般に（１７）式を示すように、Ｇｒｅｅｎ関数Ｇの持つ性質である相反性を満足せず、Ｇｒｅｅｎ関数Ｇではない。しかし、境界条件適合型応答関数は、積分計算を不要としているため、計算負荷を大幅に低減できるものである。
【００３３】
【数１７】
2004333224.tif 000019
【００３４】
このように算定された境界条件適合型応答関数を、地盤面１に関するインピーダンス境界条件を満足する関数として、Ｇｒｅｅｎ関数Ｇに代わり導入することにより、地盤面に関する積分は不要となる。したがって、残る境界面гに着目し、（７）式に対してＧｒｅｅｎ関数Ｇに代わり境界条件適合型応答関数を導入した上で、境界面гに関するインピーダンス境界条件である（６）式を代入してすると、（１８）式を得る。観測点ｐが境界面г上にある場合には、（１８）式は、境界面г上の速度ポテンシャルφ（ｐ）を未知関数とする境界積分方程式となる。
【００３５】
【数１８】
2004333224.tif 000020
【００３６】
（ステップｓ４）
境界要素法では、解析領域Ωの境界面гを面要素に分割して、（１８）式で算定した境界積分方程式を離散化し数値解を求めることから、境界面гをＭ個の面要素に分割すると、（１８）式から（１９）式が算定される。なお、面要素は、対象とする周波数に応じた大きさに適宜分割すればよい。
これら境界面г上の境界積分方程式は、（１９）式を見てもわかるように、面要素毎に音響アドミッタンス比Ａｍを設定することができるため、吸音率が大きく変化する壁面にも対応することができるものである。
【００３７】
【数１９】
2004333224.tif 000021
【００３８】
上述する（１９）式は、速度ポテンシャルφ（ｐ）に関する複素係数連立１次方程式であり、本実施の形態では、境界面гをＭ個に分割しているので、Ｍ個の複素係数連立１次方程式が得られる。そこで、右辺の第１項をΨ０（ｐｌ）とおいて行列形式に書き改めると、（２０）式となる。
【００３９】
【数２０】
2004333224.tif 000022
【００４０】
ただし、ｌ＝ｍの場合には、音源と観測点が一致して境界条件適合型応答関数が発散するため、特異点処理を行う必要がある。（２０）式の行列要素ａｌｍ  （ｌ、ｍ＝１、２、・・・、Ｍ）は、特異点処理の要、不要に応じて場合分けすると、必要な場合は（２１）式、不要な場合は（２２）式で表すことができる。
【００４１】
【数２１】
2004333224.tif 000023
【数２２】
2004333224.tif 000024
【００４２】
ここで、特異点処理を必要とする（２１）式に示した行列要素ａｌｌ  について、従来より用いられているＬａｐｌａｃｅ方程式に関するＧｒｅｅｎ関数Ｇの特異点処理法に準じた特異点処理を示す。
図４に示すように、面要素∂гｌの中心ｐｌを原点として法線ベクトルｎをｚ負方向とする局所座標系ｘ’ｙ’ｚ’を考慮すると、（２３）式が成り立ち、境界条件適合型応答関数は、第１項の寄与が支配的となり、第２項の鏡像による寄与は無視できる。
【００４３】
【数２３】
2004333224.tif 000025
【００４４】
次に、面積Δｓｌを同じ面積の半径ａの円に置き換えると、境界条件適合型応答関数に関する面積分は（２４）式に近似できる。
【００４５】
【数２４】
2004333224.tif 000026
【００４６】
境界条件適合型応答関数の法線方向微分∂Ｇ／∂ｎも同様に、鏡像による寄与は無視できることから、（２５）式に近似できる。
【００４７】
【数２５】
2004333224.tif 000027
【００４８】
従って、（２１）式に示された、特異点処理を必要とする行列要素ａｌｌは、（２４）式及び（２５）式により（２６）式として表すことができる。
【００４９】
【数２６】
2004333224.tif 000028
【００５０】
（ステップｓ５）
以上より、（２０）式を解くことにより、境界面г上の速度ポテンシャル分布を求めることができる。境界面г上の速度ポテンシャル分布を、受音点ｐが解析領域Ω内にある場合の境界積分方程式（１８）式に代入することにより、解析領域Ω内の速度ポテンシャルが得られる。
【００５１】
上述する３次元音響波動解析方法の精度を検証するため、同一条件に対して境界条件適合型応答関数、Ｇｒｅｅｎ関数Ｇ、境界条件適合型応答関数と同様に積分計算が不要なＷｅｎｚｅｌ、及びＴｈｏｍａｓｓｏｎの近似式の、４つの算定式各々から速度ポテンシャルφ（ｐ）を算出してその結果を比較した。
【００５２】
なお、上述する３次元音響波動解析方法の精度を検証するために用いたＧｒｅｅｎ関数Ｇは、３次元音響波動解析方法について詳述する際に用いた（８）（９）（１０）式を満足し、また相反性を有するもので、（２７）式で表すことができる。
【００５３】
【数２７】
2004333224.tif 000029
【００５４】
【数２８】
2004333224.tif 000030
【００５５】
また、境界条件適合型応答関数と同じ積分計算が不要なＷｅｎｚｅｌ、及びＴｈｏｍａｓｓｏｎの近似式は、この（２７）式から音源ｑと観測点ｐの位置や、境界面гの音響アドミッタンス比Ａｇの値に応じて場合分けした近似解を求めている。
【００５６】
図５に示すように、音源ｑは、地上１．５ｍに設定し、周波数は１００Ｈｚに設定した。地盤面１の吸音率はαｇは、０．１とし、地盤面１からの反射音の位相のずれは０と仮定した。すなわち、吸音率αｇ  と音響アドミッタンス比Ａｇ  は、（２９）式で表すことができる。
【００５７】
【数２９】
2004333224.tif 000031
【００５８】
その他、種々の条件は、以下に示すとおりである。
観測点ｐの高さは、ｚ＝１．５、５．０、１０．０ｍの３条件とし、音源ｑと観測点ｐ間の水平距離ｒを０〜２０ｍまで変化させた場合の速度ポテンシャルφ（ｐ）を算定した。
Ｇｒｅｅｎ関数Ｇに含まれる半無限区間に関する積分は、有限区間｛０、１０００｝の積分に置き換えてＲｏｍｂｅｒｇ法により計算した。最大反復回数は２０回とし、積分値の相対誤差が１．０×１０−１０  以下になれば計算打ち切りとした。なお、Ｗｅｎｚｅｌ及びＴｈｏｍａｓｓｏｎの両近似式で特異点となるｒ＝０は計算から除外した。なお、計算には計測装置２としてパソコン（ＣＰＵ：Ｐｅｎｔｉｕｍ（登録商標）４１．８ＧＨｚ）を使用している。
【００５９】
これら４つの算定式より得られた速度ポテンシャルφ（ｐ）の算定結果を図６（ａ）から（ｃ）に示す。図６を見ると、境界条件適合型応答関数は、観測点ｐの高さによらずＧｒｅｅｎ関数Ｇの計算値に良く一致し、両者はほとんど重なっている。図６（ａ）に示すように、ｚ＝１．５ｍの場合でも、相対誤差は最大４．２％であり、音圧レベルに換算して０．４ｄＢ以下である。
一方、境界条件適合型応答関数  と同様で積分計算を不要とするＷｅｎｚｅｌ及びＴｈｏｍａｓｓｏｎの両近似式は、音源ｑと観測点ｐ間の水平距離ｒが大きい場合や観測点が地盤面１から離れる場合には、精度が急激に悪化している様子がわかる。
【００６０】
このように、境界条件適合型応答関数を用いた３次元音響波動解析方法は、Ｇｒｅｅｎ関数Ｇを用いた場合と比較して計算精度に大きな差異が無く、境界条件適合型応答関数と同様に積分計算を不要としているＷｅｎｚｅｌ、及びＴｈｏｍａｓｓｏｎの近似式と比較して、計算精度を大幅に向上することができるものである。
【００６１】
また、Ｗｅｎｚｅｌ、及びＴｈｏｍａｓｓｏｎの近似式で実施しているような、音源ｑと観測点ｐの位置や、境界面гの音響アドミッタンス比Ａｇの値に応じた場合分けが必要がなく、境界条件適合型応答関数を導入した（１８）式に示す境界積分方程式一つで、解析対象Ωの３次元音場を一様に表現できるため、広域の３次元音場解析への適用が可能となる。
【００６２】
なお、Ｗｅｎｚｅｌ、及びＴｈｏｍａｓｓｏｎの近似式や、本発明の３次元音響波動解析方法に適用する境界条件適合型応答関数のように、積分計算を不要とする構成は、計算負荷を大幅に低減することを目的とするものである。境界条件適合型応答関数と、（２７）式にも示したように、積分計算が必要となるＧｒｅｅｎ関数Ｇの両者について、１ケース（観測点数ｐ＝２００点）の計算に要した時間を算出すると、それぞれ０．０１６ｓ、２９４．２ｓであった。つまり、境界条件適合型応答関数を用いた本発明の３次元音響波動解析方法は、境界条件適合型応答関数に積分計算を不要とすることで、Ｇｒｅｅｎ関数Ｇと比較して、計算時間を約１／２００００以下に短縮することができるものである。
【００６３】
このように、本発明の３次元音響波動解析方法は、積分計算を不要とする境界条件適合型応答関数を導入することにより、計算速度が大幅に向上するだけでなく、従来のＧｒｅｅｎ関数Ｇを適用するような厳密解と比較しても計算精度に大きな差異が生じない程度の高い計算精度を有しているとともに、音源ｑと観測点ｐの位置や、境界面гの音響アドミッタンス比Ａｇの値に応じた場合分けの必要がなく、境界条件適合型応答関数を導入した境界積分方程式一つで、解析対象Ωの３次元音場を一様に表現できるため、広域な３次元空間に対して高精度な波動音響シミュレーションによる検討が実施でき、従来より実験及び実測に必要だった費用と工数を大幅に削減することが可能となる。
【００６４】
また、波動性を考慮した波動音響解析方法を適用しているため、音波の波動性が問題となるような低周波数音問題、例えば建物屋上に設置された設備機器から発生する低周波数音が周辺地域に及ぼす影響を、高精度に解析することが可能となる。これにより、騒音実体の把握、騒音対策効果の事前評価に活用でき、騒音対策の重点化を図ることが可能となる。
【００６５】
【発明の効果】
請求項１記載の３次元音響波動解析方法によれば、境界要素法による波動音響解析法を介して、少なくとも地盤面及び建物群の外壁面に境界面を有する広域な３次元空間を対象領域とした３次元音場を、速度ポテンシャルを用いて定式化する３次元音響波動解析方法であって、速度ポテンシャルを解とする支配方程式に、波動場の波源に対する応答を表す関数であり、前記支配方程式の地盤面に関するインピーダンス境界条件と同一の境界条件を満足する関数を導入して、前記支配方程式を、前記対象領域内で一様な境界積分方程式に変換する。
【００６６】
これにより、計算速度が大幅に向上するだけでなく、従来のＧｒｅｅｎ関数Ｇを適用するような厳密解と比較しても計算精度に大きな差異が生じない程度の高い計算精度を有しているとともに、音源と観測点の位置や、境界面の音響アドミッタンス比の値に応じた場合分けが必要がなく、境界積分方程式一つで、解析対象の３次元音場を一様に表現できるため、広域な３次元空間に対して高精度な波動音響シミュレーションによる検討が実施でき、従来より実験及び実測に必要だった費用と工数を大幅に削減することが可能となる。
【００６７】
また、波動性を考慮した波動音響解析方法を適用しているため、音波の波動性が問題となるような低周波数音問題、例えば建物屋上に設置された設備機器から発生する低周波数音が周辺地域に及ぼす影響を、高精度に解析することが可能となる。これにより、騒音実体の把握、騒音対策効果の事前評価に活用でき、騒音対策の重点化を図ることが可能となる。
【図面の簡単な説明】
【図１】本発明に係る３次元音響波動解析方法の行程を示すフロー図である。
【図２】本発明に係る３次元音場を定式化する解析領域Ωを示す図である。
【図３】本発明に係る境界条件適合型応答関数導出のための鏡像の設定を示す図である。
【図４】本発明に係る特異点処理のための局所座標系を示す図である。
【図５】本発明に係る境界条件適合型応答関数の精度検証のための計算条件を示す図である。
【図６】本発明に係る速度ポテンシャル計算値の比較グラフを示す図である。
【符号の説明】
１  地盤面
２  計測装置
３  パワーアンプ
４  マイクアンプ
【図１】
2004333224.tif 000032
【図２】
2004333224.tif 000033
【図３】
2004333224.tif 000034
【図４】
2004333224.tif 000035
【図５】
2004333224.tif 000036
【図６】
2004333224.tif 000037