IP Force 特許公報掲載プロジェクト 2022.1.31 β版

知財求人 - 知財ポータルサイト「IP Force」

▶ 株式会社安藤・間の特許一覧 ▶ 国立大学法人神戸大学の特許一覧

特許7285514音場解析装置、音場解析方法及びプログラム
<>
  • 特許-音場解析装置、音場解析方法及びプログラム 図1
  • 特許-音場解析装置、音場解析方法及びプログラム 図2
  • 特許-音場解析装置、音場解析方法及びプログラム 図3
  • 特許-音場解析装置、音場解析方法及びプログラム 図4
  • 特許-音場解析装置、音場解析方法及びプログラム 図5
  • 特許-音場解析装置、音場解析方法及びプログラム 図6
  • 特許-音場解析装置、音場解析方法及びプログラム 図7
  • 特許-音場解析装置、音場解析方法及びプログラム 図8
  • 特許-音場解析装置、音場解析方法及びプログラム 図9
  • 特許-音場解析装置、音場解析方法及びプログラム 図10
  • 特許-音場解析装置、音場解析方法及びプログラム 図11
< >
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2023-05-25
(45)【発行日】2023-06-02
(54)【発明の名称】音場解析装置、音場解析方法及びプログラム
(51)【国際特許分類】
   G01H 3/00 20060101AFI20230526BHJP
   G06F 30/23 20200101ALI20230526BHJP
【FI】
G01H3/00 A
G06F30/23
【請求項の数】 7
(21)【出願番号】P 2019135950
(22)【出願日】2019-07-24
(65)【公開番号】P2021018221
(43)【公開日】2021-02-15
【審査請求日】2022-04-19
【新規性喪失の例外の表示】特許法第30条第2項適用 発行者名:一般社団法人日本音響学会、刊行物名:日本音響学会2018年秋季研究発表会講演論文集、発行年月日:平成30年8月29日 集会名:日本音響学会2018年秋季研究発表会、開催日:平成30年9月12日 発行者名:株式会社安藤・間 技術研究所、刊行物名:安藤ハザマ研究年報 Vol.6 2018、発行年月日:平成31年2月20日
(73)【特許権者】
【識別番号】303057365
【氏名又は名称】株式会社安藤・間
(73)【特許権者】
【識別番号】504150450
【氏名又は名称】国立大学法人神戸大学
(74)【代理人】
【識別番号】100095407
【弁理士】
【氏名又は名称】木村 満
(74)【代理人】
【識別番号】100098246
【弁理士】
【氏名又は名称】砂場 哲郎
(74)【代理人】
【識別番号】100132883
【弁理士】
【氏名又は名称】森川 泰司
(74)【代理人】
【識別番号】100181618
【弁理士】
【氏名又は名称】宮脇 良平
(72)【発明者】
【氏名】吉田 卓彌
(72)【発明者】
【氏名】奥園 健
(72)【発明者】
【氏名】阪上 公博
【審査官】奥野 尭也
(56)【参考文献】
【文献】特開2017-111749(JP,A)
【文献】特開2014-119882(JP,A)
【文献】米国特許出願公開第2015/0377915(US,A1)
【文献】吉田卓彌、奥園健、阪上公博,長方形要素を用いた陽的時間領域有限要素法による室内音場解析のための修正積分則,日本音響学会誌,日本,一般社団法人 日本音響学会,2017年07月01日,72巻7号,pp. 367-373,DOI:https://doi.org/10.20697/jasj.72.7_367
(58)【調査した分野】(Int.Cl.,DB名)
G01H 1/00-17/00
(57)【特許請求の範囲】
【請求項1】
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を離散化する要素の形状に依存しない数値積分点を用いて計算された質量行列及び剛性行列を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析装置であって、
前記一階常微分方程式において時間の離散化誤差を低減するために最適化された線形多段法を適用することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において前記微分ベクトルの時間一階微分を1次精度の差分近似で離散化することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する数値計算部、を備え、
前記線形多段法は、Adams-Bashforth型の時間積分法であり、
前記音圧ベクトルの時間更新式は、第n時間ステップでの前記音圧ベクトルの値を、第(n-1)時間ステップでの前記音圧ベクトルの値と、第(n-4)時間ステップから第(n-1)時間ステップの前記微分ベクトルの値を、前記時間の離散化誤差に関して4次精度を達成する値に最適化された重み係数を乗じた上で加算することにより計算されるベクトルの値と、により定める式である、
音場解析装置。
【請求項2】
前記音圧ベクトルの時間更新式は、第n時間ステップでの前記音圧ベクトルの値を、第(n-1)時間ステップでの前記音圧ベクトルの値と、第(n-4)時間ステップから前記第(n-1)時間ステップの前記微分ベクトルの値と、により定める式であり、
前記微分ベクトルの時間更新式は、前記第n時間ステップでの前記微分ベクトルの値を、前記第(n-1)時間ステップでの前記微分ベクトルの値と、前記第n時間ステップでの前記音圧ベクトルの値と、により定める式である、
請求項1に記載の音場解析装置。
【請求項3】
前記数値計算部は、
前記音圧ベクトルの時間更新式に従って、前記第n時間ステップでの前記音圧ベクトルの値を、前記第(n-1)時間ステップでの前記音圧ベクトルの値と、前記第(n-4)時間ステップから前記第(n-1)時間ステップの前記微分ベクトルの値と、から計算する音圧ベクトル計算部と、
前記微分ベクトルの時間更新式に従って、前記第n時間ステップでの前記微分ベクトルの値を、前記第(n-1)時間ステップでの前記微分ベクトルの値と、前記第n時間ステップでの前記音圧ベクトルの値と、から計算する微分ベクトル計算部と、
nの値を1ずつ増加させながら、前記音圧ベクトル計算部により前記第n時間ステップでの前記音圧ベクトルの値を計算する処理と、前記微分ベクトル計算部により前記第n時間ステップでの前記微分ベクトルの値を計算する処理と、を繰り返す繰り返し部と、を備える、
請求項2に記載の音場解析装置。
【請求項4】
前記数値計算部により計算された前記複数の時間ステップでの前記音圧ベクトルの値に基づく出力情報を出力する出力部、を更に備える、
請求項1からのいずれか1項に記載の音場解析装置。
【請求項5】
前記出力部は、前記出力情報として、前記音圧の時間変化と、前記音圧の空間分布と、のうちの少なくとも一方を表す情報を出力する、
請求項に記載の音場解析装置。
【請求項6】
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を離散化する要素の形状に依存しない数値積分点を用いて計算された質量行列及び剛性行列を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析方法であって、
前記一階常微分方程式において時間の離散化誤差を低減するために最適化された線形多段法を適用することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において前記微分ベクトルの時間一階微分を1次精度の差分近似で離散化することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算
前記線形多段法は、Adams-Bashforth型の時間積分法であり、
前記音圧ベクトルの時間更新式は、第n時間ステップでの前記音圧ベクトルの値を、第(n-1)時間ステップでの前記音圧ベクトルの値と、第(n-4)時間ステップから第(n-1)時間ステップの前記微分ベクトルの値を、前記時間の離散化誤差に関して4次精度を達成する値に最適化された重み係数を乗じた上で加算することにより計算されるベクトルの値と、により定める式である、
音場解析方法。
【請求項7】
コンピュータを、
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を離散化する要素の形状に依存しない数値積分点を用いて計算された質量行列及び剛性行列を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析装置として機能させるプログラムであって、
前記コンピュータを、
前記一階常微分方程式において時間の離散化誤差を低減するために最適化された線形多段法を適用することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において前記微分ベクトルの時間一階微分を1次精度の差分近似で離散化することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する数値計算部、として機能させ、
前記線形多段法は、Adams-Bashforth型の時間積分法であり、
前記音圧ベクトルの時間更新式は、第n時間ステップでの前記音圧ベクトルの値を、第(n-1)時間ステップでの前記音圧ベクトルの値と、第(n-4)時間ステップから第(n-1)時間ステップの前記微分ベクトルの値を、前記時間の離散化誤差に関して4次精度を達成する値に最適化された重み係数を乗じた上で加算することにより計算されるベクトルの値と、により定める式である、
プログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、音場解析装置、音場解析方法及びプログラムに関する。
【背景技術】
【0002】
数値シミュレーションにより音場を解析する手法が知られている。例えば、非特許文献1~3は、一階常微分方程式(ODE:Ordinary Differential Equation)に基づく時間領域有限要素法(TDFEM:Time-Domain Finite Element Method)を用いて音場を解析する手法を開示している。具体的に説明すると、非特許文献1は、以下の(手順1)~(手順4)に従って音場を解析する手法を開示している。
【0003】
(手順1)
剛壁境界又は振動境界により囲われた音場を対象とした場合の半離散化方程式に、対角質量行列と音圧の時間一階微分ベクトルとを導入することで得られる連立一階常微分方程式を対象とする。連立一階常微分方程式は、具体的には、解析対象の空間における音圧の分布を示すベクトルp(以下、“音圧ベクトルp”と呼ぶ。)と、音圧ベクトルpの時間一階微分と等価なベクトルv(以下、“微分ベクトルv”と呼ぶ。)と、を未知数として、下記(1)式及び(2)式のように表される。
【0004】
【数1】
【0005】
【数2】
【0006】
ここで、(1)式及び(2)式における“M”、“K”、及び“D”は、それぞれ質量行列、剛性行列、及び対角質量行列を表す。質量行列M、剛性行列K、及び対角質量行列Dは、それぞれ、要素質量行列M、要素剛性行列K、及び要素対角質量行列Dを重ね合わせて構成される全体行列である。また、(1)式及び(2)式において、 “f”は解析対象の空間に加えられる外力の分布を示す外力ベクトルを表し、“c”は音速を表し、“・”は時間に関する一階微分を表す。
【0007】
(手順2)
空気領域の有限要素として正方形又は立方体形状の線形一次要素を設定する。そして、設定された空気領域の有限要素に対して、非特許文献2に開示された修正積分則を適用して、(1)式及び(2)式における質量行列Mと剛性行列Kとを計算する。非特許文献2に開示された修正積分則では、それぞれ要素質量行列M及び要素剛性行列KをGauss-Legendre求積法により計算するための数値積分点α及びαを、従来の値である±√(1/3)から下記(3)式のように修正する。これにより、時間及び空間の離散化を起因とする誤差である離散化誤差に関して4次精度を達成する。
【0008】
【数3】
【0009】
(3)式において、“τ”はクーラン数を表す。クーラン数τは、音速cと、時間領域有限要素法における時間刻み幅Δtと、正方形又は立方体形状要素の要素長hと、を用いて、“τ=cΔt/h”と定められる。
【0010】
(手順3)
(1)式及び(2)式で表される連立一階常微分方程式を時間方向に離散化することで時間進行スキームを構築し、音波の時間進行を計算する。
【0011】
(手順4)
(手順3)で音波の時間進行を計算する際に、数値的安定性を保つため、(1)式の音圧ベクトルpの時間一階微分を1次精度の前進差分近似で離散化し、更に(2)式の微分ベクトルvの時間一階微分を1次精度の後退差分近似で離散化する。この結果として、下記(4)式及び(5)式が得られる。
【0012】
【数4】
【0013】
【数5】
【0014】
ここで、(4)式及び(5)式における“n”は時間ステップ数を表す自然数であり、“Δt”は時間刻み幅を表す。非特許文献1に開示された手法では、(4)式及び(5)式により表される時間更新式に従って、音波の時間進行を計算する。
【0015】
このような一階常微分方程式に基づく時間領域有限要素法は、計算効率に関する長所として次の2点を有する。
(I)時間進行スキームが陽的であるため、連立一次方程式の求解が不要である。そのため、時間ステップあたりの計算負荷が少ない。
(II)修正積分則の適用により、使用する要素が一次要素であるにもかかわらず、時間及び空間で4次の精度を達成することができる。
【0016】
また、非特許文献3は、一階常微分方程式に基づく時間領域有限要素法を用いた音場解析の手法において、正方形又は立方体以外の任意の形状の要素を用いた場合に離散化誤差を空間4次及び時間2次の精度で低減することが可能な要素行列の数値積分点の設定法を、理論解析及び数値実験により明らかにしている。
【先行技術文献】
【非特許文献】
【0017】
【文献】Takumi Yoshida, Takeshi Okuzono, Kimihiro Sakagami, “Numerically stable explicit time-domain finite element method for room acoustics simulation using an equivalent impedance model”, Noise Control Engineering Journal, 66(3), 176-189 (1 May2018)
【文献】B. Yue, MN. Guddati, “Dispersion-reducing finite elements for transient acoustics”, Journal of Acoustical Society of America, 118(4), 2132-2141(2005)
【文献】吉田卓彌、奥園健、阪上公博、"長方形要素を用いた陽的時間領域有限要素法による室内音場解析のための修正積分則"、日本音響学会誌72巻7号(2016), pp.367-373
【発明の概要】
【発明が解決しようとする課題】
【0018】
数値シミュレーションにより音場を解析する上で、計算効率の向上は重要である。非特許文献1に開示された一階常微分方程式に基づく時間領域有限要素法は、非特許文献2に開示された修正積分則の適用により、優れた計算効率を有する。その一方で、修正積分則を適用するためには、使用できる要素の形状が正方形又は立方体に限定されるという制約がある。そのため、解析対象の空間が複雑な形状を含む場合には、形状の近似精度に起因する誤差を生じる恐れがある。
【0019】
これに対して、非特許文献3に開示された手法では、任意の形状の要素を扱うことが可能である。その一方で、空間の離散化誤差に関して精度が4次であるのに対して時間の離散化誤差に関して精度が2次であるため、時間の精度が低いという課題がある。このような状況のもと、時間領域有限要素法を用いた音場解析において、要素の形状を限定せずに、時間及び空間の離散化誤差に関して高い精度を達成することが可能な手法が求められている。
【0020】
本発明は、以上の課題を解決するためのものであり、時間領域有限要素法を用いた音場解析において、要素の形状を限定せずに、時間及び空間の離散化誤差に関して高い精度を達成することが可能な音場解析装置、音場解析方法及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【0021】
上記目的を達成するために、本発明の第1の観点に係る音場解析装置は、
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を離散化する要素の形状に依存しない数値積分点を用いて計算された質量行列及び剛性行列を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析装置であって、
前記一階常微分方程式において時間の離散化誤差を低減するために最適化された線形多段法を適用することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において前記微分ベクトルの時間一階微分を1次精度の差分近似で離散化することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する数値計算部、を備え
前記線形多段法は、Adams-Bashforth型の時間積分法であり、
前記音圧ベクトルの時間更新式は、第n時間ステップでの前記音圧ベクトルの値を、第(n-1)時間ステップでの前記音圧ベクトルの値と、第(n-4)時間ステップから第(n-1)時間ステップの前記微分ベクトルの値を、前記時間の離散化誤差に関して4次精度を達成する値に最適化された重み係数を乗じた上で加算することにより計算されるベクトルの値と、により定める式である。
【0022】
前記音圧ベクトルの時間更新式は、第n時間ステップでの前記音圧ベクトルの値を、第(n-1)時間ステップでの前記音圧ベクトルの値と、第(n-4)時間ステップから前記第(n-1)時間ステップの前記微分ベクトルの値と、により定める式であり、
前記微分ベクトルの時間更新式は、前記第n時間ステップでの前記微分ベクトルの値を、前記第(n-1)時間ステップでの前記微分ベクトルの値と、前記第n時間ステップでの前記音圧ベクトルの値と、により定める式であっても良い。
【0023】
前記数値計算部は、
前記音圧ベクトルの時間更新式に従って、前記第n時間ステップでの前記音圧ベクトルの値を、前記第(n-1)時間ステップでの前記音圧ベクトルの値と、前記第(n-4)時間ステップから前記第(n-1)時間ステップの前記微分ベクトルの値と、から計算する音圧ベクトル計算部と、
前記微分ベクトルの時間更新式に従って、前記第n時間ステップでの前記微分ベクトルの値を、前記第(n-1)時間ステップでの前記微分ベクトルの値と、前記第n時間ステップでの前記音圧ベクトルの値と、から計算する微分ベクトル計算部と、
nの値を1ずつ増加させながら、前記音圧ベクトル計算部により前記第n時間ステップでの前記音圧ベクトルの値を計算する処理と、前記微分ベクトル計算部により前記第n時間ステップでの前記微分ベクトルの値を計算する処理と、を繰り返す繰り返し部と、を備えても良い。
【0025】
前記数値計算部により計算された前記複数の時間ステップでの前記音圧ベクトルの値に基づく出力情報を出力する出力部、を更に備えても良い。
【0026】
前記出力部は、前記出力情報として、前記音圧の時間変化と、前記音圧の空間分布と、のうちの少なくとも一方を表す情報を出力しても良い。
【0027】
上記目的を達成するために、本発明の第2の観点に係る音場解析方法は、
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を離散化する要素の形状に依存しない数値積分点を用いて計算された質量行列及び剛性行列を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析方法であって、
前記一階常微分方程式において時間の離散化誤差を低減するために最適化された線形多段法を適用することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において前記微分ベクトルの時間一階微分を1次精度の差分近似で離散化することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算し、
前記線形多段法は、Adams-Bashforth型の時間積分法であり、
前記音圧ベクトルの時間更新式は、第n時間ステップでの前記音圧ベクトルの値を、第(n-1)時間ステップでの前記音圧ベクトルの値と、第(n-4)時間ステップから第(n-1)時間ステップの前記微分ベクトルの値を、前記時間の離散化誤差に関して4次精度を達成する値に最適化された重み係数を乗じた上で加算することにより計算されるベクトルの値と、により定める式である。
【0028】
上記目的を達成するために、本発明の第3の観点に係るプログラムは、
コンピュータを、
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を離散化する要素の形状に依存しない数値積分点を用いて計算された質量行列及び剛性行列を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析装置として機能させるプログラムであって、
前記コンピュータを、
前記一階常微分方程式において時間の離散化誤差を低減するために最適化された線形多段法を適用することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において前記微分ベクトルの時間一階微分を1次精度の差分近似で離散化することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する数値計算部、として機能させ
前記線形多段法は、Adams-Bashforth型の時間積分法であり、
前記音圧ベクトルの時間更新式は、第n時間ステップでの前記音圧ベクトルの値を、第(n-1)時間ステップでの前記音圧ベクトルの値と、第(n-4)時間ステップから第(n-1)時間ステップの前記微分ベクトルの値を、前記時間の離散化誤差に関して4次精度を達成する値に最適化された重み係数を乗じた上で加算することにより計算されるベクトルの値と、により定める式である。
【発明の効果】
【0029】
本発明によれば、時間領域有限要素法を用いた音場解析において、要素の形状を限定せずに、時間及び空間の離散化誤差に関して高い精度を達成できる。
【図面の簡単な説明】
【0030】
図1】本発明の実施形態に係る音場解析装置の構成を示すブロック図である。
図2】実施形態におけるメッシュの設定例を示す図である。
図3】実施形態に係る音場解析装置における数値計算部の構成を示すブロック図である。
図4】実施形態に係る音場解析装置によって得られた音圧の時間変化の表示例を示す図である。
図5】実施形態に係る音場解析装置によって得られた音圧の空間分布の表示例を示す図である。
図6】実施形態に係る音場解析装置によって実行される処理の流れを示すフローチャートである。
図7】本発明の第1の数値実験において解析対象として使用した2次元の音場を示す図である。
図8】(a)は、第1の数値実験において、従来の手法により計算された音圧の時間波形を、参照解の時間波形と共に示す図である。(b)は、第1の数値実験において、実施形態の手法により計算された音圧の時間波形を、参照解の時間波形と共に示す図である。
図9】第1の数値実験において、一般的なAdams法により計算された音圧の時間波形を、参照解の時間波形と共に示す図である。
図10】本発明の第2の数値実験において解析対象として使用した2次元の音場を示す図である。
図11】第2の数値実験において、実施形態の手法により計算された音圧の時間波形を、参照解の時間波形と共に示す図である。
【発明を実施するための形態】
【0031】
以下、本発明の実施形態について、図面を参照して説明する。なお、図中同一又は相当する部分には同一符号を付す。
【0032】
(実施形態)
本実施形態に係る音場解析装置10は、時間領域有限要素法を用いた波動音響数値シミュレーションにより音場を解析し、これにより音響分野における材料や測定法の開発、設計支援等に活用することが可能な装置である。ここで、音場とは、音波が存在する空間を意味する。
【0033】
図1に、音場解析装置10の構成を示す。図1に示すように、音場解析装置10は、制御部11と、記憶部12と、操作部13と、表示部14と、通信部15と、を備える。
【0034】
制御部11は、音場解析装置10の各部と制御バスを介して接続されており、音場解析装置10全体の動作を統括制御する。具体的に説明すると、制御部11は、CPU(Central Processing Unit)、ROM(Read Only Memory)及びRAM(Random Access Memory)を備える。CPUは、例えばマイクロプロセッサ等であって、様々な処理及び演算を実行する中央演算処理部である。制御部11において、CPUが、ROMに記憶されている制御プログラムを読み出して、RAMをワークメモリとして用いながら、各種の処理を実行する。
【0035】
記憶部12は、フラッシュメモリ、ハードディスク等の不揮発性メモリである。記憶部12は、オペレーションシステム、アプリケーションプログラム等の、制御部11が各種処理を行うために使用するプログラム及びデータを記憶する。また、記憶部12は、制御部11が各種処理を行うことにより生成又は取得するデータを記憶する。
【0036】
操作部13は、キーボード、マウス、ボタン、タッチパッド、タッチパネル等の入力デバイスを備えており、ユーザから操作を受け付ける。ユーザは、操作部13を操作することによって、様々な指示を音場解析装置10に入力することができる。操作部13は、ユーザから入力された操作指示を受け付けると、受け付けた操作指示を制御部11に送信する。
【0037】
表示部14は、液晶ディスプレイ、有機EL(Electro Luminescence)ディスプレイ等の表示デバイスを備える。表示部14は、図示しない表示駆動回路によって駆動され、制御部11による制御のもとで様々な画像を表示する。例えば、表示部14は、数値シミュレーションにより音場を解析した結果を表す画像を表示する。
【0038】
通信部15は、音場解析装置10が外部の機器と通信するための通信インタフェースを備える。通信部15は、例えばインターネット等の広域ネットワーク、LAN(Local Area Network)、USB(Universal Serial Bus)等の有線又は無線による通信を介して、外部の機器と通信する。
【0039】
制御部11は、図1に示すように、機能的に、設定受付部110と、数値計算部120と、出力部130と、を備える。制御部11において、CPUがROMに記憶されたプログラムをRAMに読み出して実行することにより、これら各部として機能する。
【0040】
設定受付部110は、音場解析装置10において音場を解析するための設定を受け付ける。ここで、設定は、数値計算部120が数値シミュレーションにより音圧の分布を計算するための条件やパラメータ等であって、具体的には解析対象の空間、境界条件、メッシュ等に関する情報である。設定受付部110は、ユーザが所望の設定を入力するための設定画面を表示部14に表示する。ユーザは、表示部14に表示された設定画面を見ながら、操作部13を操作して、設定を入力する。設定受付部110は、ユーザにより入力された設定を受け付ける。
【0041】
第1に、設定受付部110は、音場解析の目的物である解析対象の空間の設定を受け付ける。具体的に説明すると、設定受付部110は、操作部13を介してユーザから入力された操作に従って、解析対象の空間の次元及び形状の設定を受け付ける。
【0042】
例えば、長方形状の2次元平面内における音場を解析する場合、設定受付部110は、解析対象の空間として、長方形状の2次元の空間の設定を受け付ける。或いは、3次元の室内空間内における音場を解析する場合には、設定受付部110は、解析対象の空間として、その室内空間に対応する形状の3次元の空間の設定を受け付ける。また、3次元の音響材料の内部における音場を解析する場合には、設定受付部110は、解析対象の空間として、その音響材料に対応する形状の3次元の空間の設定を受け付ける。
【0043】
第2に、設定受付部110は、境界条件の設定を受け付ける。境界条件は、解析対象の空間を囲む境界に関する条件である。設定受付部110は、操作部13を介してユーザから入力された操作に従って、境界条件として、剛壁境界と振動境界とのうちのいずれかの設定を受け付ける。
【0044】
剛壁境界とは、音響的に剛である境界であって、境界上における粒子の速度が0である境界である。剛壁境界に入射した音波は、吸収されずに完全に反射される。これに対して、振動境界とは、境界上における粒子が振動可能な境界である。
【0045】
第3に、設定受付部110は、メッシュの設定を受け付ける。ここで、メッシュとは、解析対象の空間を数値的に解析するために、解析対象の空間を四角形等の単純な形状をした複数の有限要素に分割(離散化)したものである。なお、有限要素は、単に“要素”とも呼ぶ。設定受付部110は、操作部13を介してユーザから入力された操作に従って、メッシュを構成する複数の要素の形状及びサイズの設定を受け付ける。
【0046】
図2に、2次元の解析対象の空間20を四角形状の要素で分割した例を示す。解析対象の空間20が2次元である場合、解析対象の空間20は、例えば図2に示すように、1辺の長さhの正方形状の要素によって、横方向にX個、縦方向にY個に分割される。すなわち、解析対象の空間20は、合計(X×Y)個の要素で離散化される。或いは、解析対象の空間が3次元である場合には、図示は省略するが、解析対象の空間は、立方体、直方体等の3次元形状の要素によって、横方向、縦方向及び高さ方向に分割される。
【0047】
音場解析装置10は、このように離散化された要素の単位で音圧の値を計算することにより、解析対象の空間における音圧の分布を解析する。要素のサイズを小さく設定するほど、空間の離散化に伴う誤差が小さくなるため解析精度が向上するが、計算量が増大する。
【0048】
なお、後述するように、本実施形態において、数値計算部120は、要素の形状に依存しない手法を用いて音圧ベクトルpの値を計算する。そのため、要素の形状は、正方形に限らず、任意形状の四辺形であっても良いし、任意形状の六面体であっても良い。
【0049】
第4に、設定受付部110は、時間刻み幅Δt、総ステップ数N、及び外力ベクトルfの設定を受け付ける。時間刻み幅Δtは、解析対象の空間における音波の時間進行を計算するための時間ステップの単位である。時間刻み幅Δtを小さく設定するほど、時間の離散化に伴う誤差が小さくなるため解析精度が向上するが、計算量が増大する。総ステップ数Nは、解析対象の空間における音波の時間進行を初期状態から何ステップ先まで計算するかを示す値である。数値計算部120は、Δt×Nの時間区間に亘って、音圧を計算する。外力ベクトルfは、計算開始から計算終了までの時間区間において、解析対象の空間に加えられる外力の分布を示すベクトルである。
【0050】
このように、設定受付部110は、時間領域有限要素法により音場を解析するための設定を、操作部13を介してユーザから受け付ける。設定受付部110は、制御部11が操作部13等と協働することにより実現される。
【0051】
図1に戻って、数値計算部120は、設定受付部110により受け付けられた設定のもとで、数値計算を実行する。具体的に説明すると、数値計算部120は、解析対象の空間における音圧の分布を示す音圧ベクトルpと、音圧ベクトルpの時間一階微分と等価な微分ベクトルvと、に関する一階常微分方程式から得られる音圧ベクトルpの時間更新式と微分ベクトルvの時間更新式とに従って、複数の時間ステップでの音圧ベクトルpの値を計算する。一階常微分方程式は、(1)式及び(2)式に示したように、質量行列M及び剛性行列Kを含む式である。
【0052】
ここで、数値計算部120による計算手法について説明する。本実施形態において、数値計算部120は、メッシュを構成する要素の形状に制約されずに離散化誤差の低減を図るため、非特許文献1に開示された(4)式とは異なる時間更新式に従って、数値計算を実行する。具体的に説明すると、非特許文献2に開示された修正積分則では、要素の形状が正方形又は立方体に限定される。これに対して、本実施形態では、正方形又は立方体に限定されない任意形状の四辺形または六面体要素で離散化できるようにするため、要素の形状に依存しない数値積分点を使用する。具体的には、非特許文献3に開示された手法と同様に、数値積分点α及びαを、“α=±√(4/3)”及び“α=±√(2/3)”と設定する。(1)式及び(2)式に含まれる質量行列M及び剛性行列Kは、このように要素の形状に依存しない数値積分点α及びαを用いて計算される。
【0053】
但し、非特許文献3に開示された手法では、離散化誤差を低減する精度に関して、空間精度は4次で達成できるが、時間精度は2次までしか達成できなかった。これに対して、本実施形態では、空間と時間のどちらも4次の精度で離散化誤差を低減することを目的に、時間の離散化誤差を低減するために最適化された線形多段法を適用することにより、(1)式における音圧ベクトルpの時間一階微分を時間方向に離散化する。線形多段法は、過去の複数の時間ステップでの値を用いて次の時間ステップでの値を計算する数値解法である。
【0054】
より詳細には、数値計算部120は、線形多段法として、Adams-Bashforth型の時間積分法(以下、“Adams法”と呼ぶ。)を用いる。具体的には、時間4次精度のAdams法を適用することにより、下記(6)式のように表される音圧ベクトルpの時間更新式を導出する。(6)式は、音圧ベクトルpの時間勾配を、過去4つの時間ステップにおける微分ベクトルvを用いて表現した式である。なお、Adams法の詳細は、例えば“Philipp O.J. Scherer, Computational Physics -Simulation of Classical and Quantum Systems-, Springer, ed. 3(2017)”に開示されている。
【0055】
【数6】
【0056】
ここで、(6)式における“b”(n=1,2,3,4)は、重み係数を表す。安定した計算のためには、重み係数b(n=1,2,3,4)は、下記(7)式を満たす、すなわち重み係数bの総和が1になることが必要条件となる。(6)式における重み係数b(n=1,2,3,4)を、(7)式の条件を満たした上で適切に設定することで、時間積分法を最適化することができる。
【0057】
【数7】
【0058】
(5)式及び(6)式について、下記の(8)式で表される二次元自由空間における平面波の伝搬を仮定し、非特許文献2に開示された分散誤差解析を実施する。その結果として、下記の(9)式に示すように、数値的な音速cと厳密な音速cとの相対誤差として定義される分散誤差の式が得られる。なお、(8)式において、“k”は波数を表し、“θ”は音波の伝搬方向角を表す。また、“i”は虚数単位である。“x”、“y”は二次元直交座標系における任意の座標を、“ω”は角周波数を表す。
【0059】
【数8】
【0060】
【数9】
【0061】
ここで、(9)式における係数F,F,Fは、それぞれ下記の(10)式、(11)式及び(12)式のように表される。
【0062】
【数10】
【0063】
【数11】
【0064】
【数12】
【0065】
ここで、(9)式において、“b=14/12”、“b=-5/12”、“b=4/12”及び“b=-1/12”とすると、係数F,F,Fのそれぞれを消去することができる。その結果、下記(13)式のように、離散化誤差に関して時間及び空間で4次の精度を達成することができる。
【0066】
【数13】
【0067】
従って、時空間4次精度を達成するAdams法を用いた音圧ベクトルpの時間更新式は、下記の(14)式のように表される。
【0068】
【数14】
【0069】
数値計算部120は、このように線形多段法を適用することで得られた音圧ベクトルpの時間更新式である(14)式と、微分ベクトルvの時間更新式である(5)式と、に従って、複数の時間ステップでの音圧ベクトルpの値を計算する。(5)式は、背景技術で説明したように、(2)式において微分ベクトルvの時間一階微分を1次精度の差分近似で離散化することにより得られる微分ベクトルvの時間更新式である。
【0070】
図3に、数値計算部120のより詳細な構成を示す。図3に示すように、数値計算部120は、行列設定部121と、音圧ベクトル計算部123と、微分ベクトル計算部125と、繰り返し部127と、の機能を有する。
【0071】
行列設定部121は、(14)式と(5)式とに含まれる質量行列M、剛性行列K、及び対角質量行列Dを設定する。質量行列M、剛性行列K、及び対角質量行列Dは、解析対象の空間を離散化する要素の数に対応するサイズの行列であって、各行列に含まれる要素の値は、解析対象の空間の次元及び形状と、離散化する要素の形状と、に依存して定められる。そのため、行列設定部121は、設定受付部110によりユーザから受け付けられた設定に応じて、質量行列M、剛性行列K、及び対角質量行列Dを設定する。これにより、行列設定部121は、音圧ベクトルp及び微分ベクトルvの時間更新式を確立する。
【0072】
より詳細には、行列設定部121は、非特許文献3に開示された手法と同様に、質量行列M、剛性行列K、及び対角質量行列Dを設定する。具体的に説明すると、行列設定部121は、メッシュの要素長hに依存しない数値積分点“α=±√(4/3)”及び“α=±√(2/3)”を適用する。そして、行列設定部121は、この数値積分点α及びαを用いて、Gauss-Legendre求積法により要素質量行列M及び要素剛性行列Kを計算し、質量行列M及び剛性行列Kを設定する。ここで、数値積分点α及びαから要素質量行列M及び要素剛性行列Kを計算し、更に質量行列M及び剛性行列Kを得る手法は、非特許文献1,2等に開示されている周知の手法を用いることができる。また、行列設定部121は、非特許文献1,2等に開示されている周知の手法を用いて、対角行列である要素対角質量行列Dを計算し、対角質量行列Dを設定する。
【0073】
音圧ベクトル計算部123は、音圧ベクトルpの時間更新式である(14)式に従って、第n時間ステップでの音圧ベクトルpを計算する。ここで、第n時間ステップでの音圧ベクトルpは、解析対象の空間に設定された複数の要素のそれぞれにおける、計算開始時点からn番目の時間ステップでの音圧の値を要素として有するベクトルである。
【0074】
上述したように、(14)式は、第n時間ステップでの音圧ベクトルpの値を、第(n-1)時間ステップでの音圧ベクトルpn-1の値と、第n時間ステップでの微分ベクトルvの値と、により定める式である。従って、音圧ベクトル計算部123は、第n時間ステップでの音圧ベクトルpの値を、第(n-1)時間ステップでの音圧ベクトルpn-1の値と、第(n-4)時間ステップから第(n-1)時間ステップの微分ベクトルvn-1,vn-2,vn-3,vn-4の値と、から計算する。
【0075】
具体的に説明すると、音圧ベクトル計算部123は、第(n-1)時間ステップでの音圧ベクトルpn-1の値に、第(n-4)時間ステップから第(n-1)時間ステップの微分ベクトルvn-4~vn-1により定められる第1の変化ベクトルの値を加算することにより、第n時間ステップでの音圧ベクトルpの値を計算する。第1の変化ベクトルは、第(n-1)時間ステップから第n時間ステップまでの間における音圧ベクトルpの変化量を表すベクトルである。第1の変化ベクトルは、具体的には、(14)式に示したように、第(n-4)時間ステップから第(n-1)時間ステップの微分ベクトルvn-1~vn-4を、それぞれに重み係数b~bを乗じた上で加算し、更に質量行列Mと対角質量行列Dの逆行列と時間刻み幅Δtとを乗算することにより計算される。重み係数b~bは、(9)式に示した分散誤差の1次から3次までの係数F,F,Fを0にする値、言い換えると、時間の離散化に起因する誤差に関して4次精度を達成する値に最適化されている。
【0076】
微分ベクトル計算部125は、微分ベクトルvの時間更新式である(5)式に従って、第n時間ステップでの微分ベクトルvを計算する。ここで、第n時間ステップでの微分ベクトルvは、解析対象の空間に設定された複数の要素のそれぞれにおける、計算開始時間ステップからn番目の時間ステップでの音圧の時間一階微分に相当する値を要素として有するベクトルである。
【0077】
上述したように、(5)式は、第n時間ステップでの微分ベクトルvの値を、第(n-1)時間ステップでの微分ベクトルvn-1の値と、第n時間ステップでの音圧ベクトルpの値と、により定める式である。従って、微分ベクトル計算部125は、第n時間ステップでの微分ベクトルvの値を、第(n-1)時間ステップでの微分ベクトルvn-1の値と、第n時間ステップでの音圧ベクトルpの値と、から計算する。
【0078】
具体的に説明すると、微分ベクトル計算部125は、第(n-1)時間ステップでの微分ベクトルvn-1の値に、音圧ベクトル計算部123により計算された第n時間ステップでの音圧ベクトルpにより定められる第2の変化ベクトルの値を加算することにより、第n時間ステップでの微分ベクトルvの値を計算する。第2の変化ベクトルは、第(n-1)時間ステップから第n時間ステップまでの間における微分ベクトルvの変化量を表すベクトルである。第2の変化ベクトルは、具体的には、(5)式に示したように、第n時間ステップでの外力ベクトルf、第n時間ステップでの音圧ベクトルp、剛性行列K、対角質量行列Dの逆行列、及び時間刻み幅Δtにより計算される。
【0079】
繰り返し部127は、nの値を1ずつ増加させながら、音圧ベクトル計算部123により第n時間ステップでの音圧ベクトルpの値を計算する処理と、微分ベクトル計算部125により第n時間ステップでの微分ベクトルvの値を計算する処理と、を繰り返す。具体的に説明すると、繰り返し部127は、nの値を1から総ステップ数Nに達するまで変化させながら、n=kのときの音圧ベクトルp及び微分ベクトルvの値を、n=k-1のときに計算された音圧ベクトルpk-1の値と、n=k-1~k-4のときに計算された微分ベクトルvk-1~vk-4の値と、を用いて計算する。これにより、繰り返し部127は、第1時間ステップから第N時間ステップまでの複数の時間ステップにおける音圧ベクトルp及び微分ベクトルvの値を計算する。
【0080】
なお、n=1のときの音圧ベクトルpn-1と微分ベクトルvn-1~vn-4、n=2のときの微分ベクトルvn-2~vn-4、n=3のときの微分ベクトルvn-3,vn-4、及びn=4のときの微分ベクトルvn-4は定義されていない。そのため、これらのベクトルの値は、初期条件として0等の適当な値に予め設定しておく。その上で、繰り返し部127は、(14)式及び(5)式に従って音圧ベクトルp及び微分ベクトルvを計算する。
【0081】
このようにして、数値計算部120は、解析対象の空間における音圧の分布を示す音圧ベクトルpとその時間一階微分と等価な微分ベクトルvとを1ステップずつ計算することにより、音波の時間進行を解析する。数値計算部120は、制御部11が記憶部12等と協働することにより実現される。
【0082】
図1に戻って、出力部130は、数値計算部120により計算された複数の時間ステップでの音圧ベクトルpに基づく出力情報を出力する。複数の時間ステップでの音圧ベクトルpに基づく出力情報とは、数値計算部120による数値計算により得られた第1時間ステップから第N時間ステップまでの音圧ベクトルp~pに基づいて表現することが可能な情報である。出力部130は、出力情報として、音圧の時間変化と、音圧の空間分布と、のうちの少なくとも一方を表す情報を出力する。
【0083】
例えば、図4に示すように、出力部130は、出力情報として、解析対象の空間内の特定の位置における音圧の時間変化を示す情報を表示部14に表示出力する。出力部130は、数値計算部120により計算された第1時間ステップから第N時間ステップまでの音圧ベクトルp~pに含まれる音圧の値のうちの、ユーザにより指定された位置の音圧の値を時系列順に並べる。これにより、出力部130は、特定の位置における音圧の時間変化を表す画像を生成し、図4に示すように表示部14に表示出力する。
【0084】
別の例として、図5に示すように、出力部130は、出力情報として、特定の時間ステップでの解析対象の空間における音圧の空間分布を示す情報を表示部14に表示出力する。出力部130は、数値計算部120により計算された第1時間ステップから第N時間ステップまでの音圧ベクトルp~pに含まれる音圧の値のうちの、ユーザにより指定された時間ステップにおける音圧ベクトルの値を、空間分布を表すように並べる。これにより、出力部130は、特定の時間ステップでの音圧の空間分布を表す画像を生成し、図5に示すように表示部14に表示出力する。なお、図5では、例として1次元方向(縦方向、横方向等)の分布を表す画像を表示しているが、出力部130は、音圧の2次元又は3次元の分布を表す画像を表示しても良い。
【0085】
或いは、出力部130は、音圧の時間変化と空間分布とをどちらも出力しても良い。例えば、出力部130は、解析対象の空間内を音波がどのように伝わるかを可視化するために、音圧の空間分布の複数の時間ステップにおける時間変化を表す動画像を生成し、表示部14に表示出力しても良い。ユーザは、このように表示部14に表示された音圧の時間変化や空間分布を確認することで、解析対象の空間内で音波がどのように伝搬又は分布するかの情報を得ることができる。
【0086】
更には、時間変化や空間分布以外にも、出力部130は、出力情報として、解析対象の空間の音響特性を表す音響パラメータ(例えば残響時間)を出力しても良い。出力部130に出力情報としてどのような情報を出力させるかは、ユーザが操作部13を介して適宜設定することができる。出力部130は、制御部11が表示部14等と協働することにより実現される。
【0087】
以上のように構成される音場解析装置10によって実行される音場解析処理の流れについて、図6に示すフローチャートを参照して説明する。図6に示す音場解析処理は、音場解析装置10が音場解析処理を実行可能な状態において、ユーザからの指示入力に従って適宜開始される。
【0088】
音場解析処理を開始すると、制御部11は、設定受付部110として機能し、音場解析における条件及びパラメータの設定を受け付ける(ステップS11)。具体的に説明すると、制御部11は、解析対象の空間の次元及び形状、境界条件、要素の形状及びサイズ、時間刻み幅Δt、総ステップ数N、外力ベクトルf等の設定を、操作部13を介してユーザから受け付ける。
【0089】
設定を受け付けると、制御部11は、行列設定部121として機能し、質量行列M、剛性行列K、及び対角質量行列Dを設定する(ステップS12)。制御部11は、非特許文献3と同様に、数値積分点α及びαとして、要素の形状に依存しない値である“α=±√(4/3)”及び“α=±√(2/3)”を用いる。そして、制御部11は、数値積分点α及びαからそれぞれ要素質量行列M及び要素剛性行列Kを計算し、質量行列M及び剛性行列Kを設定する。
【0090】
質量行列M、剛性行列K、及び対角質量行列Dを設定すると、制御部11は、nの値を1に初期化する(ステップS13)。そして、制御部11は、nの値を増加させながら、音圧ベクトルp及び微分ベクトルvを1ステップ毎に計算する処理に移行する。
【0091】
第1に、制御部11は、音圧ベクトル計算部123として機能し、(14)式に従って音圧ベクトルpを計算する(ステップS14)。具体的に説明すると、制御部11は、第(n-1)時間ステップでの音圧ベクトルpn-1と、第(n-4)時間ステップから第(n-1)時間ステップでの微分ベクトルvn-4~vn-1と、を(14)式に代入することにより、第n時間ステップでの音圧ベクトルpを計算する。
【0092】
音圧ベクトルpを計算すると、制御部11は、第2に、微分ベクトル計算部125として機能し、(5)式に従って微分ベクトルvを計算する(ステップS15)。具体的に説明すると、制御部11は、第(n-1)時間ステップでの微分ベクトルvn-1と、第n時間ステップでの外力ベクトルfと、ステップS14で計算された第n時間ステップでの音圧ベクトルpと、を(5)式に代入することにより、第n時間ステップでの微分ベクトルvを計算する。
【0093】
微分ベクトルvを計算すると、制御部11は、nの値が総ステップ数Nよりも小さいか否かを判定する(ステップS16)。nの値が総ステップ数Nよりも小さい場合(ステップS16;YES)、制御部11は、nの値をインクリメントする(ステップS17)。そして、制御部11は、繰り返し部127として機能し、処理をステップS14に戻してステップS14~S16の処理を再び実行する。
【0094】
具体的に説明すると、制御部11は、インクリメント前のnの値で既に計算された音圧ベクトルpと微分ベクトルvを用いて、インクリメント後のnの値における音圧ベクトルpと微分ベクトルvを計算する。そして、制御部11は、nの値が総ステップ数Nに達したか否かを判定し、nの値が総ステップ数Nに達していない場合、nの値をインクリメントして再びステップS14~S16の処理を実行する。このようにして、制御部11は、nの値が総ステップ数Nに達するまで、nの値を1ずつ増加させながら音圧ベクトルpと微分ベクトルvを計算することで、解析対象の空間における音波の時間進行を計算する。
【0095】
最終的に、nの値が総ステップ数Nに達すると(ステップS16;NO)、制御部11は、出力部130として機能し、計算結果を出力する(ステップS18)。例えば、制御部11は、ステップS14~S16の処理を繰り返すことで計算された音圧ベクトルp~pにより表現される音波の時間波形や音圧分布を示す画像を生成する。そして、制御部11は、図4及び図5に示したように、音圧の時間変化や空間分布を示す画像を表示部14に表示する。
【0096】
以上説明したように、本実施形態に係る音場解析装置10は、解析対象の空間を離散化する要素の形状に依存しない数値積分点を用いて計算された質量行列M及び剛性行列Kを含む一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する装置であって、一階常微分方程式において時間の離散化誤差を低減するために最適化された線形多段法を適用することにより得られる音圧ベクトルpの時間更新式と、微分ベクトルvの時間一階微分を1次精度の差分近似で離散化することにより得られる微分ベクトルvの時間更新式と、に従って、複数の時間ステップでの音圧ベクトルpの値を計算する。このように、質量行列M及び剛性行列Kが要素の形状に依存しない数値積分点を用いて計算されたものであり、且つ、時間の離散化誤差を低減するために最適化された線形多段法を適用するため、本実施形態に係る音場解析装置10は、空間を離散化する要素の形状を限定せずに、時間及び空間の離散化誤差に関して高い精度を達成する。
【0097】
特に、時間領域有限要素法では、時間及び空間の離散化に起因する誤差である離散化誤差が生じる。そのため、離散化誤差を抑えて計算結果を妥当なものにするためには、解析対象となる音波の波長に比べて十分に細かいメッシュが必要となる。しかしながら、可聴域の音波を解析対象とする場合、非常に細かい計算メッシュが必要となり、演算負荷が大きくなる。従って、音響分野で計算力学的手法を活用するためには、比較的粗い離散化でも誤差を抑えることができ、且つ、妥当な計算が実施できる計算アルゴリズムを開発することで、計算効率を向上させることが重要となる。非特許文献2に開示された修正積分則は、時間及び空間の離散化誤差に関して4次精度を達成するため計算効率に優れるが、使用できる要素の形状が正方形又は立方体に限定されるという制約がある。
【0098】
この打開策として、非特許文献3に開示された手法では、任意の形状の要素を扱うことができるが、時間の精度が低いため、長時間の時間応答の計算が要求される解析では不利である。従って、要素の形状の制約は、時間領域有限要素法における汎用性を低下させる欠点となる。本実施形態に係る音場解析装置10は、要素の形状の制約を受けない時間領域有限要素法において、時間及び空間の離散化誤差を4次の精度で低減するために最適化された高次の時間積分法を適用する。非特許文献3に開示された従来の手法では、時間2次精度であるため、時間刻み幅Δtを1/2にすると誤差が1/4に低減されるのに対して、実施形態の手法では、時間4次精度であるため、時間刻み幅Δtを1/2にすると誤差が1/16に低減することができる。その結果、任意の形状の要素を用いることができ、且つ、時間及び空間の離散化誤差に関して高い精度を達成することが可能な計算アルゴリズムを構築することができる。
【0099】
(数値実験)
次に、上記実施形態に係る音場解析装置10によって音圧をどの程度精度良く推定できるかを数値実験により評価した。これにより、上記実施形態に係る音場解析装置10による精度改善の効果を検証した。
【0100】
<数値実験1>
図7に、第1の数値実験において解析対象の空間として使用したモデルを示す。第1の数値実験では、図7に示すように、4m×3mの剛壁境界で囲われた2次元のモデル内における音場を解析対象として、従来の手法と実施形態の手法とを用いて解析した。
【0101】
ここで、従来の手法として、非特許文献3に開示された離散化誤差を空間4次及び時間2次の精度で低減可能な手法を使用した。解析には、要素長h=0.01[m]の正方形形状の要素を使用した。また、音速をc=343.7[m/s]と設定し、空気密度ρ=1.205[kg/m]と設定した。
【0102】
音源点を(x,y)=(2.0,1.0)に設置し、外力ベクトルfとして、時刻0で音源点に上限周波数4kHzのガウシアンパルスを発生させた。そして、受音点を(x,y)=(3.0,2.0)に設置し、受音点における時間応答を従来の手法と実施形態の手法で計算し、それぞれ参照解と比較した。
【0103】
参照解は、時間積分法として、非特許文献2に開示されたFox-Goodwin法を用いた陰的時間領域有限要素法により計算した。なお、参照解の計算においても、従来及び実施形態の手法と同様の有限要素を用いた。時間刻み幅Δtは、クーラン数τ(=cΔt/h)を0.75から0.01刻みで小さくしていき、安定計算が可能であった最大の値を採用した。具体的に、採用したクーラン数は、従来の手法でτ=0.64、実施形態の手法でτ=0.45、一般的なAdams法でτ=0.33であった。
【0104】
図8(a)に、従来の手法により計算された音圧の時間波形を、参照解の時間波形と比較した結果を示す。また、図8(b)に、実施形態の手法により計算された音圧の時間波形を、参照解の時間波形と比較した結果を示す。図8(a)において、濃い実線は、参照解(Reference)の時間波形を表し、薄い破線は、従来の手法(Conventional)により計算された音圧の時間波形を表す。また、図8(b)において、濃い実線は、参照解(Reference)の時間波形を表し、薄い破線は、実施形態の手法(Proposed)により計算された音圧の時間波形を表す。
【0105】
図8(a)では、従来の手法による音圧の時間波形は、境界で反射された最初の反射音が観測された約0.01sの時間ステップから参照解との乖離が見られる。そして、時間と共に誤差が蓄積していくのが確認される。一方で、図8(b)では、実施形態の手法による音圧の時間波形は、参照解と良好に一致している。そのため、実施形態の手法、すなわち高次の時間積分法の適用により、音圧の時間波形を高精度に計算可能であることが確認できた。
【0106】
更に、図9に、一般的なAdams法により計算された音圧の時間波形を、参照解の時間波形と比較した結果を示す。図9において、濃い実線は、参照解(Reference)の時間波形を表し、薄い破線は、一般的なAdams法(Classical)により計算された音圧の時間波形を表す。
【0107】
ここで、一般的なAdams法は、例えば“Philipp O.J. Scherer, Computational Physics -Simulation of Classical and Quantum Systems-, Springer, ed. 3(2017)”に開示された時間3次精度のAdams法を用いて音圧分布の時間進行を計算する手法である。具体的には、一般的なAdams法では、実施形態で説明した(14)式の代わりに、下記の(15)式を用いて音圧ベクトルpを計算する。
【0108】
【数15】
【0109】
図9に示すように、一般的なAdams法を用いた手法は、誤差が非常に大きく、ほとんど参照解を再現できていないことが分かる。これにより、時間の離散化誤差を低減するために最適化された線形多段法で音圧分布の時間進行を計算する実施形態の手法が有効であることを確認できた。
【0110】
以上のような第1の数値実験により、剛壁境界で囲われた音場において正方形形状の要素を使用した場合、実施形態の手法が従来の手法や一般的なAdams法を用いた場合よりも高い精度を持つことが確認された。
【0111】
<数値実験2>
次に、第2の数値実験として、実施形態の手法の、任意の形状の有限要素を用いた解析への適用性を検証するため、拡散体が設置された不整形な音場の解析を行った。第2の数値実験では、解析対象の空間として、図10に示すような10.8m×10.2mの剛壁境界に囲われた2次元の長方形状の音場を使用した。このような長方形状の音場内において、4つの境界のうちのx=10.8mの境界の近傍に、それぞれ音波を全反射する反射板によって構成される複数の拡散体を設置した。
【0112】
図10の下部に、音場内における破線で囲った部分の拡大図を示す。図10の下部に示すように、0.13m~0.27mの長さの拡散体を、x=10.8mの境界に対して15°~35°の角度で設置した。そして、このような長さ及び角度で配置された複数の拡散体のセットを、x=10.8mの境界に沿って、繰り返し設置した。
【0113】
このように様々な長さの拡散体が様々な角度で設置されているため、矩形形状の有限要素のみを用いては、音場を詳細に離散化することができない。そこで、矩形形状ではなく、最大要素長が0.025m以下となるように設定された四節点4辺形の要素を用いて、音場を離散化した。この音場の中心に、外力ベクトルfとして上限周波数3kHzのガウシアンパルスを与えた。なお、時間刻み幅は、Δt=1/150000[s]と設定した。
【0114】
図11に、計算開始から0.1sの間に位置(x,y)=(7.8,6.1)において観測された音圧の時間波形を、参照解と実施形態の手法とで比較した結果を示す。図11において、濃い実線は、参照解(Reference)の時間波形を表し、薄い破線は、実施形態の手法(Proposed)により計算された音圧の時間波形を表す。
【0115】
図11に示すように、実施形態の手法による時間波形と参照解の波形とは、良好に一致していることが分かる。そのため、第2の数値実験により、実施形態の手法を用いることで、矩形形状以外の形状のメッシュで離散化された音場を高精度で解析可能であることが確認された。
【0116】
(変形例)
以上に本発明の実施形態について説明したが、上記実施形態は一例であり、本発明の適用範囲はこれに限られない。すなわち、本発明の実施形態は種々の応用が可能であり、あらゆる実施の形態が本発明の範囲に含まれる。
【0117】
例えば、上記実施形態では、設定受付部110は、境界条件として、剛壁境界と振動境界とのいずれかの設定を受け付けた。しかしながら、本発明において、設定受付部110は、境界条件として、インピーダンス境界の設定を受け付けても良い。ここで、インピーダンス境界とは、音響インピーダンスを有する境界である。インピーダンス境界に入射した音波は、その境界の音響インピーダンス値zに応じた吸収率で吸収される。そのため、音波は、インピーダンス境界で反射することにより減衰する。
【0118】
解析対象の空間を囲む境界のうちの少なくとも一部にインピーダンス境界が含まれる場合、インピーダンス境界で反射された音波の減衰を評価するため、(2)式に示した一階常微分方程式には、減衰行列Cを含む減衰項が必要になる。そのため、音圧ベクトル計算部123は、(5)式の代わりに、減衰項が含まれる形式に変形された式に従って、微分ベクトルvを計算する。
【0119】
上記実施形態では、音場解析装置10は、表示部14を備えており、出力部130は、数値計算部120により計算された複数の時間ステップでの音圧ベクトルpに基づく出力情報を表示部14に表示出力した。しかしながら、本発明において、出力部130は、通信部15による通信を介して、出力情報を音場解析装置10の外部の装置に出力しても良い。そして、外部の装置が、例えば図4又は図5に示したように、出力部130から出力された出力情報を表示するようにしても良い。このように出力情報が外部の装置に出力される場合には、音場解析装置10は、表示部14を備えていなくても良い。
【0120】
上記実施形態では、音場解析装置10の制御部11において、CPUがROM又は記憶部12に記憶されたプログラムを実行することによって、設定受付部110、数値計算部120及び出力部130のそれぞれとして機能した。しかしながら、本発明において、制御部11は、専用のハードウェアであってもよい。専用のハードウェアとは、例えば単一回路、複合回路、プログラム化されたプロセッサ、ASIC(Application Specific Integrated Circuit)、FPGA(Field-Programmable Gate Array)、これらの組み合わせ等である。制御部11が専用のハードウェアである場合、各部の機能それぞれを個別のハードウェアで実現してもよいし、各部の機能をまとめて単一のハードウェアで実現してもよい。また、各部の機能のうち、一部を専用のハードウェアによって実現し、他の一部をソフトウェア又はファームウェアによって実現してもよい。このように、制御部11は、ハードウェア、ソフトウェア、ファームウェア、又は、これらの組み合わせによって、上述の各機能を実現することができる。
【0121】
本発明に係る音場解析装置の動作を規定する動作プログラムを既存のパーソナルコンピュータ、情報端末装置等に適用することで、当該パーソナルコンピュータ又は情報端末装置等を、本発明に係る音場解析装置として機能させることも可能である。
【0122】
また、このようなプログラムの配布方法は任意であり、例えば、CD-ROM(Compact Disk ROM)、DVD(Digital Versatile Disk)、MO(Magneto Optical Disk)、メモリカード等のコンピュータ読み取り可能な記録媒体に格納して配布してもよいし、インターネット等の通信ネットワークを介して配布してもよい。
【0123】
本発明は、本発明の広義の精神と範囲を逸脱することなく、様々な実施形態及び変形が可能とされるものである。また、上述した実施形態は、この発明を説明するためのものであり、本発明の範囲を限定するものではない。すなわち、本発明の範囲は、実施形態ではなく、特許請求の範囲によって示される。そして特許請求の範囲内及びそれと同等の発明の意義の範囲内で施される様々な変形が、この発明の範囲内とみなされる。
【符号の説明】
【0124】
10 音場解析装置
11 制御部
12 記憶部
13 操作部
14 表示部
15 通信部
20 空間
110 設定受付部
120 数値計算部
121 行列設定部
123 音圧ベクトル計算部
125 微分ベクトル計算部
127 繰り返し部
130 出力部
図1
図2
図3
図4
図5
図6
図7
図8
図9
図10
図11