JP7285513B2 - 音場解析装置、音場解析方法及びプログラム - Google Patents

Description

特許法第３０条第２項適用 発行者名：一般社団法人日本音響学会、刊行物名：日本音響学会２０１８年秋季研究発表会講演論文集、発行年月日：平成３０年８月２９日 集会名：日本音響学会２０１８年秋季研究発表会、開催日：平成３０年９月１２日

本発明は、音場解析装置、音場解析方法及びプログラムに関する。

数値シミュレーションにより音場を解析する手法が知られている。例えば、非特許文献１，２は、一階常微分方程式（ＯＤＥ：Ordinary Differential Equation）に基づく時間領域有限要素法（ＴＤＦＥＭ：Time-Domain Finite Element Method）を用いて音場を解析する手法を開示している。具体的に説明すると、非特許文献１は、以下の（手順１）～（手順４）に従って音場を解析する手法を開示している。

（手順１）
剛壁境界、振動境界又はインピーダンス境界により囲われた音場を対象とした場合の半離散化方程式に、対角質量行列と音圧の時間一階微分と等価なベクトルとを導入することで得られる連立一階常微分方程式を対象とする。連立一階常微分方程式は、具体的には、解析対象の空間における音圧の分布を示すベクトルｐ（以下、“音圧ベクトルｐ”と呼ぶ。）と、音圧ベクトルｐの時間一階微分と等価なベクトルｖ（以下、“微分ベクトルｖ”と呼ぶ。）と、を未知数として、下記（１）式及び（２）式のように表される。

ここで、（１）式及び（２）式における“Ｍ”、“Ｋ”、“Ｃ”、及び“Ｄ”は、それぞれ質量行列、剛性行列、減衰行列、及び対角質量行列を表す。質量行列Ｍ、剛性行列Ｋ、減衰行列Ｃ、及び対角質量行列Ｄは、それぞれ、要素質量行列Ｍ_ｅ、要素剛性行列Ｋ_ｅ、要素減衰行列Ｃ_ｅ、及び要素対角質量行列Ｄ_ｅを重ね合わせて構成される全体行列である。また、（１）式及び（２）式において、“ｆ”は解析対象の空間に加えられる外力の分布を示す外力ベクトルを表し、“ｃ_０”は音速を表し、“・”は時間に関する一階微分を表す。

（手順２）
空気領域の有限要素として正方形又は立方体形状の線形一次要素を設定する。そして、設定された空気領域の有限要素に対して、非特許文献２に開示された修正積分則を適用して、（１）式及び（２）式における質量行列Ｍと剛性行列Ｋとを計算する。非特許文献２に開示された修正積分則では、それぞれ要素質量行列Ｍ_ｅ及び要素剛性行列Ｋ_ｅをＧａｕｓｓ－Ｌｅｇｅｎｄｒｅ求積法により計算するための数値積分点α_ｍ及びα_ｋを、従来の値である±√（１／３）から下記（３）式のように修正する。これにより、時間及び空間の離散化を起因とする誤差である離散化誤差に関して４次精度を達成する。

（３）式において、“τ”はクーラン数を表す。クーラン数τは、音速ｃ_０と、時間領域有限要素法における時間刻み幅Δｔと、正方形又は立方体形状要素の要素長ｈと、を用いて、“τ＝ｃ_０Δｔ／ｈ”と定められる。

（手順３）
（１）式及び（２）式で表される連立一階常微分方程式を時間方向に離散化することで時間進行スキームを構築し、音波の時間進行を計算する。

（手順４）
（手順３）で音波の時間進行を計算する際に、数値的安定性を保つため、（１）式の音圧ベクトルｐの時間一階微分を１次精度の前進差分近似で離散化し、更に（２）式の微分ベクトルｖの時間一階微分を１次精度の後退差分近似で離散化する。この結果として、下記（４）式及び（５）式が得られる。

ここで、（４）式及び（５）式における“ｎ”は時間ステップ数を表す自然数であり、“Δｔ”は時間刻み幅を表す。非特許文献１に開示された手法では、（４）式及び（５）式により表される時間更新式に従って、音波の時間進行を計算する。なお、（５）式は、ｖ^ｎに関して陰的な式であるが、減衰行列Ｃに集中化を施すことによって、ｖ^ｎに関して陽的な式に変形することができる。

このような一階常微分方程式に基づく時間領域有限要素法は、計算効率に関する長所として次の２点を有する。
（Ｉ）時間進行スキームが陽的であるため、連立一次方程式の求解が不要である。そのため、時間ステップあたりの計算負荷が少ない。
（ＩＩ）修正積分則の適用により、使用する要素が一次要素であるにもかかわらず、時間及び空間で４次の精度を達成することができる。

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)

数値シミュレーションにより音場を解析する上で、計算効率の向上は重要である。非特許文献１に開示された一階常微分方程式に基づく時間領域有限要素法は、非特許文献２に開示された修正積分則の適用により、空気領域の計算に関して優れた計算効率を有する。その一方で、減衰項の時間離散化精度が低いため、境界上での吸音の予測精度が低いという課題がある。

本発明は、以上の課題を解決するためのものであり、時間領域有限要素法を用いた音場解析において、数値的安定性を悪化させることなく、吸音の予測精度を向上させることが可能な音場解析装置、音場解析方法及びプログラムを提供することを目的とする。

上記目的を達成するために、本発明の第１の観点に係る音場解析装置は、
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を囲む境界での音波の減衰を評価するための減衰項を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析装置であって、
前記一階常微分方程式において前記音圧ベクトルの時間一階微分を１次精度の差分近似で離散化することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において前記微分ベクトルの時間一階微分を１次精度の差分近似で離散化し、更に前記減衰項における前記音圧ベクトルの時間一階微分を２次精度の差分近似で離散化することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する数値計算部、を備え、
前記微分ベクトルの時間更新式は、第ｎ時間ステップでの前記微分ベクトルに関して陰的な式であり、
前記数値計算部は、連立一次方程式の解法として直接法または反復法を局所的に用いて前記微分ベクトルの時間更新式を解くことにより、前記第ｎ時間ステップでの前記微分ベクトルの値を計算する。

前記音圧ベクトルの時間更新式は、第ｎ時間ステップでの前記音圧ベクトルの値を、第（ｎ－１）時間ステップでの前記音圧ベクトルの値と、前記第ｎ時間ステップでの前記微分ベクトルの値と、により定める式であり、
前記微分ベクトルの時間更新式は、前記第ｎ時間ステップでの前記微分ベクトルの値を、第（ｎ－２）時間ステップ及び前記第（ｎ－１）時間ステップでの前記音圧ベクトルの値と、前記第（ｎ－１）時間ステップでの前記微分ベクトルの値と、により定める式であっても良い。

前記数値計算部は、
前記微分ベクトルの時間更新式に従って、前記第ｎ時間ステップでの前記微分ベクトルの値を、前記第（ｎ－２）時間ステップ及び前記第（ｎ－１）時間ステップでの前記音圧ベクトルの値と、前記第（ｎ－１）時間ステップでの前記微分ベクトルの値と、から計算する微分ベクトル計算部と、
前記音圧ベクトルの時間更新式に従って、前記第ｎ時間ステップでの前記音圧ベクトルの値を、前記第（ｎ－１）時間ステップでの前記音圧ベクトルの値と、前記第ｎ時間ステップでの前記微分ベクトルの値と、から計算する音圧ベクトル計算部と、
ｎの値を１ずつ増加させながら、前記微分ベクトル計算部により前記第ｎ時間ステップでの前記微分ベクトルの値を計算する処理と、前記音圧ベクトル計算部により前記第ｎ時間ステップでの前記音圧ベクトルの値を計算する処理と、を繰り返す繰り返し部と、を備えても良い。

前記境界のうちの少なくとも一部にはインピーダンス境界が含まれ、
前記数値計算部は、前記第ｎ時間ステップでの前記微分ベクトルの値のうちの、前記インピーダンス境界上の節点における値を前記直接法または前記反復法で計算し、前記第ｎ時間ステップでの前記微分ベクトルの値のうちの、前記インピーダンス境界上以外の節点における値を陽的に計算しても良い。

前記空間は、正方形又は立方体形状の一次要素で離散化され、
前記一階常微分方程式は、時間及び前記空間の離散化に起因する誤差に関して４次精度を達成する数値積分点を用いて計算された質量行列及び剛性行列を含んでも良い。

前記数値計算部により計算された前記複数の時間ステップでの前記音圧ベクトルの値に基づく出力情報を出力する出力部、を更に備えても良い。

前記出力部は、前記出力情報として、前記音圧の時間変化と、前記音圧の空間分布と、のうちの少なくとも一方を表す情報を出力しても良い。

上記目的を達成するために、本発明の第２の観点に係る音場解析方法は、
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を囲む境界での音波の減衰を評価するための減衰項を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析方法であって、
前記一階常微分方程式において前記音圧ベクトルの時間一階微分を１次精度の差分近似で離散化することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において前記微分ベクトルの時間一階微分を１次精度の差分近似で離散化し、更に前記減衰項における前記音圧ベクトルの時間一階微分を２次精度の差分近似で離散化することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算し、
前記微分ベクトルの時間更新式は、第ｎ時間ステップでの前記微分ベクトルに関して陰的な式であり、
前記第ｎ時間ステップでの前記微分ベクトルの値は、連立一次方程式の解法として直接法または反復法を局所的に用いて前記微分ベクトルの時間更新式を解くことにより、計算される。

上記目的を達成するために、本発明の第３の観点に係るプログラムは、
コンピュータを、
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を囲む境界での音波の減衰を評価するための減衰項を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析装置として機能させるプログラムであって、
前記コンピュータを、
前記一階常微分方程式において前記音圧ベクトルの時間一階微分を１次精度の差分近似で離散化することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において前記微分ベクトルの時間一階微分を１次精度の差分近似で離散化し、更に前記減衰項における前記音圧ベクトルの時間一階微分を２次精度の差分近似で離散化することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する数値計算部、として機能させ、
前記微分ベクトルの時間更新式は、第ｎ時間ステップでの前記微分ベクトルに関して陰的な式であり、
前記数値計算部は、連立一次方程式の解法として直接法または反復法を局所的に用いて前記微分ベクトルの時間更新式を解くことにより、前記第ｎ時間ステップでの前記微分ベクトルの値を計算する。

本発明によれば、時間領域有限要素法を用いた音場解析において、数値的安定性を悪化させることなく、吸音の予測精度を向上させることができる。

本発明の実施形態に係る音場解析装置の構成を示すブロック図である。 実施形態におけるメッシュの設定例を示す図である。 実施形態に係る音場解析装置における数値計算部の構成を示すブロック図である。 実施形態に係る音場解析装置によって得られた音圧の時間変化の表示例を示す図である。 実施形態に係る音場解析装置によって得られた音圧の空間分布の表示例を示す図である。 実施形態に係る音場解析装置によって実行される処理の流れを示すフローチャートである。 本発明の数値実験において解析対象として使用した３次元の音響管モデルを示す図である。 数値実験において使用されたメッシュの要素サイズ、空間解像度、及び時間刻み幅を示す図である。 数値実験において計算された複素音圧分布の理論解と数値解とを示す図である。 数値実験において、実施形態の手法で計算された数値解と理論解との間の相対誤差ノルムと、従来の手法で計算された数値解と理論解との間の相対誤差ノルムと、を比較して示す図である。

以下、本発明の実施形態について、図面を参照して説明する。なお、図中同一又は相当する部分には同一符号を付す。

（実施形態）
本実施形態に係る音場解析装置１０は、時間領域有限要素法を用いた波動音響数値シミュレーションにより音場を解析し、これにより音響分野における材料や測定法の開発、設計支援等に活用することが可能な装置である。ここで、音場とは、音波が存在する空間を意味する。

図１に、音場解析装置１０の構成を示す。図１に示すように、音場解析装置１０は、制御部１１と、記憶部１２と、操作部１３と、表示部１４と、通信部１５と、を備える。

制御部１１は、音場解析装置１０の各部と制御バスを介して接続されており、音場解析装置１０全体の動作を統括制御する。具体的に説明すると、制御部１１は、ＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）及びＲＡＭ（Random Access Memory）を備える。ＣＰＵは、例えばマイクロプロセッサ等であって、様々な処理及び演算を実行する中央演算処理部である。制御部１１において、ＣＰＵが、ＲＯＭに記憶されている制御プログラムを読み出して、ＲＡＭをワークメモリとして用いながら、各種の処理を実行する。

記憶部１２は、フラッシュメモリ、ハードディスク等の不揮発性メモリである。記憶部１２は、オペレーションシステム、アプリケーションプログラム等の、制御部１１が各種処理を行うために使用するプログラム及びデータを記憶する。また、記憶部１２は、制御部１１が各種処理を行うことにより生成又は取得するデータを記憶する。

操作部１３は、キーボード、マウス、ボタン、タッチパッド、タッチパネル等の入力デバイスを備えており、ユーザから操作を受け付ける。ユーザは、操作部１３を操作することによって、様々な指示を音場解析装置１０に入力することができる。操作部１３は、ユーザから入力された操作指示を受け付けると、受け付けた操作指示を制御部１１に送信する。

表示部１４は、液晶ディスプレイ、有機ＥＬ（Electro Luminescence）ディスプレイ等の表示デバイスを備える。表示部１４は、図示しない表示駆動回路によって駆動され、制御部１１による制御のもとで様々な画像を表示する。例えば、表示部１４は、数値シミュレーションにより音場を解析した結果を表す画像を表示する。

通信部１５は、音場解析装置１０が外部の機器と通信するための通信インタフェースを備える。通信部１５は、例えばインターネット等の広域ネットワーク、ＬＡＮ（Local Area Network）、ＵＳＢ（Universal Serial Bus）等の有線又は無線による通信を介して、外部の機器と通信する。

制御部１１は、図１に示すように、機能的に、設定受付部１１０と、数値計算部１２０と、出力部１３０と、を備える。制御部１１において、ＣＰＵがＲＯＭに記憶されたプログラムをＲＡＭに読み出して実行することにより、これら各部として機能する。

設定受付部１１０は、音場解析装置１０において音場を解析するための設定を受け付ける。ここで、設定は、数値計算部１２０が数値シミュレーションにより音圧の分布を計算するための条件やパラメータ等であって、具体的には解析対象の空間、境界条件、メッシュ等に関する情報である。設定受付部１１０は、ユーザが所望の設定を入力するための設定画面を表示部１４に表示する。ユーザは、表示部１４に表示された設定画面を見ながら、操作部１３を操作して、設定を入力する。設定受付部１１０は、ユーザにより入力された設定を受け付ける。

第１に、設定受付部１１０は、音場解析の目的物である解析対象の空間の設定を受け付ける。具体的に説明すると、設定受付部１１０は、操作部１３を介してユーザから入力された操作に従って、解析対象の空間の次元及び形状の設定を受け付ける。

例えば、長方形状の２次元平面内における音場を解析する場合、設定受付部１１０は、解析対象の空間として、長方形状の２次元の空間の設定を受け付ける。或いは、３次元の室内空間内における音場を解析する場合には、設定受付部１１０は、解析対象の空間として、その室内空間に対応する形状の３次元の空間の設定を受け付ける。また、３次元の音響材料の内部における音場を解析する場合には、設定受付部１１０は、解析対象の空間として、その音響材料に対応する形状の３次元の空間の設定を受け付ける。

第２に、設定受付部１１０は、境界条件の設定を受け付ける。境界条件は、解析対象の空間を囲む境界に関する条件である。設定受付部１１０は、操作部１３を介してユーザから入力された操作に従って、境界条件として、剛壁境界と振動境界とインピーダンス境界とのうちのいずれかの設定を受け付ける。

剛壁境界とは、音響的に剛である境界であって、境界上における粒子の速度が０である境界である。剛壁境界に入射した音波は、吸収されずに完全に反射される。これに対して、振動境界とは、境界上における粒子が振動可能な境界である。また、インピーダンス境界とは、音響インピーダンスを有する境界である。インピーダンス境界に入射した音波は、その境界の音響インピーダンス値ｚ_ｎに応じた吸収率で吸収される。そのため、音波は、インピーダンス境界で反射することにより減衰する。

解析対象の空間を囲む境界のうちの少なくとも一部にインピーダンス境界が含まれる場合、インピーダンス境界で反射された音波の減衰を評価するため、一階常微分方程式には、（２）式に示したように減衰行列Ｃを含む減衰項が必要になる。これに対して、解析対象の空間を囲む境界にインピーダンス境界が含まれない場合、（２）式における減衰行列Ｃを含む減衰項は不要になる。以下では、解析対象の空間を囲む境界の少なくとも一部にインピーダンス境界が含まれる場合を例にとって説明する。境界条件としてインピーダンス境界を設定する場合、設定受付部１１０は、そのインピーダンス境界における音響インピーダンス値ｚ_ｎの設定をユーザから受け付ける。

第３に、設定受付部１１０は、メッシュの設定を受け付ける。ここで、メッシュは、解析対象の空間を数値的に解析するために、解析対象の空間を四角形等の単純な形状をした複数の有限要素に分割（離散化）したものである。なお、有限要素は、単に“要素”とも呼ぶ。設定受付部１１０は、操作部１３を介してユーザから入力された操作に従って、メッシュを構成する複数の要素の形状及びサイズの設定を受け付ける。

図２に、２次元の解析対象の空間２０を四角形状の要素で分割した例を示す。解析対象の空間２０が２次元である場合、解析対象の空間２０は、例えば図２に示すように、１辺の長さｈの正方形状の要素によって、横方向にＸ個、縦方向にＹ個に分割される。すなわち、解析対象の空間２０は、合計（Ｘ×Ｙ）個の要素で離散化される。或いは、解析対象の空間が３次元である場合には、図示は省略するが、解析対象の空間は、立方体、直方体等の３次元形状の要素によって、横方向、縦方向及び高さ方向に分割される。

音場解析装置１０は、このように離散化された要素の単位で音圧の値を計算することにより、解析対象の空間における音圧の分布を解析する。要素のサイズを小さく設定するほど、空間の離散化に伴う誤差が小さくなるため解析精度が向上するが、計算量が増大する。

第４に、設定受付部１１０は、時間刻み幅Δｔ、総ステップ数Ｎ、及び外力ベクトルｆの設定を受け付ける。時間刻み幅Δｔは、解析対象の空間における音波の時間進行を計算するための時間ステップの単位である。時間刻み幅Δｔを小さく設定するほど、時間の離散化に伴う誤差が小さくなるため解析精度が向上するが、計算量が増大する。総ステップ数Ｎは、解析対象の空間における音波の時間進行を初期状態から何ステップ先まで計算するかを示す値である。数値計算部１２０は、Δｔ×Ｎの時間区間に亘って、音圧を計算する。外力ベクトルｆは、計算開始から計算終了までの時間区間において、解析対象の空間に加えられる外力の分布を示すベクトルである。

このように、設定受付部１１０は、時間領域有限要素法により音場を解析するための設定を、操作部１３を介してユーザから受け付ける。設定受付部１１０は、制御部１１が操作部１３等と協働することにより実現される。

図１に戻って、数値計算部１２０は、設定受付部１１０により受け付けられた設定のもとで、数値計算を実行する。具体的に説明すると、数値計算部１２０は、解析対象の空間における音圧の分布を示す音圧ベクトルｐと、音圧ベクトルｐの時間一階微分と等価な微分ベクトルｖと、に関する一階常微分方程式から得られる音圧ベクトルｐの時間更新式と微分ベクトルｖの時間更新式と、に従って、複数の時間ステップでの音圧ベクトルｐの値を計算する。一階常微分方程式は、（１）式及び（２）式に示したように、境界での音波の減衰を評価するための減衰項を含む式である。

ここで、数値計算部１２０による計算手法について説明する。本実施形態において、数値計算部１２０は、数値的安定性を悪化させることなく、吸音予測精度の改善を図るために、非特許文献１に開示された（４）式及び（５）式とは異なる時間更新式に従って、数値計算を実行する。具体的に説明すると、数値計算部１２０により使用される時間更新式は、背景技術で説明した（手順１）～（手順４）のうちの（手順１）～（手順３）は共通とし、（手順４）を変えることで得られる。

まず、（１）式の音圧ベクトルｐの時間一階微分を１次精度の後退差分近似で離散化する。これにより、（１）式を下記（６）式のように変形する。更に、（２）式の微分ベクトルｖの時間一階微分を１次精度の前進差分近似で離散化する。これにより、（２）式を下記（７）式のように変形する。

次に、（７）式の減衰項、すなわち減衰行列Ｃを含む項における音圧ベクトルｐ^ｎ－１の時間一階微分を、時間２次精度の中央差分近似により離散化する。そして、音圧ベクトルｐ^ｎに（６）式を代入すると、下記（８）式が得られる。

更に、（８）式により表される音圧ベクトルｐの時間一階微分を（７）式の減衰項に代入すると、以下の（９）式が得られる。

（５）式が、（２）式における音圧ベクトルｐの時間一階微分を１次精度の差分近似で離散化することにより計算されたものであるのに対して、（９）式は、（７）式における音圧ベクトルｐの時間一階微分を２次精度の差分近似で離散化することにより計算されたものである。これは、非特許文献１に開示された手法で吸音の予測精度が低くなる原因を、（２）式の減衰項における時間離散化の精度が低いためであると仮定したためである。そのため、本実施形態では、（７）式の減衰項における音圧ベクトルｐの時間一階微分を、１次精度ではなく、２次精度で有限差分近似して得られた（９）式を使用する。数値計算部１２０は、このようにして得られた音圧ベクトルｐの時間更新式である（６）式と微分ベクトルｖの時間更新式である（９）式と、に従って、複数の時間ステップでの音圧ベクトルｐの値を計算する。

図３に、数値計算部１２０のより詳細な構成を示す。図３に示すように、数値計算部１２０は、行列設定部１２１と、微分ベクトル計算部１２３と、音圧ベクトル計算部１２５と、繰り返し部１２７と、の機能を有する。

行列設定部１２１は、（９）式と（６）式とに含まれる質量行列Ｍ、剛性行列Ｋ、減衰行列Ｃ、及び対角質量行列Ｄを設定する。質量行列Ｍ、剛性行列Ｋ、減衰行列Ｃ、及び対角質量行列Ｄは、解析対象の空間を離散化する要素の数に対応するサイズの行列であって、各行列に含まれる要素の値は、解析対象の空間の次元及び形状と、離散化する要素の形状と、に依存して定められる。そのため、行列設定部１２１は、設定受付部１１０によりユーザから受け付けられた設定に応じて、質量行列Ｍ、剛性行列Ｋ、減衰行列Ｃ、及び対角質量行列Ｄを設定する。これにより、行列設定部１２１は、音圧ベクトルｐ及び微分ベクトルｖの時間更新式を確立する。

より詳細には、行列設定部１２１は、非特許文献２に開示された修正積分則を適用する。すなわち、行列設定部１２１は、（３）式に示した数値積分点α_ｍ及びα_ｋを用いて、質量行列Ｍ及び剛性行列Ｋを計算する。（３）式に示した数値積分点α_ｍ及びα_ｋは、解析対象の空間が正方形又は立方体形状の一次要素で離散化された場合において、一次要素の要素長ｈに依存して定められる数値積分点であって、時間及び空間の離散化に起因する誤差に関して４次精度を達成する。

具体的に説明すると、行列設定部１２１は、設定受付部１１０により設定された時間刻み幅Δｔ及び要素長ｈからクーラン数τを計算し、（３）式に従って数値積分点α_ｍ及びα_ｋを設定する。そして、行列設定部１２１は、設定した数値積分点α_ｍ及びα_ｋを用いてＧａｕｓｓ－Ｌｅｇｅｎｄｒｅ求積法により要素質量行列Ｍ_ｅ及び要素剛性行列Ｋ_ｅを計算し、質量行列Ｍ及び剛性行列Ｋを設定する。ここで、数値積分点α_ｍ及びα_ｋから要素質量行列Ｍ_ｅ及び要素剛性行列Ｋ_ｅを計算し、更に質量行列Ｍ及び剛性行列Ｋを得る手法は、非特許文献１，２等に開示されている周知の手法を用いることができる。また、行列設定部１２１は、非特許文献１，２等に開示されている周知の手法を用いて、対角行列である要素対角質量行列Ｄ_ｅを計算し、対角質量行列Ｄを設定する。

更に、行列設定部１２１は、インピーダンス境界において設定された音響インピーダンス値ｚ_ｎを用いて要素減衰行列Ｃ_ｅを計算し、減衰行列Ｃを設定する。音響インピーダンス値ｚ_ｎから要素減衰行列Ｃ_ｅを計算し、更に減衰行列Ｃを得る手法は、非特許文献１等に開示されている周知の手法を用いることができる。

微分ベクトル計算部１２３は、微分ベクトルｖの時間更新式である（９）式に従って、第ｎ時間ステップでの微分ベクトルｖ^ｎを計算する。ここで、第ｎ時間ステップでの微分ベクトルｖ^ｎは、解析対象の空間に設定された複数の要素のそれぞれにおける、計算開始時点からｎ番目の時間ステップでの音圧の時間一階微分に相当する値を要素として有するベクトルである。

上述したように、（９）式は、第ｎ時間ステップでの微分ベクトルｖ^ｎの値を、第（ｎ－２）時間ステップ及び第（ｎ－１）時間ステップでの音圧ベクトルｐ^ｎ－２，ｐ^ｎ－１の値と、第（ｎ－１）時間ステップでの微分ベクトルの値ｖ^ｎ－１と、により定める式である。従って、微分ベクトル計算部１２３は、第ｎ時間ステップでの微分ベクトルｖ^ｎの値を、第（ｎ－２）時間ステップ及び第（ｎ－１）時間ステップでの音圧ベクトルｐ^ｎ－２，ｐ^ｎ－１の値と、第（ｎ－１）時間ステップでの微分ベクトルｖ^ｎ－１の値と、から計算する。

具体的に説明すると、微分ベクトル計算部１２３は、（９）式の右辺において、第（ｎ－１）時間ステップでの微分ベクトルｖ^ｎ－１及び対角質量行列Ｄにより定められる第１項と、第ｎ時間ステップでの外力ベクトルｆ^ｎ、第（ｎ－１）時間ステップでの音圧ベクトルｐ^ｎ－１及び剛性行列Ｋにより定められる第２項と、第（ｎ－２）時間ステップでの音圧ベクトルｐ^ｎ－２と第（ｎ－１）時間ステップでの音圧ベクトルｐ^ｎ－１との差分及び減衰行列Ｃにより第３項と、をそれぞれ計算する。これにより、微分ベクトル計算部１２３は、（９）式の左辺に含まれる第ｎ時間ステップでの微分ベクトルｖ^ｎを計算する。

ここで、（９）式は、第ｎ時間ステップでの微分ベクトルｖ^ｎに関して陰的な式である。そのため、微分ベクトルｖ^ｎを得るには、連立一次方程式の求解が必要になる。数値計算による連立一次方程式の解法には、直接法と反復法とがある。直接法は、連立一次方程式を、式変形により直接的に解く手法である。これに対して、反復法は、誤差が予め定められた収束判定値まで収束するまで反復計算により解を計算する手法である。

連立一次方程式の求解を行う場合、陽的計算の場合と比べて、演算量が増大する。これを回避するため、本実施形態では、連立一次方程式の解法を局所的に適用することで（９）式を計算する。

より詳細に説明すると、微分ベクトルｖ^ｎの係数行列は、インピーダンス境界上の節点に対応する行のみに非対角成分を有し、インピーダンス境界上以外の節点に対応する行では非対角成分を有さない。従って、（９）式において、インピーダンス境界上の節点については連立一次方程式の求解を行い、インピーダンス境界上以外の節点については陽的に計算することで、全節点に対して連立一次方程式の求解を行う場合と比べて、大幅に演算量を削減できる。そのため、微分ベクトル計算部１２３は、第ｎ時間ステップでの微分ベクトルｖ^ｎの値のうちの、インピーダンス境界上の節点における値を直接法または反復法で計算する。これに対して、微分ベクトル計算部１２３は、第ｎ時間ステップでの微分ベクトルｖ^ｎの値のうちの、インピーダンス境界上以外の節点における値を陽的に計算する。言い換えると、微分ベクトル計算部１２３は、解析対象の空間に設定されたインピーダンス境界上の節点については、（９）式を用いて陰的に計算し、インピーダンス境界上以外の節点については、（９）式から減衰項を除いた以下の（１０）式を用いて陽的に計算する。

このように、連立一次方程式の解法として直接法または反復法を局所的に適用することで、微分ベクトルｖ^ｎを計算するのに必要な演算量を低減することができる。一般的に、インピーダンス境界上の節点数は全体の節点数に比べて非常に少ないので、結果として、実施形態に手法は陽的な計算と同程度の演算量で第ｎ時間ステップでの微分ベクトルｖ^ｎの値を計算することができるようになる。

音圧ベクトル計算部１２５は、音圧ベクトルｐの時間更新式である（６）式に従って、第ｎ時間ステップでの音圧ベクトルｐ^ｎを計算する。ここで、第ｎ時間ステップでの音圧ベクトルｐ^ｎは、解析対象の空間に設定された複数の要素のそれぞれにおける、計算開始時間ステップからｎ番目の時間ステップでの音圧の値を要素として有するベクトルである。

上述したように、（６）式は、第ｎ時間ステップでの音圧ベクトルｐ^ｎの値を、第（ｎ－１）時間ステップでの音圧ベクトルｐ^ｎ－１の値と、第ｎ時間ステップでの微分ベクトルｖ^ｎの値と、により定める式である。従って、音圧ベクトル計算部１２５は、第ｎ時間ステップでの音圧ベクトルｐ^ｎの値を、第（ｎ－１）時間ステップでの音圧ベクトルｐ^ｎ－１の値と、微分ベクトル計算部１２３により計算された第ｎ時間ステップでの微分ベクトルｖ^ｎの値と、から計算する。

具体的に説明すると、音圧ベクトル計算部１２５は、第（ｎ－１）時間ステップでの音圧ベクトルｐ^ｎ－１に、微分ベクトル計算部１２３により計算された第ｎ時間ステップでの微分ベクトルｖ^ｎにより定められる変化ベクトルを加算することにより、第ｎ時間ステップでの音圧ベクトルｐ^ｎを計算する。変化ベクトルは、第（ｎ－１）時間ステップから第ｎ時間ステップまでの間における音圧ベクトルｐの変化量を表すベクトルである。変化ベクトルは、具体的には、（６）式に示したように、微分ベクトルｖ^ｎに対して、質量行列Ｍと対角質量行列Ｄの逆行列と時間刻み幅Δｔとを乗算することにより得られる。

繰り返し部１２７は、ｎの値を１ずつ増加させながら、微分ベクトル計算部１２３により第ｎ時間ステップでの微分ベクトルｖ^ｎの値を計算する処理と、音圧ベクトル計算部１２５により第ｎ時間ステップでの音圧ベクトルｐ^ｎの値を計算する処理と、を繰り返す。具体的に説明すると、繰り返し部１２７は、ｎの値を１から総ステップ数Ｎに達するまで変化させながら、ｎ＝ｋのときの微分ベクトルｖ^ｋ及び音圧ベクトルｐ^ｋの値を、ｎ＝ｋ－１，ｋ－２のときに計算された音圧ベクトルｐ^ｋ－１，ｐ^ｋ－２の値と、ｎ＝ｋ－１のときに計算された微分ベクトルｖ^ｋ－１の値と、を用いて計算する。これにより、繰り返し部１２７は、第１時間ステップから第Ｎ時間ステップまでの複数の時間ステップにおける音圧ベクトルｐ^ｎ及び微分ベクトルｖ^ｎの値を計算する。

なお、ｎ＝１のときの音圧ベクトルｐ^ｎ－１，ｐ^ｎ－２と微分ベクトルｖ^ｎ－１、及びｎ＝２のときの音圧ベクトルｐ^ｎ－２は定義されていない。そのため、これらのベクトルの値は、初期条件として０等の適当な値に予め設定しておく。その上で、繰り返し部１２７は、（６）式及び（９）式に従って音圧ベクトルｐ^ｎ及び微分ベクトルｖ^ｎを計算する。

このようにして、数値計算部１２０は、解析対象の空間における音圧の分布を示す音圧ベクトルｐとその時間一階微分と等価な微分ベクトルｖとを１ステップずつ計算することにより、音波の時間進行を計算する。数値計算部１２０は、制御部１１が記憶部１２等と協働することにより実現される。

図１に戻って、出力部１３０は、数値計算部１２０により計算された複数の時間ステップでの音圧ベクトルｐに基づく出力情報を出力する。複数の時間ステップでの音圧ベクトルｐに基づく出力情報とは、数値計算部１２０による数値計算により得られた第１時間ステップから第Ｎ時間ステップまでの音圧ベクトルｐ^１～ｐ^Ｎに基づいて表現することが可能な情報である。出力部１３０は、出力情報として、音圧の時間変化と、音圧の空間分布と、のうちの少なくとも一方を表す情報を出力する。

例えば、図４に示すように、出力部１３０は、出力情報として、解析対象の空間内の特定の位置における音圧の時間変化を示す情報を表示部１４に表示出力する。出力部１３０は、数値計算部１２０により計算された第１時間ステップから第Ｎ時間ステップまでの音圧ベクトルｐ^１～ｐ^Ｎに含まれる音圧の値のうちの、ユーザにより指定された位置の音圧の値を時系列順に並べる。これにより、出力部１３０は、特定の位置における音圧の時間変化を表す画像を生成し、図４に示すように表示部１４に表示出力する。

別の例として、図５に示すように、出力部１３０は、出力情報として、特定の時間ステップでの解析対象の空間における音圧の空間分布を示す情報を表示部１４に表示出力する。出力部１３０は、数値計算部１２０により計算された第１時間ステップから第Ｎ時間ステップまでの音圧ベクトルｐ^１～ｐ^Ｎに含まれる音圧の値のうちの、ユーザにより指定された時間ステップにおける音圧ベクトルの値を、空間分布を表すように並べる。これにより、出力部１３０は、特定の時間ステップでの音圧の空間分布を表す画像を生成し、図５に示すように表示部１４に表示出力する。なお、図５では、例として１次元方向（縦方向、横方向等）の分布を表す画像を表示しているが、出力部１３０は、音圧の２次元又は３次元の分布を表す画像を表示しても良い。

或いは、出力部１３０は、音圧の時間変化と空間分布とをどちらも出力しても良い。例えば、出力部１３０は、解析対象の空間内を音波がどのように伝わるかを可視化するために、音圧の空間分布の複数の時間ステップにおける時間変化を表す動画像を生成し、表示部１４に表示出力しても良い。ユーザは、このように表示部１４に表示された音圧の時間変化や空間分布を確認することで、解析対象の空間内で音波がどのように伝搬又は分布するかの情報を得ることができる。

更には、時間変化や空間分布以外にも、出力部１３０は、出力情報として、解析対象の空間の音響特性を表す音響パラメータ（例えば残響時間）を出力しても良い。出力部１３０に出力情報としてどのような情報を出力させるかは、ユーザが操作部１３を介して適宜設定することができる。出力部１３０は、制御部１１が表示部１４等と協働することにより実現される。

以上のように構成される音場解析装置１０によって実行される音場解析処理の流れについて、図６に示すフローチャートを参照して説明する。図６に示す音場解析処理は、音場解析装置１０が音場解析処理を実行可能な状態において、ユーザからの指示入力に従って適宜開始される。

音場解析処理を開始すると、制御部１１は、設定受付部１１０として機能し、音場解析における条件及びパラメータの設定を受け付ける（ステップＳ１１）。具体的に説明すると、制御部１１は、解析対象の空間の次元及び形状、境界条件、要素の形状及びサイズ、時間刻み幅Δｔ、総ステップ数Ｎ、外力ベクトルｆ、インピーダンス境界における音響インピーダンス値ｚ_ｎ等の設定を、操作部１３を介してユーザから受け付ける。

設定を受け付けると、制御部１１は、行列設定部１２１として機能し、質量行列Ｍ、剛性行列Ｋ、減衰行列Ｃ、及び対角質量行列Ｄを設定する（ステップＳ１２）。具体的に説明すると、制御部１１は、ステップＳ１１で設定された時間刻み幅Δｔとメッシュの要素長ｈとによりクーラン数τを計算する。そして、制御部１１は、（３）式に示した数値積分点α_ｍから要素質量行列Ｍ_ｅを計算し、質量行列Ｍを設定する。また、制御部１１は、（３）式に示した数値積分点α_ｋから要素剛性行列Ｋ_ｅを計算し、剛性行列Ｋを設定する。更に、制御部１１は、ステップＳ１１で設定された音響インピーダンス値ｚ_ｎから要素減衰行列Ｃ_ｅを計算し、減衰行列Ｃを設定する。

質量行列Ｍ、剛性行列Ｋ、減衰行列Ｃ、及び対角質量行列Ｄを設定すると、制御部１１は、ｎの値を１に初期化する（ステップＳ１３）。そして、制御部１１は、ｎの値を増加させながら、音圧ベクトルｐ及び微分ベクトルｖを１ステップ毎に計算する処理に移行する。

第１に、制御部１１は、微分ベクトル計算部１２３として機能し、（９）式に従って微分ベクトルｖ^ｎを計算する（ステップＳ１４）。具体的に説明すると、制御部１１は、第（ｎ－２）時間ステップ及び第（ｎ－１）時間ステップでの音圧ベクトルｐ^ｎ－２，ｐ^ｎ－１と、第（ｎ－１）時間ステップでの微分ベクトルｖ^ｎ－１と、第（ｎ－１）時間ステップでの外力ベクトルｆ^ｎ－１と、を（９）式に代入することにより、第ｎ時間ステップでの微分ベクトルｖ^ｎを計算する。このとき、（９）式は、微分ベクトルｖ^ｎに関して陰的な式であるため、制御部１１は、反復法を局所的に用いて（９）式を解くことにより、第ｎ時間ステップでの微分ベクトルｖ^ｎを計算する。

微分ベクトルｖ^ｎを計算すると、制御部１１は、第２に、音圧ベクトル計算部１２５として機能し、（６）式に従って音圧ベクトルｐ^ｎを計算する（ステップＳ１５）。具体的に説明すると、制御部１１は、第（ｎ－１）時間ステップでの音圧ベクトルｐ^ｎ－１と、ステップＳ１４で計算された第ｎ時間ステップでの微分ベクトルｖ^ｎと、を（６）式に代入することにより、第ｎ時間ステップでの音圧ベクトルｐ^ｎを計算する。

音圧ベクトルｐ^ｎを計算すると、制御部１１は、ｎの値が総ステップ数Ｎよりも小さいか否かを判定する（ステップＳ１６）。ｎの値が総ステップ数Ｎよりも小さい場合（ステップＳ１６；ＹＥＳ）、制御部１１は、ｎの値をインクリメントする（ステップＳ１７）。そして、制御部１１は、繰り返し部１２７として機能し、処理をステップＳ１４に戻してステップＳ１４～Ｓ１６の処理を再び実行する。

具体的に説明すると、制御部１１は、インクリメント前のｎの値で既に計算された音圧ベクトルｐ^ｎと微分ベクトルｖ^ｎを用いて、インクリメント後のｎの値における音圧ベクトルｐ^ｎと微分ベクトルｖ^ｎを計算する。そして、制御部１１は、ｎの値が総ステップ数Ｎに達したか否かを判定し、ｎの値が総ステップ数Ｎに達していない場合、ｎの値をインクリメントして再びステップＳ１４～Ｓ１６の処理を実行する。このようにして、制御部１１は、ｎの値が総ステップ数Ｎに達するまで、ｎの値を１ずつ増加させながら音圧ベクトルｐ^ｎと微分ベクトルｖ^ｎを計算することで、解析対象の空間における音波の時間進行を計算する。

最終的に、ｎの値が総ステップ数Ｎに達すると（ステップＳ１６；ＮＯ）、制御部１１は、出力部１３０として機能し、計算結果を出力する（ステップＳ１８）。例えば、制御部１１は、ステップＳ１４～Ｓ１６の処理を繰り返すことで計算された音圧ベクトルｐ^１～ｐ^Ｎにより表現される音波の時間波形や音圧分布を示す画像を生成する。そして、制御部１１は、図４及び図５に示したように、音圧の時間変化や空間分布を示す画像を表示部１４に表示する。

以上説明したように、本実施形態に係る音場解析装置１０は、一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する装置であって、音圧ベクトルｐの時間一階微分を１次精度の差分近似で離散化することにより得られる音圧ベクトルｐの時間更新式と、微分ベクトルｖの時間一階微分を１次精度の差分近似で離散化し、更に減衰項における音圧ベクトルｐの時間一階微分を２次精度の差分近似で離散化することにより得られる微分ベクトルｖの時間更新式と、に従って、複数の時間ステップでの音圧ベクトルｐの値を計算する。このように、減衰項における音圧ベクトルｐの時間一階微分を２次精度の差分近似で離散化するため、本実施形態に係る音場解析装置１０は、吸音の予測精度を向上させることができる。

特に、時間領域有限要素法では、時間及び空間の離散化に起因する誤差である離散化誤差が生じる。そのため、離散化誤差を抑えて計算結果を妥当なものにするためには、解析対象となる音波の波長に比べて十分に細かいメッシュが必要となる。しかしながら、可聴域の音波を解析対象とする場合、非常に細かい計算メッシュが必要となり、演算負荷が大きくなる。従って、音響分野で計算力学的手法を活用するためには、比較的粗い離散化でも誤差を抑えることができ、且つ、妥当な計算が実施できる計算アルゴリズムを開発することで、計算効率を向上させることが重要となる。本実施形態に係る音場解析装置１０は、境界上での吸音の予測精度が低い原因が、減衰項に関する時間微分に用いる有限差分近似の精度にあると仮定し、減衰項における離散化精度を向上させる。これにより、インピーダンス境界上での吸音の予測精度を改善することができる。その結果、吸音を扱う場合であっても、離散化誤差を低減させることが可能な計算アルゴリズムを構築することができる。

（数値実験）
次に、上記実施形態に係る音場解析装置１０によって音圧をどの程度精度良く推定できるかを数値実験により評価した。これにより、上記実施形態に係る音場解析装置１０による精度改善の効果を検証した。

図７に、数値実験において解析対象の空間として使用したモデルを示す。数値実験では、図７に示すように長さが１ｍで断面が０．０５ｍ×０．０５ｍの細長状の３次元音響管モデルを用意した。そして、音響管モデル内の空間を解析対象として、その複素音圧分布を、従来の手法と実施形態の手法とを用いて計算した。

境界条件として、音響管モデルにおけるｘ＝０．０の面を振動境界と設定し、対向するｘ＝１．０の面をインピーダンス境界と設定した。また、音響管の側面を剛壁境界と設定した。振動境界には、外力ベクトルｆとして、上限周波数３ｋＨｚのガウシアンパルスを与えた。また、インピーダンス境界には、周波数に非依存に境界の垂直入射吸音率α_ｎを０．１～１．０の範囲において０．１刻みで変化させられるように、音響インピーダンス値ｚ_ｎを実数で与えた。

また、インピーダンス境界上の節点のみに対して反復法を局所的に適用して（９）式を計算した。使用する反復法はＣＲ法（共役残差法）とし、収束判定値は１０^－６と設定した。

このような音響管モデルにおいて、要素サイズが一様な複数の立方体形状の要素によって構成されるメッシュを設定した。具体的には、要素サイズが異なる３種類のメッシュをそれぞれ設定し、３種類のメッシュをそれぞれについて、実施形態で説明した手法に従って複素音圧分布を計算した。

図８に、音響管モデル内に設定された３種類のメッシュであるＭｅｓｈ１、Ｍｅｓｈ２及びＭｅｓｈ３の要素サイズ“ｈ”（単位：ｍ）、空間解像度“λ／ｈ”、及び時間刻み幅“Δｔ”（単位：秒）を示す。図８において、“λ”は、解析時における上限周波数における音波長を表す。図８に示すように、空間解像度λ／ｈは、Ｍｅｓｈ１からＭｅｓｈ２にかけて２倍となり、更にＭｅｓｈ２からＭｅｓｈ３にかけて２倍となるように設定した。

図９に、複素音圧分布の計算例として、インピーダンス境界の垂直入射吸音率がα_ｎ＝０．９である場合における複素音圧分布の理論解と数値解とを比較した結果を示す。図９において、横軸は図７に示した音響管モデル内に設定されたｘ方向における座標（単位：ｍ）を表し、縦軸は音圧（単位：Ｐａ）を表している。また、図９において、濃い実線と薄い実線は、それぞれ理論解の実部（Ｔｈｅｏｒｙ（Ｒｅａｌ））と虚部（Ｔｈｅｏｒｙ（Ｉｍａｇｉｎａｒｙ））を表しており、濃い破線と薄い一点鎖線は、それぞれ数値解の実部（Ｐｒｏｐｏｓｅｄ（Ｒｅａｌ））と虚部（Ｐｒｏｐｏｓｅｄ（Ｉｍａｇｉｎａｒｙ））を表している。

図９に示した理論解は、下記（１１）式に示す前進波と後退波との重ね合わせの式と、下記（１２）式に示す境界条件の式とを用いて計算される解である。なお、（１１）式における“ｋ”は波数を表し、“Ａ”と“Ａ’”はそれぞれ前進波と後退波の振幅を表す。また、（１２）式における“ω”は角振動数を、“ρ_０”は空気密度を、“ｖ_ｓ”は振動速度をそれぞれ表し、“∂／∂ｎ”は、境界面上の音圧ベクトルｐの、境界面に垂直な方向であるｘ方向への微分を表す。理論解の計算において、振動速度ｖ_ｓとして、数値解との比較のため、数値解析に用いたガウシアンパルスを離散フーリエ変換した値を与えた。

これに対して、図９に示した数値解は、３種類のメッシュのうちのＭｅｓｈ２を用いた場合に、実施形態において説明した（９）式及び（６）式に従って音圧ベクトルｐの時間進行を計算して得られた解である。

理論解と数値解とを比較した結果、図９に示すように、実部と虚部のどちらにおいても、理論解と数値解とで振幅及び位相が精度良く一致していることが分かる。すなわち、実施形態の手法を用いることで、理論解と比べて妥当な計算が可能であることが確認された。また、実施形態の手法は従来の手法と同一の時間刻み幅Δｔで計算が可能であり、高次の差分近似を適用したにもかかわらず、数値的安定性が保たれることが確認された。すなわち、従来の手法及び実施形態の手法は、同一条件の計算を行う場合の計算負荷が同程度である。

更に、従来の手法及び実施形態の手法による数値解の精度を定量的に評価するため、音圧の空間分布に関して、数値解と理論解との間の相対誤差ノルム“ｅ”を、下記（１３）式のように定義した。なお、（１３）式における“ｐ_{ｔｈｅоｒｙ}”と“ｐ_ｎ”は、それぞれ音圧の理論解と数値解である。

このように定義された相対誤差ノルムｅを、従来の手法と実施形態の手法とのそれぞれについて、０．５～３ｋＨｚの範囲において０．５ｋＨｚ刻みで計算した。ここで、従来の手法とは、非特許文献１に開示された（４）式及び（５）式に従って音圧ベクトルｐの時間進行を計算する手法である。

図１０に、垂直入射吸音率がα_ｎ＝０．２、０．５、０．９である場合のそれぞれについて、相対誤差ノルムｅを、従来の手法と実施形態の手法とで比較して示す。図１０において、四角のデータ点は、Ｍｅｓｈ１～Ｍｅｓｈ３のそれぞれについて、実施形態の手法（Ｐｒｏｐｏｓｅｄ）で計算された数値解の相対誤差ノルムｅを表している。また、円のデータ点は、Ｍｅｓｈ１～Ｍｅｓｈ３のそれぞれについて、従来の手法（Ｃｏｎｖｅｎｔｉｏｎａｌ）で計算された数値解の相対誤差ノルムｅを表している。

図１０から、吸音条件によらず、実施形態の手法は従来の手法に比べて誤差が小さく、より高精度であることが確認できる。特に、低周波数帯域において、両手法間の精度に大きな差が生じている。更には、全周波数帯域において、吸音が大きいほど、実施形態の手法が従来の手法よりも精度が優れることが確認される。すなわち、数値実験により、解析周波数及び吸音条件に依らず、実施形態の手法は従来の手法より高精度であり、特に周波数が低いほど、そして吸音率が高いほど、実施形態の手法と従来の手法との差は大きくなることが確認された。

また、図１０における垂直入射吸音率がα_ｎ＝０．９である場合の図からは、Ｍｅｓｈ２を用いた実施形態の手法は、０．５～３ｋＨｚの全範囲に亘って、Ｍｅｓｈ３を用いた従来の手法よりも誤差が小さく高精度であることが確認される。すなわち、吸音率が高い場合には、一波長当たりの要素数が９程度であるＭｅｓｈ２を用いた実施形態の手法の誤差は、一波長当たりの要素数がその倍の１８程度であるＭｅｓｈ３を用いた従来の手法の誤差よりも小さいことが分かった。３次元の時間領域有限要素法では、一波長当たりの要素数が１／２になると、自由度はその３乗の１／８となるため、主要演算量もそれに比例して１／８程度となる。従って、実施形態の手法と従来の手法は同等の数値的安定性をもつので、吸音率が高い場合には、実施形態の手法を用いることで、従来の手法の１／８程度の演算量で、且つ、従来の手法よりも小さい誤差で音圧分布を計算可能であるとの効果が得られる。

（変形例）
以上に本発明の実施形態について説明したが、上記実施形態は一例であり、本発明の適用範囲はこれに限られない。すなわち、本発明の実施形態は種々の応用が可能であり、あらゆる実施の形態が本発明の範囲に含まれる。

例えば、上記実施形態では、音場解析装置１０は、非特許文献２に開示された修正積分則を適用した。すなわち、上記実施形態では、行列設定部１２１は、（３）式に従って数値積分点α_ｍ及びα_ｋを計算し、質量行列Ｍ及び剛性行列Ｋを設定した。しかしながら、本発明において、行列設定部１２１は、上記実施形態以外の手法で、すなわち（３）式以外の数値積分点α_ｍ及びα_ｋを用いて質量行列Ｍ及び剛性行列Ｋを設定しても良い。

上記実施形態では、音場解析装置１０は、表示部１４を備えており、出力部１３０は、数値計算部１２０により計算された複数の時間ステップでの音圧ベクトルｐに基づく出力情報を表示部１４に表示出力した。しかしながら、本発明において、出力部１３０は、通信部１５による通信を介して、出力情報を音場解析装置１０の外部の装置に出力しても良い。そして、外部の装置において、例えば図４又は図５に示したように、出力部１３０から出力された出力情報が表示されるようにしても良い。このように出力情報が外部の装置に出力される場合には、音場解析装置１０は、表示部１４を備えていなくても良い。

上記実施形態では、音場解析装置１０の制御部１１において、ＣＰＵがＲＯＭ又は記憶部１２に記憶されたプログラムを実行することによって、設定受付部１１０、数値計算部１２０及び出力部１３０のそれぞれとして機能した。しかしながら、本発明において、制御部１１は、専用のハードウェアであってもよい。専用のハードウェアとは、例えば単一回路、複合回路、プログラム化されたプロセッサ、ＡＳＩＣ（Application Specific Integrated Circuit）、ＦＰＧＡ（Field-Programmable Gate Array）、これらの組み合わせ等である。制御部１１が専用のハードウェアである場合、各部の機能それぞれを個別のハードウェアで実現してもよいし、各部の機能をまとめて単一のハードウェアで実現してもよい。また、各部の機能のうち、一部を専用のハードウェアによって実現し、他の一部をソフトウェア又はファームウェアによって実現してもよい。このように、制御部１１は、ハードウェア、ソフトウェア、ファームウェア、又は、これらの組み合わせによって、上述の各機能を実現することができる。

本発明に係る音場解析装置の動作を規定する動作プログラムを既存のパーソナルコンピュータ、情報端末装置等に適用することで、当該パーソナルコンピュータ又は情報端末装置等を、本発明に係る音場解析装置として機能させることも可能である。

また、このようなプログラムの配布方法は任意であり、例えば、ＣＤ－ＲＯＭ（Compact Disk ROM）、ＤＶＤ（Digital Versatile Disk）、ＭＯ（Magneto Optical Disk）、メモリカード等のコンピュータ読み取り可能な記録媒体に格納して配布してもよいし、インターネット等の通信ネットワークを介して配布してもよい。

本発明は、本発明の広義の精神と範囲を逸脱することなく、様々な実施形態及び変形が可能とされるものである。また、上述した実施形態は、この発明を説明するためのものであり、本発明の範囲を限定するものではない。すなわち、本発明の範囲は、実施形態ではなく、特許請求の範囲によって示される。そして特許請求の範囲内及びそれと同等の発明の意義の範囲内で施される様々な変形が、この発明の範囲内とみなされる。

１０ 音場解析装置
１１ 制御部
１２ 記憶部
１３ 操作部
１４ 表示部
１５ 通信部
２０ 空間
１１０ 設定受付部
１２０ 数値計算部
１２１ 行列設定部
１２３ 微分ベクトル計算部
１２５ 音圧ベクトル計算部
１２７ 繰り返し部
１３０ 出力部

Claims (9)

1. 解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を囲む境界での音波の減衰を評価するための減衰項を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析装置であって、
   前記一階常微分方程式において前記音圧ベクトルの時間一階微分を１次精度の差分近似で離散化することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において前記微分ベクトルの時間一階微分を１次精度の差分近似で離散化し、更に前記減衰項における前記音圧ベクトルの時間一階微分を２次精度の差分近似で離散化することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する数値計算部、を備え、
   前記微分ベクトルの時間更新式は、第ｎ時間ステップでの前記微分ベクトルに関して陰的な式であり、
   前記数値計算部は、連立一次方程式の解法として直接法または反復法を局所的に用いて前記微分ベクトルの時間更新式を解くことにより、前記第ｎ時間ステップでの前記微分ベクトルの値を計算する、
   音場解析装置。
2. 前記音圧ベクトルの時間更新式は、第ｎ時間ステップでの前記音圧ベクトルの値を、第（ｎ－１）時間ステップでの前記音圧ベクトルの値と、前記第ｎ時間ステップでの前記微分ベクトルの値と、により定める式であり、
   前記微分ベクトルの時間更新式は、前記第ｎ時間ステップでの前記微分ベクトルの値を、第（ｎ－２）時間ステップ及び前記第（ｎ－１）時間ステップでの前記音圧ベクトルの値と、前記第（ｎ－１）時間ステップでの前記微分ベクトルの値と、により定める式である、
   請求項１に記載の音場解析装置。
3. 前記数値計算部は、
   前記微分ベクトルの時間更新式に従って、前記第ｎ時間ステップでの前記微分ベクトルの値を、前記第（ｎ－２）時間ステップ及び前記第（ｎ－１）時間ステップでの前記音圧ベクトルの値と、前記第（ｎ－１）時間ステップでの前記微分ベクトルの値と、から計算する微分ベクトル計算部と、
   前記音圧ベクトルの時間更新式に従って、前記第ｎ時間ステップでの前記音圧ベクトルの値を、前記第（ｎ－１）時間ステップでの前記音圧ベクトルの値と、前記第ｎ時間ステップでの前記微分ベクトルの値と、から計算する音圧ベクトル計算部と、
   ｎの値を１ずつ増加させながら、前記微分ベクトル計算部により前記第ｎ時間ステップでの前記微分ベクトルの値を計算する処理と、前記音圧ベクトル計算部により前記第ｎ時間ステップでの前記音圧ベクトルの値を計算する処理と、を繰り返す繰り返し部と、を備える、
   請求項２に記載の音場解析装置。
4. 前記境界のうちの少なくとも一部にはインピーダンス境界が含まれ、
   前記数値計算部は、前記第ｎ時間ステップでの前記微分ベクトルの値のうちの、前記インピーダンス境界上の節点における値を前記直接法または前記反復法で計算し、前記第ｎ時間ステップでの前記微分ベクトルの値のうちの、前記インピーダンス境界上以外の節点における値を直接法で陽的に計算する、
   請求項１から３のいずれか１項に記載の音場解析装置。
5. 前記空間は、正方形又は立方体形状の一次要素で離散化され、
   前記一階常微分方程式は、時間及び前記空間の離散化に起因する誤差に関して４次精度を達成する数値積分点を用いて計算された質量行列及び剛性行列を含む、
   請求項１から４のいずれか１項に記載の音場解析装置。
6. 前記数値計算部により計算された前記複数の時間ステップでの前記音圧ベクトルの値に基づく出力情報を出力する出力部、を更に備える、
   請求項１から４のいずれか１項に記載の音場解析装置。
7. 前記出力部は、前記出力情報として、前記音圧の時間変化と、前記音圧の空間分布と、のうちの少なくとも一方を表す情報を出力する、
   請求項６に記載の音場解析装置。
8. 解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を囲む境界での音波の減衰を評価するための減衰項を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析方法であって、
   前記一階常微分方程式において前記音圧ベクトルの時間一階微分を１次精度の差分近似で離散化することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において前記微分ベクトルの時間一階微分を１次精度の差分近似で離散化し、更に前記減衰項における前記音圧ベクトルの時間一階微分を２次精度の差分近似で離散化することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算し、
   前記微分ベクトルの時間更新式は、第ｎ時間ステップでの前記微分ベクトルに関して陰的な式であり、
   前記第ｎ時間ステップでの前記微分ベクトルの値は、連立一次方程式の解法として直接法または反復法を局所的に用いて前記微分ベクトルの時間更新式を解くことにより、計算される、
   音場解析方法。
9. コンピュータを、
   解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を囲む境界での音波の減衰を評価するための減衰項を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析装置として機能させるプログラムであって、
   前記コンピュータを、
   前記一階常微分方程式において前記音圧ベクトルの時間一階微分を１次精度の差分近似で離散化することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において前記微分ベクトルの時間一階微分を１次精度の差分近似で離散化し、更に前記減衰項における前記音圧ベクトルの時間一階微分を２次精度の差分近似で離散化することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する数値計算部、として機能させ、
   前記微分ベクトルの時間更新式は、第ｎ時間ステップでの前記微分ベクトルに関して陰的な式であり、
   前記数値計算部は、連立一次方程式の解法として直接法または反復法を局所的に用いて前記微分ベクトルの時間更新式を解くことにより、前記第ｎ時間ステップでの前記微分ベクトルの値を計算する、
   プログラム。

Images