JP2022112624A - Sound field analysis apparatus, sound field analysis method, and program - Google Patents
Sound field analysis apparatus, sound field analysis method, and program Download PDFInfo
- Publication number
- JP2022112624A JP2022112624A JP2021008481A JP2021008481A JP2022112624A JP 2022112624 A JP2022112624 A JP 2022112624A JP 2021008481 A JP2021008481 A JP 2021008481A JP 2021008481 A JP2021008481 A JP 2021008481A JP 2022112624 A JP2022112624 A JP 2022112624A
- Authority
- JP
- Japan
- Prior art keywords
- time
- sound pressure
- vector
- differential
- sound
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 64
- 239000013598 vector Substances 0.000 claims abstract description 208
- 238000000034 method Methods 0.000 claims abstract description 106
- 238000004364 calculation method Methods 0.000 claims abstract description 59
- 238000009826 distribution Methods 0.000 claims abstract description 26
- 230000006870 function Effects 0.000 claims description 53
- 230000010354 integration Effects 0.000 claims description 20
- 230000001419 dependent effect Effects 0.000 claims description 19
- 238000012892 rational function Methods 0.000 claims description 9
- 239000011159 matrix material Substances 0.000 description 76
- 238000010521 absorption reaction Methods 0.000 description 18
- 230000008569 process Effects 0.000 description 13
- 238000004891 communication Methods 0.000 description 9
- 238000013016 damping Methods 0.000 description 9
- 230000002123 temporal effect Effects 0.000 description 9
- 238000002474 experimental method Methods 0.000 description 8
- 238000004088 simulation Methods 0.000 description 8
- 238000010586 diagram Methods 0.000 description 7
- 239000006096 absorbing agent Substances 0.000 description 6
- 230000008859 change Effects 0.000 description 6
- 238000012545 processing Methods 0.000 description 6
- 238000007796 conventional method Methods 0.000 description 5
- 238000012986 modification Methods 0.000 description 4
- 230000004048 modification Effects 0.000 description 4
- 238000005401 electroluminescence Methods 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 239000011148 porous material Substances 0.000 description 2
- 230000002238 attenuated effect Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000004141 dimensional analysis Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000009472 formulation Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
Images
Landscapes
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
Description
本発明は、音場解析装置、音場解析方法及びプログラムに関する。 The present invention relates to a sound field analysis device, sound field analysis method and program.
数値シミュレーションにより音場を解析する手法が知られている。例えば、非特許文献1,2は、一階常微分方程式(ODE:Ordinary Differential Equation)に基づく時間領域有限要素法(TDFEM:Time-Domain Finite Element Method)を用いて音場を解析する手法を開示している。具体的に説明すると、非特許文献1は、以下の(手順1)~(手順4)に従って音場を解析する手法を開示している。
A method of analyzing a sound field by numerical simulation is known. For example,
(手順1)
剛壁境界、振動境界又はインピーダンス境界により囲われた音場を対象とした場合の半離散化方程式に、対角質量行列と音圧の時間一階微分と等価なベクトルとを導入することで得られる連立一階常微分方程式を対象とする。連立一階常微分方程式は、具体的には、解析対象の空間における音圧の分布を示すベクトルp(以下、“音圧ベクトルp”と呼ぶ。)と、音圧ベクトルpの時間一階微分と等価なベクトルv(以下、“微分ベクトルv”と呼ぶ。)と、を未知数として、下記(1)式及び(2)式のように表される。
(Step 1)
It is obtained by introducing the diagonal mass matrix and the vector equivalent to the time-first derivative of the sound pressure into the semi-discretized equation for the sound field surrounded by the rigid wall boundary, the vibration boundary, or the impedance boundary. We are interested in a system of first-order ordinary differential equations Specifically, the simultaneous first-order ordinary differential equations are a vector p (hereinafter referred to as "sound pressure vector p") indicating the distribution of sound pressure in the space to be analyzed, and the time first-order differential of the sound pressure vector p A vector v equivalent to (hereinafter referred to as a "differential vector v") is expressed as the following equations (1) and (2) as unknowns.
ここで、(1)式及び(2)式における“M”、“K”、“C”、及び“D”は、それぞれ質量行列、剛性行列、減衰行列、及び対角質量行列を表す。質量行列M、剛性行列K、減衰行列C、及び対角質量行列Dは、それぞれ、要素質量行列Me、要素剛性行列Ke、要素減衰行列Ce、及び要素対角質量行列Deを重ね合わせて構成される全体行列である。また、(1)式及び(2)式において、“f”は解析対象の空間に加えられる外力の分布を示す外力ベクトルを表し、“c0”は音速を表し、“・”は時間に関する一階微分を表す。また、ynは、周波数に非依存な比音響アドミタンス比を表す。ynは、比音響インピーダンス比の逆数である。 Here, "M", "K", "C", and "D" in formulas (1) and (2) represent mass matrix, stiffness matrix, damping matrix, and diagonal mass matrix, respectively. The mass matrix M, stiffness matrix K, damping matrix C, and diagonal mass matrix D superimpose the element mass matrix M e , the element stiffness matrix K e , the element damping matrix C e , and the element diagonal mass matrix D e , respectively. It is the overall matrix constructed together. Further, in the equations (1) and (2), “f” represents an external force vector indicating the distribution of the external force applied to the space to be analyzed, “c 0 ” represents the speed of sound, and “·” represents time-related represents the derivative. Also, y n represents the frequency-independent specific acoustic admittance ratio. y n is the reciprocal of the specific acoustic impedance ratio.
(手順2)
空気領域の有限要素として正方形又は立方体形状の線形一次要素を設定する。そして、設定された空気領域の有限要素に対して、非特許文献2に開示された離散化誤差の低減法である修正積分則を適用して、(1)式及び(2)式における質量行列Mと剛性行列Kとを計算する。
(Step 2)
A square or cubic linear linear element is set as the finite element of the air region. Then, the modified integral rule, which is a discretization error reduction method disclosed in Non-Patent Document 2, is applied to the set finite elements of the air region, and the mass matrix in equations (1) and (2) Compute M and the stiffness matrix K.
(手順3)
(1)式及び(2)式で表される連立一階常微分方程式を時間方向に離散化することで時間進行スキームを構築し、音波の時間進行を計算する。
(Step 3)
A time progression scheme is constructed by discretizing the simultaneous first-order ordinary differential equations represented by equations (1) and (2) in the time direction, and the time progression of sound waves is calculated.
(手順4)
(手順3)で音波の時間進行を計算する際に、数値的安定性を保つため、(1)式の音圧ベクトルpの時間一階微分を1次精度の前進差分近似で離散化し、更に(2)式の微分ベクトルvの時間一階微分を1次精度の後退差分近似で離散化する。この結果として、下記(3)式及び(4)式が得られる。
(Step 4)
In order to maintain numerical stability when calculating the time progression of sound waves in (procedure 3), the time first derivative of the sound pressure vector p in equation (1) is discretized by forward difference approximation with first-order accuracy, and (2) Discretize the first-order time differential of the differential vector v of the equation (2) by backward difference approximation with first-order accuracy. As a result, the following formulas (3) and (4) are obtained.
ここで、(3)式及び(4)式における“n”は時間ステップ数を表す自然数であり、“Δt”は時間刻み幅を表す。非特許文献1に開示された手法では、(3)式及び(4)式により表される時間更新式に従って、音波の時間進行を計算する。なお、(4)式は、vnに関して陰的な式であるが、減衰行列Cに集中化を施すことによって、vnに関して陽的な式に変形することができる。
Here, "n" in equations (3) and (4) is a natural number representing the number of time steps, and "Δt" represents the time step width. In the method disclosed in
このような一階常微分方程式に基づく時間領域有限要素法は、計算効率に関する長所として次の2点を有する。
(I)時間進行スキームが陽的であるため、連立一次方程式の求解が不要である。そのため、時間ステップあたりの計算負荷が少ない。
(II)修正積分則の適用により、使用する要素が一次要素であるにもかかわらず、時間及び空間で4次の精度を達成することができる。
The time-domain finite element method based on such a first-order ordinary differential equation has the following two advantages regarding computational efficiency.
(I) Since the time marching scheme is explicit, no simultaneous linear equation solution is required. Therefore, the computational load per time step is small.
(II) By application of the modified integral rule, fourth order accuracy in time and space can be achieved despite the fact that the elements used are first order elements.
陽的な時間領域有限要素法を含む計算力学的手法を用いた音響シミュレーションでは、境界面における吸音の取り扱い方法が予測結果を大きく左右する。特に、境界の吸音における周波数依存性を適切に取り扱うことは、実際の音響現象を模擬するために重要である。しかし、陽的な時間領域有限要素法による音響解析手法では、実装の容易さと計算負荷の少なさから、境界の吸音を周波数に非依存な表面インピーダンスモデルで考慮しており、吸音の周波数依存性を近似的に扱っている。これは、時間領域解析での高精度な吸音モデルの定式化には、周波数依存特性の考慮に伴い、計算負荷の大きい畳み込みが発生するためである。そのため、時間領域有限要素法を用いた音場解析において、周波数に依存するインピーダンス境界を効率的に扱うことが可能な手法が求められている。 In acoustic simulations using computational mechanics methods including explicit time-domain finite element methods, how to handle sound absorption at the interface greatly influences the prediction results. In particular, proper handling of frequency dependence in boundary sound absorption is important for simulating actual acoustic phenomena. However, the acoustic analysis method based on the explicit time-domain finite element method considers the sound absorption at the boundary with a frequency-independent surface impedance model because of its ease of implementation and low computational load. is treated as an approximation. This is because the formulation of a highly accurate sound absorption model in time domain analysis involves convolution with a large computational load due to consideration of frequency dependent characteristics. Therefore, in sound field analysis using the time-domain finite element method, a method that can efficiently handle frequency-dependent impedance boundaries is required.
本発明は、以上の課題を解決するためのものであり、時間領域有限要素法を用いた音場解析において、周波数に依存するインピーダンス境界を効率的に扱うことが可能な音場解析装置、音場解析方法及びプログラムを提供することを目的とする。 SUMMARY OF THE INVENTION The present invention is intended to solve the above problems. An object of the present invention is to provide a field analysis method and program.
上記目的を達成するために、本発明の第1の観点に係る音場解析装置は、
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を囲む境界での音波の減衰を評価するための減衰項を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析装置であって、
前記一階常微分方程式において、前記音圧ベクトルの時間一階微分を差分近似で離散化することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において、前記減衰項に含まれる周波数に依存する比音響アドミタンス比の逆フーリエ変換値と前記音圧ベクトルの時間一階微分との畳み込みをAuxiliary Differential Equations法で近似することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する数値計算部、を備える。
In order to achieve the above object, the sound field analysis device according to the first aspect of the present invention includes:
A first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in a space to be analyzed and a differential vector equivalent to the time first order differential of said sound pressure vector, wherein the sound wave at the boundary surrounding said space A sound field analysis device for analyzing a sound field using a time domain finite element method based on the first-order ordinary differential equation including an attenuation term for evaluating the attenuation of
In the first-order ordinary differential equation, the time update formula of the sound pressure vector obtained by discretizing the time first-order differential of the sound pressure vector by difference approximation, and in the first-order ordinary differential equation, the attenuation term A time update formula for the differential vector obtained by approximating the convolution of the included frequency-dependent inverse Fourier transform value of the specific acoustic admittance ratio with the time first derivative of the sound pressure vector using the Auxiliary Differential Equations method, according to , a numerical calculator for calculating values of the sound pressure vector at a plurality of time steps.
前記微分ベクトルの時間更新式は、前記比音響アドミタンス比を有理関数近似することにより得られる補助関数を含み、
前記数値計算部は、前記補助関数の時間発展を時間積分法により計算しても良い。
The differential vector time update formula includes an auxiliary function obtained by rational function approximation of the specific acoustic admittance ratio,
The numerical calculation unit may calculate the time evolution of the auxiliary function by a time integration method.
前記数値計算部は、前記時間積分法として4次精度の陽的ルンゲクッタ法を用いて、前記補助関数の時間発展を計算しても良い。 The numerical calculation unit may calculate the time evolution of the auxiliary function using an explicit Runge-Kutta method with fourth-order precision as the time integration method.
前記音圧ベクトルの時間更新式は、第n時間ステップでの前記音圧ベクトルの値を、第(n-1)時間ステップでの前記音圧ベクトルの値と、第(n-1)時間ステップでの前記微分ベクトルの値と、により定める式であり、
前記微分ベクトルの時間更新式は、第n時間ステップでの前記微分ベクトルの値を、第(n-1)時間ステップ及び第n時間ステップでの前記音圧ベクトルの値と、第(n-1)時間ステップでの前記微分ベクトルの値と、第n時間ステップでの前記補助関数の値と、により定める式であっても良い。
The time update formula for the sound pressure vector is obtained by combining the value of the sound pressure vector at the nth time step with the value of the sound pressure vector at the (n−1)th time step and the (n−1)th time step. and the value of the differential vector at
The time update formula for the differential vector is such that the value of the differential vector at the nth time step is combined with the value of the sound pressure vector at the (n−1)th time step and the nth time step, and the (n−1)th ) A formula determined by the value of the differential vector at the time step and the value of the auxiliary function at the n-th time step.
前記微分ベクトルの時間更新式は、前記微分ベクトルの時間一階微分を一次精度の後退差分近似で離散化し、更に前記音圧ベクトルの時間一階微分を二次精度の中央差分近似で離散化することにより得られても良い。 The time update formula for the differential vector discretizes the time first derivative of the differential vector by backward difference approximation with first-order accuracy, and further discretizes the time first derivative of the sound pressure vector by central difference approximation with second-order accuracy. It may be obtained by
上記目的を達成するために、本発明の第2の観点に係る音場解析方法は、
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を囲む境界での音波の減衰を評価するための減衰項を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析方法であって、
前記一階常微分方程式において、前記音圧ベクトルの時間一階微分を差分近似で離散化することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において、前記減衰項に含まれる周波数に依存する比音響アドミタンス比の逆フーリエ変換値と前記音圧ベクトルの時間一階微分との畳み込みをAuxiliary Differential Equations法で近似することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する。
In order to achieve the above object, a sound field analysis method according to a second aspect of the present invention comprises:
A first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in a space to be analyzed and a differential vector equivalent to the time first order differential of said sound pressure vector, wherein the sound wave at the boundary surrounding said space A sound field analysis method for analyzing a sound field using a time domain finite element method based on the first-order ordinary differential equation including an attenuation term for evaluating the attenuation of
In the first-order ordinary differential equation, the time update formula of the sound pressure vector obtained by discretizing the time first-order differential of the sound pressure vector by difference approximation, and in the first-order ordinary differential equation, the attenuation term A time update formula for the differential vector obtained by approximating the convolution of the included frequency-dependent inverse Fourier transform value of the specific acoustic admittance ratio with the time first derivative of the sound pressure vector using the Auxiliary Differential Equations method, according to , to calculate the values of the sound pressure vector at multiple time steps.
上記目的を達成するために、本発明の第3の観点に係るプログラムは、
コンピュータを、
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を囲む境界での音波の減衰を評価するための減衰項を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析装置として機能させるプログラムであって、
前記コンピュータを、
前記一階常微分方程式において、前記音圧ベクトルの時間一階微分を差分近似で離散化することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において、前記減衰項に含まれる周波数に依存する比音響アドミタンス比の逆フーリエ変換値と前記音圧ベクトルの時間一階微分との畳み込みをAuxiliary Differential Equations法で近似することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する数値計算部、として機能させる。
In order to achieve the above object, a program according to the third aspect of the present invention,
the computer,
A first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in a space to be analyzed and a differential vector equivalent to the time first order differential of said sound pressure vector, wherein the sound wave at the boundary surrounding said space A program that functions as a sound field analysis device that analyzes a sound field using a time domain finite element method based on the first-order ordinary differential equation including an attenuation term for evaluating the attenuation of
said computer,
In the first-order ordinary differential equation, the time update formula of the sound pressure vector obtained by discretizing the time first-order differential of the sound pressure vector by difference approximation, and in the first-order ordinary differential equation, the attenuation term A time update formula for the differential vector obtained by approximating the convolution of the included frequency-dependent inverse Fourier transform value of the specific acoustic admittance ratio with the time first derivative of the sound pressure vector using the Auxiliary Differential Equations method, according to , a numerical calculation unit that calculates the values of the sound pressure vector at a plurality of time steps.
本発明によれば、時間領域有限要素法を用いた音場解析において、周波数に依存するインピーダンス境界を効率的に扱うことができる。 According to the present invention, frequency-dependent impedance boundaries can be efficiently handled in sound field analysis using the time domain finite element method.
以下、本発明の実施形態について、図面を参照して説明する。なお、図中同一又は相当する部分には同一符号を付す。 BEST MODE FOR CARRYING OUT THE INVENTION Hereinafter, embodiments of the present invention will be described with reference to the drawings. The same reference numerals are given to the same or corresponding parts in the drawings.
(実施形態)
本実施形態に係る音場解析装置10は、時間領域有限要素法を用いた波動音響数値シミュレーションにより音場を解析し、これにより音響分野における材料や測定法の開発、設計支援等に活用することが可能な装置である。ここで、音場とは、音波が存在する空間を意味する。
(embodiment)
The sound
図1に、音場解析装置10の構成を示す。図1に示すように、音場解析装置10は、制御部11と、記憶部12と、操作部13と、表示部14と、通信部15と、を備える。
FIG. 1 shows the configuration of the sound
制御部11は、音場解析装置10の各部と制御バスを介して接続されており、音場解析装置10全体の動作を統括制御する。具体的に説明すると、制御部11は、CPU(Central Processing Unit)、ROM(Read Only Memory)及びRAM(Random Access Memory)を備える。CPUは、例えばマイクロプロセッサ等であって、様々な処理及び演算を実行する中央演算処理部である。制御部11において、CPUが、ROMに記憶されている制御プログラムを読み出して、RAMをワークメモリとして用いながら、各種の処理を実行する。
The
記憶部12は、フラッシュメモリ、ハードディスク等の不揮発性メモリである。記憶部12は、オペレーションシステム、アプリケーションプログラム等の、制御部11が各種処理を行うために使用するプログラム及びデータを記憶する。また、記憶部12は、制御部11が各種処理を行うことにより生成又は取得するデータを記憶する。
The
操作部13は、キーボード、マウス、ボタン、タッチパッド、タッチパネル等の入力デバイスを備えており、ユーザから操作を受け付ける。ユーザは、操作部13を操作することによって、様々な指示を音場解析装置10に入力することができる。操作部13は、ユーザから入力された操作指示を受け付けると、受け付けた操作指示を制御部11に送信する。
The
表示部14は、液晶ディスプレイ、有機EL(Electro Luminescence)ディスプレイ等の表示デバイスを備える。表示部14は、図示しない表示駆動回路によって駆動され、制御部11による制御のもとで様々な画像を表示する。例えば、表示部14は、数値シミュレーションにより音場を解析した結果を表す画像を表示する。
The
通信部15は、音場解析装置10が外部の機器と通信するための通信インタフェースを備える。通信部15は、例えばインターネット等の広域ネットワーク、LAN(Local Area Network)、USB(Universal Serial Bus)等の有線又は無線による通信を介して、外部の機器と通信する。
The
制御部11は、図1に示すように、機能的に、設定受付部110と、数値計算部120と、出力部130と、を備える。制御部11において、CPUがROMに記憶されたプログラムをRAMに読み出して実行することにより、これら各部として機能する。
The
設定受付部110は、音場解析装置10において音場を解析するための設定を受け付ける。ここで、設定は、数値計算部120が数値シミュレーションにより音圧の分布を計算するための条件やパラメータ等であって、具体的には解析対象の空間、境界条件、メッシュ等に関する情報である。設定受付部110は、ユーザが所望の設定を入力するための設定画面を表示部14に表示する。ユーザは、表示部14に表示された設定画面を見ながら、操作部13を操作して、設定を入力する。設定受付部110は、ユーザにより入力された設定を受け付ける。
The
第1に、設定受付部110は、音場解析の目的物である解析対象の空間の設定を受け付ける。具体的に説明すると、設定受付部110は、操作部13を介してユーザから入力された操作に従って、解析対象の空間の次元及び形状の設定を受け付ける。
First, the
例えば、長方形状の2次元平面内における音場を解析する場合、設定受付部110は、解析対象の空間として、長方形状の2次元の空間の設定を受け付ける。或いは、3次元の室内空間内における音場を解析する場合には、設定受付部110は、解析対象の空間として、その室内空間に対応する形状の3次元の空間の設定を受け付ける。
For example, when analyzing a sound field in a rectangular two-dimensional plane, the
第2に、設定受付部110は、境界条件の設定を受け付ける。境界条件は、解析対象の空間を囲む境界に関する条件である。設定受付部110は、操作部13を介してユーザから入力された操作に従って、境界条件として、剛壁境界と振動境界とインピーダンス境界とのうちのいずれかの設定を受け付ける。
Second, the
剛壁境界とは、音響的に剛である境界であって、境界上における粒子の速度が0である境界である。剛壁境界に入射した音波は、吸収されずに完全に反射される。これに対して、振動境界とは、境界上における粒子が振動可能な境界である。また、インピーダンス境界とは、音響インピーダンスを有する境界である。インピーダンス境界に入射した音波は、その境界の音響インピーダンス値znに応じた吸収率で吸収される。そのため、音波は、インピーダンス境界で反射することにより減衰する。 A rigid wall boundary is a boundary that is acoustically rigid and the velocity of particles on the boundary is zero. A sound wave incident on a rigid wall boundary is completely reflected without being absorbed. In contrast, an oscillating boundary is a boundary on which particles can oscillate. Also, an impedance boundary is a boundary having an acoustic impedance. A sound wave incident on an impedance boundary is absorbed with an absorption rate corresponding to the acoustic impedance value zn of the boundary. Therefore, sound waves are attenuated by reflection at impedance boundaries.
解析対象の空間を囲む境界のうちの少なくとも一部にインピーダンス境界が含まれる場合、インピーダンス境界で反射された音波の減衰を評価するため、一階常微分方程式には、(2)式に示したように減衰行列Cを含む減衰項が必要になる。これに対して、解析対象の空間を囲む境界にインピーダンス境界が含まれない場合、(2)式における減衰行列Cを含む減衰項は不要になる。以下では、解析対象の空間を囲む境界の少なくとも一部にインピーダンス境界が含まれる場合を例にとって説明する。境界条件としてインピーダンス境界を設定する場合、設定受付部110は、そのインピーダンス境界における音響インピーダンス値znの設定をユーザから受け付ける。
When an impedance boundary is included in at least a part of the boundaries surrounding the space to be analyzed, the first-order ordinary differential equation includes the following equation (2) in order to evaluate the attenuation of the sound wave reflected by the impedance boundary. A damping term containing the damping matrix C is required as follows. On the other hand, if the boundary surrounding the space to be analyzed does not include an impedance boundary, the attenuation term including the attenuation matrix C in equation (2) becomes unnecessary. In the following, an example will be described in which an impedance boundary is included in at least a part of the boundary surrounding the space to be analyzed. When setting an impedance boundary as a boundary condition, the
第3に、設定受付部110は、メッシュの設定を受け付ける。ここで、メッシュは、解析対象の空間を数値的に解析するために、解析対象の空間を四角形等の単純な形状をした複数の有限要素に分割(離散化)したものである。なお、有限要素は、単に“要素”とも呼ぶ。設定受付部110は、操作部13を介してユーザから入力された操作に従って、メッシュを構成する複数の要素の形状及びサイズの設定を受け付ける。
Third, the
図2に、2次元の解析対象の空間20を四角形状の要素で分割した例を示す。解析対象の空間20が2次元である場合、解析対象の空間20は、例えば図2に示すように、1辺の長さhの正方形状の要素によって、横方向にX個、縦方向にY個に分割される。すなわち、解析対象の空間20は、合計(X×Y)個の要素で離散化される。或いは、解析対象の空間が3次元である場合には、図示は省略するが、解析対象の空間は、立方体、直方体等の3次元形状の要素によって、横方向、縦方向及び高さ方向に分割される。
FIG. 2 shows an example in which a two-dimensional
音場解析装置10は、このように離散化された要素の単位で音圧の値を計算することにより、解析対象の空間における音圧の分布を解析する。要素のサイズを小さく設定するほど、空間の離散化に伴う誤差が小さくなるため解析精度が向上するが、計算量が増大する。
The sound
第4に、設定受付部110は、時間刻み幅Δt、総ステップ数Ns、及び外力ベクトルfの設定を受け付ける。時間刻み幅Δtは、解析対象の空間における音波の時間進行を計算するための時間ステップの単位である。時間刻み幅Δtを小さく設定するほど、時間の離散化に伴う誤差が小さくなるため解析精度が向上するが、計算量が増大する。総ステップ数Nsは、解析対象の空間における音波の時間進行を初期状態から何ステップ先まで計算するかを示す値である。数値計算部120は、Δt×Nsの時間区間に亘って、音圧を計算する。外力ベクトルfは、計算開始から計算終了までの時間区間において、解析対象の空間に加えられる外力の分布を示すベクトルである。
Fourth, the
このように、設定受付部110は、時間領域有限要素法により音場を解析するための設定を、操作部13を介してユーザから受け付ける。設定受付部110は、制御部11が操作部13等と協働することにより実現される。
In this manner, the
図1に戻って、数値計算部120は、設定受付部110により受け付けられた設定のもとで、数値計算を実行する。具体的に説明すると、数値計算部120は、解析対象の空間における音圧の分布を示す音圧ベクトルpと、音圧ベクトルpの時間一階微分と等価な微分ベクトルvと、に関する一階常微分方程式から得られる音圧ベクトルpの時間更新式と微分ベクトルvの時間更新式と、に従って、複数の時間ステップでの音圧ベクトルpの値を計算する。一階常微分方程式は、(1)式及び(2)式に示した式であって、境界での音波の減衰を評価するための減衰項を含む。
Returning to FIG. 1,
ここで、数値計算部120による計算手法について説明する。本実施形態において、数値計算部120は、音圧ベクトルpの時間更新式として、背景技術で説明した(3)式を用いる。(3)式は、背景技術で説明したように、(1)式において音圧ベクトルpの時間一階微分を1次精度の差分近似で離散化することにより得られる音圧ベクトルpの時間更新式である。一方で、数値計算部120は、境界での吸音の周波数依存性を取り扱うために、微分ベクトルvの時間更新式として、背景技術で説明した(4)式とは異なる時間更新式を用いる。
Here, a calculation method by the
以下、本実施形態における微分ベクトルvの時間更新式の導出について説明する。周波数に依存するインピーダンス境界を扱う場合、背景技術で説明した(2)式は、以下の(5)式のように書き換えられる。ここで、(5)式における“*”は、以下の(6)式で定められる畳み込みを表す。また、(6)式において、“^”付きのyは、周波数に依存する比音響アドミタンス比y(ω)の逆フーリエ変換値を表す。 Derivation of the time update formula for the differential vector v in this embodiment will be described below. When dealing with frequency-dependent impedance boundaries, equation (2) described in Background Art is rewritten as equation (5) below. Here, "*" in the equation (5) represents convolution defined by the following equation (6). Also, in the equation (6), y with "^" represents the inverse Fourier transform value of the frequency-dependent specific acoustic admittance ratio y(ω).
ここで、一般的に、時間領域有限要素法による音響解析は、計算時間ステップ数が大きくなる。よって、(5)式における畳み込みの計算コストは膨大となるため、畳み込みを伴う手法の実用的な使用は困難である。従って、時間領域有限要素法で周波数に依存するインピーダンス境界を扱うには、効率的な畳み込み手法が必要となる。 Here, in general, acoustic analysis by the time domain finite element method requires a large number of calculation time steps. Therefore, since the computational cost of convolution in equation (5) is enormous, it is difficult to practically use a method involving convolution. Therefore, efficient convolution techniques are needed to deal with frequency-dependent impedance boundaries in the time-domain finite element method.
そこで、本実施形態では、計算効率に優れた再帰的な畳み込み手法であるADE(Auxiliary Differential Equations)法により、時間領域有限要素法に周波数依存性を効率的に組み込む。ここで、ADE法の詳細は、例えば“D. Dragna, A generalized recursive convolution method for time-domain propagation in porous media, Journal of Acoustical Society of America, 1030-1042 (2015) ”に開示されている。 Therefore, in this embodiment, the ADE (Auxiliary Differential Equations) method, which is a recursive convolution method with excellent computational efficiency, is used to efficiently incorporate frequency dependence into the time domain finite element method. Details of the ADE method are disclosed, for example, in "D. Dragna, A generalized recursive convolution method for time-domain propagation in porous media, Journal of Acoustical Society of America, 1030-1042 (2015)".
ADE法を用いる場合、周波数に依存する比音響アドミタンス比y(ω)を、以下の(7)式のように有理関数近似する。(7)式において、ω,jは、それぞれ角周波数と虚数単位j2=-1である。また、nrp,ncpは、それぞれ有理関数近似に使用される実数の極λiの数及び複素共役な極の対αi±jβiの数である。y∞,Ai,Bi,Ciは、係数である。 When using the ADE method, the frequency-dependent specific acoustic admittance ratio y(ω) is approximated by a rational function as shown in Equation (7) below. (7), ω and j are the angular frequency and the imaginary unit j 2 =−1, respectively. Also, n rp , n cp are the number of real poles λ i and the number of complex conjugate pole pairs α i ±jβ i used for rational function approximation, respectively. y ∞ , A i , B i , C i are coefficients.
(7)式における極λi,αi,βi及び係数y∞,Ai,Bi,Ciは、全て実数であり、その値はデータフィッティングの手法を用いることで得られる。例えば、上述の極及び係数に適当な初期値を与えた(7)式を用いて複数の周波数におけるy(ω)を計算した計算値と、その理論値又は実験値とを、比較する。そして、(7)式を用いた計算値と理論値又は実験値とが高い精度で一致するように、データフィッティングの手法を用いて、極λi,αi,βi及び係数y∞,Ai,Bi,Ciを決定する。なお、因果性を満たすため、λi,αiは正の値となる。 The poles λ i , α i , β i and the coefficients y ∞ , A i , B i , C i in the equation (7) are all real numbers, and their values are obtained by using a data fitting method. For example, a calculated value obtained by calculating y(ω) at a plurality of frequencies using equation (7) in which appropriate initial values are given to the poles and coefficients described above is compared with the theoretical or experimental value. Then, using a data fitting technique, poles λ i , α i , β i and coefficients y ∞ , A Determine i , B i , C i . Note that λ i and α i are positive values in order to satisfy causality.
有理関数近似された比音響アドミタンス比y(ω)を用いることで、(5)式の畳み込みは、以下の(8)式のように、補助関数φi,ψi (1),ψi (2)を用いて近似できる。補助関数φi,ψi (1),ψi (2)は、アキュムレータと呼ばれる。 By using the rational function-approximated specific acoustic admittance ratio y(ω), the convolution of the equation (5) becomes the auxiliary functions φ i , ψ i (1) , ψ i ( 2) can be used for approximation. The auxiliary functions φ i , φ i (1) , φ i (2) are called accumulators.
ここで、(8)式における補助関数φi,ψi (1),ψi (2)の時間発展は、それぞれ以下の(9)、(10)、(11)式により表される一階常微分方程式により計算される。 Here, the time evolution of the auxiliary functions φ i , ψ i (1) and ψ i (2) in the equation (8) are represented by the following equations (9), (10) and (11) respectively. It is calculated by an ordinary differential equation.
数値計算部120は、補助関数φi,ψi
(1),ψi
(2)の時間発展を、時間積分法により計算する。時間積分法として、例えば、ルンゲクッタ法、オイラー法、クランク・ニコルソン法、ホイン法等の周知の手法を用いることができる。
The
なお、時間積分法の精度と安定性は、計算全体の精度と安定性に影響を及ぼす。そのため、効率的なシミュレーションのためには、適切な時間積分法を選択する必要がある。後述の数値実験で示すように、陽的な時間領域有限要素法では、4次精度の陽的ルンゲクッタ法を用いることで(8)式を高精度かつ安定的に計算できる。そのため、以下では、時間積分法として、4次精度の陽的ルンゲクッタ法を用いる場合を例にとって説明する。なお、時間領域有限要素法における4次精度の陽的ルンゲクッタ法の詳細は、例えば“P. O. J. Scherer, Computational Physics: Simulation of Classical and Quantum Systems, Third ed. (Springer Nature,Berlin/Heidelberg, 2017).”に開示されている。 The accuracy and stability of the time integration method affect the accuracy and stability of the entire calculation. Therefore, it is necessary to choose an appropriate time integration method for efficient simulation. As shown in the numerical experiments described later, in the explicit time-domain finite element method, the explicit Runge-Kutta method with fourth-order accuracy can be used to calculate the equation (8) with high accuracy and stability. Therefore, the case where the explicit Runge-Kutta method with fourth-order precision is used as the time integration method will be described below as an example. For details of the fourth-order precision explicit Runge-Kutta method in the time-domain finite element method, see, for example, "P.O.J. Scherer, Computational Physics: Simulation of Classical and Quantum Systems, Third ed. (Springer Nature, Berlin/Heidelberg, 2017)." disclosed in
続いて、周波数に依存するインピーダンス境界を扱う場合の陽的スキームを構築する。具体的には、(8)式を(5)式に代入することで、以下の(12)式が導出される。 We then construct an explicit scheme for dealing with frequency-dependent impedance boundaries. Specifically, the following equation (12) is derived by substituting equation (8) into equation (5).
(12)式を安定的に解くため、微分ベクトルvの時間一階微分を一次精度の後退差分近似で離散化する。更に、周波数に依存するインピーダンス境界を高精度に計算するため、音圧ベクトルpの時間一階微分を2次精度の中央差分近似で離散化する。その際、音圧ベクトルpの時間一階微分は、(3)式を踏まえて以下の(13)式のように離散化される。 In order to stably solve the equation (12), the time first derivative of the differential vector v is discretized by the backward difference approximation of the first order accuracy. Furthermore, in order to calculate the frequency-dependent impedance boundary with high accuracy, the first-order time derivative of the sound pressure vector p is discretized by a second-order accuracy central difference approximation. At that time, the first-order time differential of the sound pressure vector p is discretized as shown in the following equation (13) based on the equation (3).
更に、完全な陽的スキームとするため、(13)式の質量行列Mを集中化して対角質量行列Dに変換する。結果として、微分ベクトルvの時間更新式は、以下の(14)式のように表される。(14)式は、周波数に依存する比音響アドミタンス比y(ω)の逆フーリエ変換値と音圧ベクトルpの時間一階微分との畳み込みをADE法により近似することにより得られた微分ベクトルvの時間更新式である。 Furthermore, for a fully explicit scheme, the mass matrix M in Eq. (13) is centralized and transformed into a diagonal mass matrix D. As a result, the time update formula for the differential vector v is expressed as the following formula (14). Equation (14) is a differential vector v is the time update formula of
数値計算部120は、このようにして導出された微分ベクトルvの時間更新式である(14)式と、音圧ベクトルpの時間更新式である(3)式と、に従って、複数の時間ステップでの音圧ベクトルpの値を計算する。
図3に、数値計算部120のより詳細な構成を示す。図3に示すように、数値計算部120は、行列設定部121と、音圧ベクトル計算部123と、補助関数計算部124と、微分ベクトル計算部125と、繰り返し部127と、の機能を有する。
FIG. 3 shows a more detailed configuration of the
行列設定部121は、(3)式と(14)式とに含まれる質量行列M、剛性行列K、減衰行列C、及び対角質量行列Dを設定する。質量行列M、剛性行列K、減衰行列C、及び対角質量行列Dは、解析対象の空間を離散化する要素の数に対応するサイズの行列であって、各行列に含まれる要素の値は、解析対象の空間の次元及び形状と、離散化する要素の形状と、に依存して定められる。そのため、行列設定部121は、設定受付部110によりユーザから受け付けられた設定に応じて、質量行列M、剛性行列K、減衰行列C、及び対角質量行列Dを設定する。これにより、行列設定部121は、音圧ベクトルp及び微分ベクトルvの時間更新式を確立する。
The
より詳細には、行列設定部121は、非特許文献2に開示された修正積分則を適用して、質量行列M及び剛性行列Kを計算する。計算方法としては、例えば、非特許文献2に習い、行列設定部121において、設定受付部110により設定された時間刻み幅Δt及び要素長hからクーラン数τを計算し、数値積分点αm及びαkを設定する。そして、行列設定部121は、設定した数値積分点αm及びαkを用いてGauss-Legendre求積法により要素質量行列Me及び要素剛性行列Keを計算し、質量行列M及び剛性行列Kを設定する。
More specifically, the
ここで、数値積分点αm及びαkから要素質量行列Me及び要素剛性行列Keを計算し、更に質量行列M及び剛性行列Kを得る手法は、非特許文献1,2等に開示されている周知の手法を用いることができる。また、行列設定部121は、非特許文献1,2等に開示されている周知の手法を用いて、対角行列である要素対角質量行列Deを計算し、対角質量行列Dを設定する。
Here, the method of calculating the element mass matrix M e and the element stiffness matrix K e from the numerical integration points α m and α k and further obtaining the mass matrix M and the stiffness matrix K are disclosed in
更に、行列設定部121は、インピーダンス境界において設定された音響インピーダンス値znを用いて要素減衰行列Ceを計算し、減衰行列Cを設定する。音響インピーダンス値znから要素減衰行列Ceを計算し、更に減衰行列Cを得る手法は、非特許文献1等に開示されている周知の手法を用いることができる。
Furthermore, the
音圧ベクトル計算部123は、音圧ベクトルpの時間更新式である(3)式に従って、第n時間ステップでの音圧ベクトルpnを計算する。ここで、第n時間ステップでの音圧ベクトルpnは、解析対象の空間に設定された複数の要素のそれぞれにおける、計算開始時点からn番目の時間ステップでの音圧の値を要素として有するベクトルである。
The sound
上述したように、(3)式は、第n時間ステップでの音圧ベクトルpnの値を、第(n-1)時間ステップでの音圧ベクトルpn-1の値と、第(n-1)時間ステップの微分ベクトルvn-1の値と、により定める式である。従って、音圧ベクトル計算部123は、第n時間ステップでの音圧ベクトルpnの値を、第(n-1)時間ステップでの音圧ベクトルpn-1の値と、第(n-1)時間ステップの微分ベクトルvn-1の値と、から計算する。
As described above, the equation (3) expresses the value of the sound pressure vector p_n at the n-th time step as the value of the sound pressure vector p_n-1 at the (n-1)th time step and the value of the sound pressure vector p_n- 1 at the (n-1)th time step. -1) the value of the differential vector v n-1 of the time step; Therefore, the sound
補助関数計算部124は、補助関数φi,ψi
(1),ψi
(2)の時間更新式である(9)~(11)式に従って、補助関数φi,ψi
(1),ψi
(2)の時間発展を計算する。上述したように、微分ベクトルvの時間更新式である(14)式は、比音響アドミタンス比y(ω)を有理関数近似することにより得られる補助関数φi,ψi
(1),ψi
(2)を含んでいる。補助関数計算部124は、4次精度の陽的ルンゲクッタ法を用いて、補助関数φi,ψi
(1),ψi
(2)の時間発展を計算する。
The auxiliary
具体的に説明すると、補助関数計算部124は、第(n-1)時間ステップでの補助関数φi,ψi
(1),ψi
(2)の値と、第(n-1)時間ステップでの音圧ベクトルpn-1の時間一階微分と、4次精度の陽的ルンゲクッタ法とを用いて(9)~(11)式の常微分方程式を解くことにより、第n時間ステップでの補助関数φi,ψi
(1),ψi
(2)の値を計算する。第(n-1)時間ステップでの音圧ベクトルpn-1の時間一階微分としては第(n-1)時間ステップでの微分ベクトルvn-1を用いる。
Specifically, the
微分ベクトル計算部125は、微分ベクトルvの時間更新式である(14)式に従って、第n時間ステップでの微分ベクトルvnを計算する。ここで、第n時間ステップでの微分ベクトルvnは、解析対象の空間に設定された複数の要素のそれぞれにおける、計算開始時間ステップからn番目の時間ステップでの音圧の時間一階微分に相当する値を要素として有するベクトルである。
The
上述したように、(14)式は、第n時間ステップでの微分ベクトルvnの値を、第(n-1)時間ステップ及び第n時間ステップでの音圧ベクトルpn-1,pnの値と、第(n-1)時間ステップでの微分ベクトルvn-1の値と、第n時間ステップでの補助関数φi,ψi
(1),ψi
(2)の値と、により定める式である。従って、微分ベクトル計算部125は、第n時間ステップでの微分ベクトルvnの値を、音圧ベクトル計算部123により計算された音圧ベクトルpn-1,pnの値と、補助関数計算部124により計算された補助関数φi,ψi
(1),ψi
(2)の値と、第(n-1)時間ステップでの微分ベクトルvn-1の値と、から計算する。
As described above, the equation (14) expresses the value of the differential vector v n at the n-th time step as the sound pressure vectors p n-1 and p n at the (n-1)-th time step and the n-th time step. , the value of the differential vector v n-1 at the (n-1)th time step, the values of the auxiliary functions φ i , ψ i (1) , ψ i (2) at the n-th time step, It is a formula determined by Therefore, the
繰り返し部127は、nの値を1ずつ増加させながら、音圧ベクトル計算部123により第n時間ステップでの音圧ベクトルpnの値を計算する処理と、補助関数計算部124により第n時間ステップでの補助関数φi,ψi
(1),ψi
(2)の値を計算する処理と、微分ベクトル計算部125により第n時間ステップでの微分ベクトルvnの値を計算する処理と、を繰り返す。具体的に説明すると、繰り返し部127は、nの値を1から総ステップ数Nsに達するまで変化させながら、n=kのときの音圧ベクトルpk、補助関数φi,ψi
(1),ψi
(2)及び微分ベクトルvkの値を、n=k-1のときに計算された音圧ベクトルpk-1、補助関数φi,ψi
(1),ψi
(2)及び微分ベクトルvk-1の値を用いて計算する。これにより、繰り返し部127は、第1時間ステップから第Ns時間ステップまでの複数の時間ステップにおける音圧ベクトルp1~pNs及び微分ベクトルv1~vNsの値を計算する。
The
なお、n=0のときの音圧ベクトルpnと微分ベクトルvn、及び補助関数φi,ψi
(1),ψi
(2)の値は定義されていない。そのため、これらのベクトルの値は、初期条件として0等の適当な値に予め設定しておく。その上で、繰り返し部127は、(3)式及び(14)式に従って音圧ベクトルpn及び微分ベクトルvnを計算する。
Note that the values of the sound pressure vector p n and differential vector v n when n=0 and the auxiliary functions φ i , φ i (1) and φ i (2) are not defined. Therefore, the values of these vectors are preset to appropriate values such as 0 as initial conditions. Then, the
このようにして、数値計算部120は、解析対象の空間における音圧の分布を示す音圧ベクトルpとその時間一階微分と等価な微分ベクトルvとを1ステップずつ計算することにより、音波の時間進行を解析する。数値計算部120は、制御部11が記憶部12等と協働することにより実現される。
In this way, the
図1に戻って、出力部130は、数値計算部120により計算された複数の時間ステップでの音圧ベクトルpに基づく出力情報を出力する。複数の時間ステップでの音圧ベクトルpに基づく出力情報とは、数値計算部120による数値計算により得られた第1時間ステップから第Ns時間ステップまでの音圧ベクトルp1~pNsに基づいて表現することが可能な情報である。出力部130は、出力情報として、音圧の時間変化と、音圧の空間分布と、のうちの少なくとも一方を表す情報を出力する。
Returning to FIG. 1 , the
例えば、図4に示すように、出力部130は、出力情報として、解析対象の空間内の特定の位置における音圧の時間変化を示す情報を表示部14に表示出力する。出力部130は、数値計算部120により計算された第1時間ステップから第Ns時間ステップまでの音圧ベクトルp1~pNsに含まれる音圧の値のうちの、ユーザにより指定された位置の音圧の値を時系列順に並べる。これにより、出力部130は、特定の位置における音圧の時間変化を表す画像を生成し、図4に示すように表示部14に表示出力する。
For example, as shown in FIG. 4, the
別の例として、図5に示すように、出力部130は、出力情報として、特定の時間ステップでの解析対象の空間における音圧の空間分布を示す情報を表示部14に表示出力する。出力部130は、数値計算部120により計算された第1時間ステップから第Ns時間ステップまでの音圧ベクトルp1~pNsに含まれる音圧の値のうちの、ユーザにより指定された時間ステップにおける音圧ベクトルの値を、空間分布を表すように並べる。これにより、出力部130は、特定の時間ステップでの音圧の空間分布を表す画像を生成し、図5に示すように表示部14に表示出力する。なお、図5では、例として1次元方向(縦方向、横方向等)の分布を表す画像を表示しているが、出力部130は、音圧の2次元又は3次元の分布を表す画像を表示しても良い。
As another example, as shown in FIG. 5, the
或いは、出力部130は、音圧の時間変化と空間分布とをどちらも出力しても良い。例えば、出力部130は、解析対象の空間内を音波がどのように伝わるかを可視化するために、音圧の空間分布の複数の時間ステップにおける時間変化を表す動画像を生成し、表示部14に表示出力しても良い。ユーザは、このように表示部14に表示された音圧の時間変化や空間分布を確認することで、解析対象の空間内で音波がどのように伝搬又は分布するかの情報を得ることができる。
Alternatively, the
更には、時間変化や空間分布以外にも、出力部130は、出力情報として、解析対象の空間の音響特性を表す音響パラメータ(例えば残響時間)を出力しても良い。出力部130に出力情報としてどのような情報を出力させるかは、ユーザが操作部13を介して適宜設定することができる。出力部130は、制御部11が表示部14等と協働することにより実現される。
Furthermore, the
以上のように構成される音場解析装置10によって実行される音場解析処理の流れについて、図6に示すフローチャートを参照して説明する。図6に示す音場解析処理は、音場解析装置10が音場解析処理を実行可能な状態において、ユーザからの指示入力に従って適宜開始される。
The flow of sound field analysis processing executed by the sound
音場解析処理を開始すると、制御部11は、設定受付部110として機能し、音場解析における条件及びパラメータの設定を受け付ける(ステップS11)。具体的に説明すると、制御部11は、解析対象の空間の次元及び形状、境界条件、要素の形状及びサイズ、時間刻み幅Δt、総ステップ数Ns、外力ベクトルf、インピーダンス境界における音響インピーダンス値zn等の設定を、操作部13を介してユーザから受け付ける。
When the sound field analysis process is started, the
設定を受け付けると、制御部11は、行列設定部121として機能し、質量行列M、剛性行列K、減衰行列C、及び対角質量行列Dを設定する(ステップS12)。具体的に説明すると、制御部11は、ステップS11で設定された時間刻み幅Δtとメッシュの要素長hとから数値積分点αm及びαkを設定する。そして、制御部11は、数値積分点αm及びαkからそれぞれ要素質量行列Me及び要素剛性行列Keを計算し、質量行列M及び剛性行列Kを設定する。更に、制御部11は、ステップS11で設定された音響インピーダンス値znから要素減衰行列Ceを計算し、減衰行列Cを設定する。
Upon receiving the settings, the
行列を設定すると、制御部11は、nの値を1に初期化する(ステップS13)。そして、制御部11は、nの値を増加させながら、音圧ベクトルp及び微分ベクトルvを1ステップ毎に計算する処理に移行する。
After setting the matrix, the
第1に、制御部11は、音圧ベクトル計算部123として機能し、(3)式に従って音圧ベクトルpnを計算する(ステップS14)。具体的に説明すると、制御部11は、第(n-1)時間ステップでの音圧ベクトルpn-1と、第(n-1)時間ステップでの微分ベクトルvn-1と、を(3)式に代入することにより、第n時間ステップでの音圧ベクトルpnを計算する。
First, the
音圧ベクトルpnを計算すると、制御部11は、第2に、補助関数計算部124として機能し、第n時間ステップでの補助関数φi,ψi
(1),ψi
(2)を計算する(ステップS15)。具体的に説明すると、制御部11は、第(n-1)時間ステップでの補助関数φi,ψi
(1),ψi
(2)の値と、第(n-1)時間ステップでの音圧ベクトルpn-1の時間一階微分と、4次精度の陽的ルンゲクッタ法とを用いて(9)~(11)式の常微分方程式を解くことにより第n時間ステップでの補助関数φi,ψi
(1),ψi
(2)の値を計算する。第(n-1)時間ステップでの音圧ベクトルpn-1の時間一階微分としては第(n-1)時間ステップでの微分ベクトルvn-1を用いる。
After calculating the sound pressure vector p n , the
補助関数φi,ψi
(1),ψi
(2)を計算すると、制御部11は、第3に、微分ベクトル計算部125として機能し、(14)式に従って微分ベクトルvnを計算する(ステップS16)。具体的に説明すると、制御部11は、第(n-1)時間ステップでの微分ベクトルvn-1と、第n時間ステップでの外力ベクトルfnと、ステップS14で計算された第n時間ステップでの音圧ベクトルpnと、ステップS15で計算された第n時間ステップでの補助関数φi,ψi
(1),ψi
(2)とを(14)式に代入することにより、第n時間ステップでの微分ベクトルvnを計算する。
After calculating the auxiliary functions φ i , φ i (1) and φ i (2) , the
微分ベクトルvnを計算すると、制御部11は、nの値が総ステップ数Nsよりも小さいか否かを判定する(ステップS17)。nの値が総ステップ数Nsよりも小さい場合(ステップS17;YES)、制御部11は、nの値をインクリメントする(ステップS18)。そして、制御部11は、繰り返し部127として機能し、処理をステップS14に戻してステップS14~S17の処理を再び実行する。
After calculating the differential vector vn, the
具体的に説明すると、制御部11は、インクリメント前のnの値で既に計算された音圧ベクトルpnと補助関数φi,ψi
(1),ψi
(2)と微分ベクトルvnを用いて、インクリメント後のnの値における音圧ベクトルpnと補助関数φi,ψi
(1),ψi
(2)と微分ベクトルvnを計算する。そして、制御部11は、nの値が総ステップ数Nsに達したか否かを判定し、nの値が総ステップ数Nsに達していない場合、nの値をインクリメントして再びステップS14~S17の処理を実行する。このようにして、制御部11は、nの値が総ステップ数Nsに達するまで、nの値を1ずつ増加させながら音圧ベクトルpnと微分ベクトルvnを計算することで、解析対象の空間における音波の時間進行を計算する。
Specifically, the
最終的に、nの値が総ステップ数Nsに達すると(ステップS17;NO)、制御部11は、出力部130として機能し、計算結果を出力する(ステップS19)。例えば、制御部11は、ステップS14~S17の処理を繰り返すことで計算された音圧ベクトルp1~pNsにより表現される音波の時間波形や音圧分布を示す画像を生成する。そして、制御部11は、図4及び図5に示したように、音圧の時間変化や空間分布を示す画像を表示部14に表示する。
Finally, when the value of n reaches the total number of steps Ns (step S17; NO), the
以上説明したように、本実施形態に係る音場解析装置10は、一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する装置であって、音圧ベクトルpの時間一階微分を差分近似で離散化することにより得られる音圧ベクトルpの時間更新式と、周波数に依存する比音響アドミタンス比y(ω)の逆フーリエ変換値と音圧ベクトルpの時間一階微分との畳み込みをADE法で近似することにより得られる微分ベクトルvの時間更新式と、に従って、複数の時間ステップでの音圧ベクトルpの値を計算する。このように、計算負荷の大きい畳み込みをADE法で近似するため、本実施形態に係る音場解析装置10は、周波数に依存するインピーダンス境界を効率的に扱うことができる。
As described above, the sound
特に、周波数に非依存な表面インピーダンスモデルを用いた従来法では、吸音の周波数依存性を反映するために複数の帯域に分けての計算が必要であった。これに対して、本実施形態に係る音場解析装置10は、一度の計算でより高精度に吸音の周波依存性が反映することができる。
In particular, in the conventional method using the frequency-independent surface impedance model, it was necessary to divide the calculation into multiple bands in order to reflect the frequency dependence of sound absorption. On the other hand, the sound
また、本実施形態に係る音場解析装置10は、ADE法に付随する常微分方程式である(9)、(10)、(11)式を4次精度の陽的ルンゲクッタ法で解くことで、高精度かつ安定的に吸音の周波数依存性を反映することができる。
In addition, the sound
(数値実験)
次に、上記実施形態に係る音場解析装置10によって音圧をどの程度精度良く推定できるかを数値実験により評価した。これにより、上記実施形態に係る音場解析装置10による精度改善の効果を検証した。
(Numerical experiment)
Next, numerical experiments were conducted to evaluate how accurately the sound pressure can be estimated by the sound
図7に、数値実験において解析対象の空間として使用した音響管モデルの断面を示す。数値実験では、図7に示すように、管端に、周波数に依存するインピーダンス境界をもつ幅0.01m×0.01m、長さ0.5mの音響管内の音場を解析対象とした。 FIG. 7 shows a cross section of the acoustic tube model used as the space to be analyzed in the numerical experiments. In the numerical experiment, as shown in FIG. 7, the sound field inside an acoustic tube with a width of 0.01 m×0.01 m and a length of 0.5 m having a frequency-dependent impedance boundary at the end of the tube was analyzed.
周波数に依存するインピーダンス境界には、Mikiの式でモデル化した、背後剛壁の多孔質吸音体の表面インピーダンスを与えた。なお、Mikiの式の詳細は、“Y. Miki, "Acoustical properties of porous materials - Modification of Delany-Bazley models -," J.Acoust.Soc.Jpn.(E),11,19-24(1990).”に記載されている。多孔質材の厚みと単位厚さ流れ抵抗は、それぞれ0.025m、R=13,000Pa s/m2とした。この多孔質吸音体の比音響アドミタンス比y(ω)は、nrp=ncp=4で有理関数近似した。 The frequency-dependent impedance boundary was given the surface impedance of the back rigid-wall porous sound absorber modeled by Miki's equation. For details of Miki's formula, see "Y. Miki, "Acoustical properties of porous materials - Modification of Delany-Bazley models -," J. Acoust. Soc. Jpn. (E), 11, 19-24 (1990). ."It is described in. The thickness and unit thickness flow resistance of the porous material were set to 0.025 m, respectively, and R=13,000 Pa s/m 2 . The specific acoustic admittance ratio y(ω) of this porous sound absorber was approximated by a rational function with n rp =n cp =4.
図8に、有理関数モデルで計算した多孔質吸音体の吸音率と理論値との比較を示す。図8では、縦軸に吸音率を表し、横軸に周波数を表している。図8に示すように、破線で示す有理関数モデルの計算値と、実線で示すその理論値とは、良好に一致している。 FIG. 8 shows a comparison between the sound absorption coefficient of the porous sound absorber calculated by the rational function model and the theoretical value. In FIG. 8, the vertical axis represents the sound absorption coefficient, and the horizontal axis represents the frequency. As shown in FIG. 8, the calculated values of the rational function model indicated by the dashed line and the theoretical values indicated by the solid line are in good agreement.
なお、音響管モデルは、要素長が0.005mの6節点立方体要素で離散化した。解析の初期条件として、音響管の入口への平面波入射を仮定した。解析上限周波数は10kHzとした。また、時間刻み幅Δtは、Δt=l・ΔtLIMITとした。ここで、ΔtLIMITは、立方体要素使用時の空気領域の計算における陽的な時間領域有限要素法の時間刻み幅Δtの安定限界値である(非特許文献1)。また、lは、数値的に求めた係数である。l=1の場合、安定性は空気領域と同等であり、l<1の場合、安定性が空気領域より悪化したことを表す。 The acoustic tube model was discretized with 6-node cubic elements with an element length of 0.005 m. As an initial condition of the analysis, we assumed a plane wave incident on the entrance of the acoustic tube. The analysis upper limit frequency was set to 10 kHz. Also, the time step width Δt was set to Δt=l·Δt LIMIT . Here, Δt LIMIT is the stability limit value of the time step size Δt of the explicit time domain finite element method in the calculation of the air domain when using cubic elements (Non-Patent Document 1). Also, l is a coefficient obtained numerically. For l=1, the stability was comparable to the air region, and for l<1, the stability was worse than the air region.
図9に、実施形態の手法で計算した多孔質吸音体の吸音率と理論値との比較を示す。図9では、図8と同様に、縦軸に吸音率を表し、横軸に周波数を表している。図9に示すように、破線で示す実施形態の手法の計算値と、実線で示すその理論値とは、良好に一致している。 FIG. 9 shows a comparison between the sound absorption coefficient of the porous sound absorber calculated by the method of the embodiment and the theoretical value. In FIG. 9, as in FIG. 8, the vertical axis represents the sound absorption coefficient, and the horizontal axis represents the frequency. As shown in FIG. 9, there is good agreement between the calculated values of the method of the embodiment, indicated by the dashed line, and its theoretical values, indicated by the solid line.
また、図9には、実施形態の手法の有効性を示すため、従来法による計算結果の例もあわせて一点鎖線で載せている。従来法では、境界の吸音を周波数に非依存な表面インピーダンスモデルを用いて、250Hzから8kHzの6つのオクターブバンド毎に分けて解析した場合を示している。従来法を用いる場合、6回の計算を要する上に、吸音の周波数特性が適切に反映されていないことが分かる。これに対して、実施形態の手法を使用することで、広帯域の吸音特性を1度の計算で高精度に反映できることが確認できる。また、実施形態の手法はl=0.85で安定的に計算可能であり、安定性に関しても優れている。 In addition, in FIG. 9, in order to show the effectiveness of the method of the embodiment, an example of the calculation result by the conventional method is also shown by a dashed line. The conventional method shows a case where boundary sound absorption is analyzed by dividing it into six octave bands from 250 Hz to 8 kHz using a frequency-independent surface impedance model. It can be seen that the conventional method requires six calculations and does not properly reflect the frequency characteristics of sound absorption. On the other hand, by using the method of the embodiment, it can be confirmed that the broadband sound absorption characteristics can be reflected with high accuracy in one calculation. Moreover, the method of the embodiment can be stably calculated with l=0.85, and is excellent in terms of stability.
(変形例)
以上に本発明の実施形態について説明したが、上記実施形態は一例であり、本発明の適用範囲はこれに限られない。すなわち、本発明の実施形態は種々の応用が可能であり、あらゆる実施の形態が本発明の範囲に含まれる。
(Modification)
Although the embodiment of the present invention has been described above, the above embodiment is an example, and the scope of application of the present invention is not limited to this. That is, the embodiments of the present invention can be applied in various ways, and all embodiments are included in the scope of the present invention.
例えば、上記実施形態では、(3)式に示した音圧ベクトルpの時間更新式は、音圧ベクトルpの時間一階微分を1次精度の前進差分近似で離散化ことにより得られた。また、(14)式に示した微分ベクトルvの時間更新式は、微分ベクトルvの時間一階微分を一次精度の後退差分近似で離散化し、更に音圧ベクトルpの時間一階微分を二次精度の中央差分近似で離散化することにより得られた。しかしながら、音圧ベクトルp及び微分ベクトルvの時間一階微分は、これ以外の手法で離散化されても良い。 For example, in the above-described embodiment, the time update formula for the sound pressure vector p shown in formula (3) was obtained by discretizing the time first-order differential of the sound pressure vector p by forward difference approximation with first-order precision. In addition, the time update formula for the differential vector v shown in equation (14) discretizes the time first derivative of the differential vector v by backward difference approximation with first-order precision, and furthermore, the time first derivative of the sound pressure vector p is quadratically The accuracy was obtained by discretizing with the median difference approximation. However, the time-first derivative of the sound pressure vector p and the differential vector v may be discretized by other methods.
上記実施形態では、行列設定部121は、非特許文献2に開示された修正積分則に習って、時間刻み幅Δt及び要素長hからクーラン数τを計算し、数値積分点αm,αkを設定した。しかしながら、行列設定部121は、他の手法で数値積分点αm,αkを設定しても良い。例えば、行列設定部121は、“Takumi Yoshida, Takeshi Okuzono, Kimihiro Sakagami, “Time domain room acoustic solver with fourth-order explicit FEM using modified time integration," Applied Sciences, 10, 3750(2020).”に開示された、解析対象の空間を離散化する要素の形状に依存しない数値積分点αm,αkを使用しても良い。この場合、数値積分点αm,αkは、クーラン数τを用いずに、“αm=±√(4/3)”及び“αk=±√(2/3)”と設定される。
In the above embodiment, the
上記実施形態では、音場解析装置10は、表示部14を備えており、出力部130は、数値計算部120により計算された複数の時間ステップでの音圧ベクトルpに基づく出力情報を表示部14に表示出力した。しかしながら、出力部130は、通信部15による通信を介して、出力情報を音場解析装置10の外部の装置に出力しても良い。そして、外部の装置が、例えば図4又は図5に示したように、出力部130から出力された出力情報を表示するようにしても良い。このように出力情報が外部の装置に出力される場合には、音場解析装置10は、表示部14を備えていなくても良い。
In the above embodiment, the sound
上記実施形態では、音場解析装置10の制御部11において、CPUがROM又は記憶部12に記憶されたプログラムを実行することによって、設定受付部110、数値計算部120及び出力部130のそれぞれとして機能した。しかしながら、制御部11は、専用のハードウェアであってもよい。専用のハードウェアとは、例えば単一回路、複合回路、プログラム化されたプロセッサ、ASIC(Application Specific Integrated Circuit)、FPGA(Field-Programmable Gate Array)、これらの組み合わせ等である。制御部11が専用のハードウェアである場合、各部の機能それぞれを個別のハードウェアで実現してもよいし、各部の機能をまとめて単一のハードウェアで実現してもよい。また、各部の機能のうち、一部を専用のハードウェアによって実現し、他の一部をソフトウェア又はファームウェアによって実現してもよい。このように、制御部11は、ハードウェア、ソフトウェア、ファームウェア、又は、これらの組み合わせによって、上述の各機能を実現することができる。
In the above-described embodiment, in the
本発明に係る音場解析装置の動作を規定する動作プログラムを既存のパーソナルコンピュータ、情報端末装置等に適用することで、当該パーソナルコンピュータ又は情報端末装置等を、本発明に係る音場解析装置として機能させることも可能である。 By applying an operation program that defines the operation of the sound field analysis device according to the present invention to an existing personal computer, information terminal device, or the like, the personal computer, information terminal device, or the like can be used as the sound field analysis device according to the present invention. It is also possible to make it work.
また、このようなプログラムの配布方法は任意であり、例えば、CD-ROM(Compact Disk ROM)、DVD(Digital Versatile Disk)、MO(Magneto Optical Disk)、メモリカード等のコンピュータ読み取り可能な記録媒体に格納して配布してもよいし、インターネット等の通信ネットワークを介して配布してもよい。 Any method of distributing such programs may be used. For example, computer-readable recording media such as CD-ROMs (Compact Disk ROMs), DVDs (Digital Versatile Disks), MOs (Magneto Optical Disks), and memory cards may be used. It may be stored and distributed, or may be distributed via a communication network such as the Internet.
本発明は、本発明の広義の精神と範囲を逸脱することなく、様々な実施形態及び変形が可能とされるものである。また、上述した実施形態は、この発明を説明するためのものであり、本発明の範囲を限定するものではない。すなわち、本発明の範囲は、実施形態ではなく、特許請求の範囲によって示される。そして特許請求の範囲内及びそれと同等の発明の意義の範囲内で施される様々な変形が、この発明の範囲内とみなされる。 The present invention is capable of various embodiments and modifications without departing from the broader spirit and scope of the invention. Moreover, the above-described embodiments are for explaining the present invention, and do not limit the scope of the present invention. That is, the scope of the present invention is indicated by the claims rather than the embodiments. Various modifications made within the scope of the claims and within the meaning of the invention equivalent thereto are considered to be within the scope of the present invention.
10 音場解析装置
11 制御部
12 記憶部
13 操作部
14 表示部
15 通信部
20 空間
110 設定受付部
120 数値計算部
121 行列設定部
123 音圧ベクトル計算部
124 補助関数計算部
125 微分ベクトル計算部
127 繰り返し部
130 出力部
10 Sound
Claims (7)
前記一階常微分方程式において、前記音圧ベクトルの時間一階微分を差分近似で離散化することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において、前記減衰項に含まれる周波数に依存する比音響アドミタンス比の逆フーリエ変換値と前記音圧ベクトルの時間一階微分との畳み込みをAuxiliary Differential Equations法で近似することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する数値計算部、を備える、
音場解析装置。 A first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in a space to be analyzed and a differential vector equivalent to the time first order differential of said sound pressure vector, wherein the sound wave at the boundary surrounding said space A sound field analysis device for analyzing a sound field using a time domain finite element method based on the first-order ordinary differential equation including an attenuation term for evaluating the attenuation of
In the first-order ordinary differential equation, the time update formula of the sound pressure vector obtained by discretizing the time first-order differential of the sound pressure vector by difference approximation, and in the first-order ordinary differential equation, the attenuation term A time update formula for the differential vector obtained by approximating the convolution of the included frequency-dependent inverse Fourier transform value of the specific acoustic admittance ratio with the time first derivative of the sound pressure vector using the Auxiliary Differential Equations method, according to , a numerical calculator that calculates the values of the sound pressure vector at a plurality of time steps;
Sound field analyzer.
前記数値計算部は、前記補助関数の時間発展を時間積分法により計算する、
請求項1に記載の音場解析装置。 The differential vector time update formula includes an auxiliary function obtained by rational function approximation of the specific acoustic admittance ratio,
The numerical calculation unit calculates the time evolution of the auxiliary function by a time integration method,
The sound field analysis device according to claim 1.
請求項2に記載の音場解析装置。 The numerical calculation unit calculates the time evolution of the auxiliary function using the Runge-Kutta method, the Euler method, the Crank-Nicholson method, or the Huyn method as the time integration method.
The sound field analysis device according to claim 2.
前記微分ベクトルの時間更新式は、第n時間ステップでの前記微分ベクトルの値を、第(n-1)時間ステップ及び第n時間ステップでの前記音圧ベクトルの値と、第(n-1)時間ステップでの前記微分ベクトルの値と、第n時間ステップでの前記補助関数の値と、により定める式である、
請求項2又は3に記載の音場解析装置。 The time update formula for the sound pressure vector is obtained by combining the value of the sound pressure vector at the nth time step with the value of the sound pressure vector at the (n−1)th time step and the (n−1)th time step. and the value of the differential vector at
The time update formula for the differential vector is such that the value of the differential vector at the nth time step is combined with the value of the sound pressure vector at the (n−1)th time step and the nth time step, and the (n−1)th ) is defined by the value of the differential vector at the time step and the value of the auxiliary function at the nth time step,
The sound field analysis device according to claim 2 or 3.
請求項1から4のいずれか1項に記載の音場解析装置。 The time update formula for the differential vector discretizes the time first derivative of the differential vector by backward difference approximation with first-order accuracy, and further discretizes the time first derivative of the sound pressure vector by central difference approximation with second-order accuracy. obtained by
The sound field analysis device according to any one of claims 1 to 4.
前記一階常微分方程式において、前記音圧ベクトルの時間一階微分を差分近似で離散化することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において、前記減衰項に含まれる周波数に依存する比音響アドミタンス比の逆フーリエ変換値と前記音圧ベクトルの時間一階微分との畳み込みをAuxiliary Differential Equations法で近似することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する、
音場解析方法。 A first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in a space to be analyzed and a differential vector equivalent to the time first order differential of said sound pressure vector, wherein the sound wave at the boundary surrounding said space A sound field analysis method for analyzing a sound field using a time domain finite element method based on the first-order ordinary differential equation including an attenuation term for evaluating the attenuation of
In the first-order ordinary differential equation, the time update formula of the sound pressure vector obtained by discretizing the time first-order differential of the sound pressure vector by difference approximation, and in the first-order ordinary differential equation, the attenuation term A time update formula for the differential vector obtained by approximating the convolution of the included frequency-dependent inverse Fourier transform value of the specific acoustic admittance ratio with the time first derivative of the sound pressure vector using the Auxiliary Differential Equations method, according to , calculating values of said sound pressure vector at multiple time steps;
Sound field analysis method.
解析対象の空間における音圧の分布を示す音圧ベクトルと、前記音圧ベクトルの時間一階微分と等価な微分ベクトルと、に関する一階常微分方程式であって、前記空間を囲む境界での音波の減衰を評価するための減衰項を含む前記一階常微分方程式に基づく時間領域有限要素法を用いて音場を解析する音場解析装置として機能させるプログラムであって、
前記コンピュータを、
前記一階常微分方程式において、前記音圧ベクトルの時間一階微分を差分近似で離散化することにより得られる前記音圧ベクトルの時間更新式と、前記一階常微分方程式において、前記減衰項に含まれる周波数に依存する比音響アドミタンス比の逆フーリエ変換値と前記音圧ベクトルの時間一階微分との畳み込みをAuxiliary Differential Equations法で近似することにより得られる前記微分ベクトルの時間更新式と、に従って、複数の時間ステップでの前記音圧ベクトルの値を計算する数値計算部、として機能させる、
プログラム。 the computer,
A first order ordinary differential equation relating to a sound pressure vector indicating the distribution of sound pressure in a space to be analyzed and a differential vector equivalent to the time first order differential of said sound pressure vector, wherein the sound wave at the boundary surrounding said space A program that functions as a sound field analysis device that analyzes a sound field using a time domain finite element method based on the first-order ordinary differential equation including an attenuation term for evaluating the attenuation of
said computer,
In the first-order ordinary differential equation, the time update formula of the sound pressure vector obtained by discretizing the time first-order differential of the sound pressure vector by difference approximation, and in the first-order ordinary differential equation, the attenuation term A time update formula for the differential vector obtained by approximating the convolution of the included frequency-dependent inverse Fourier transform value of the specific acoustic admittance ratio with the time first derivative of the sound pressure vector using the Auxiliary Differential Equations method, according to , a numerical calculator that calculates the values of the sound pressure vector at multiple time steps;
program.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2021008481A JP7496972B2 (en) | 2021-01-22 | 2021-01-22 | Sound field analysis device, sound field analysis method and program |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2021008481A JP7496972B2 (en) | 2021-01-22 | 2021-01-22 | Sound field analysis device, sound field analysis method and program |
Publications (2)
Publication Number | Publication Date |
---|---|
JP2022112624A true JP2022112624A (en) | 2022-08-03 |
JP7496972B2 JP7496972B2 (en) | 2024-06-10 |
Family
ID=82657268
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2021008481A Active JP7496972B2 (en) | 2021-01-22 | 2021-01-22 | Sound field analysis device, sound field analysis method and program |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP7496972B2 (en) |
Family Cites Families (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP4059478B2 (en) | 2002-02-28 | 2008-03-12 | パイオニア株式会社 | Sound field control method and sound field control system |
JP2007040734A (en) | 2005-08-01 | 2007-02-15 | Nec Corp | Transmission mode recommend system, sonar system, transmission mode recommend method, and its program and storage medium |
JP2007188164A (en) | 2006-01-11 | 2007-07-26 | Vanderplaats Design Optimization Consulting Inc | Acoustic structure compound optimal design analysis method and its optimal design system and its analysis program and recording medium with its analysis program recorded |
US9135699B2 (en) | 2012-03-15 | 2015-09-15 | Siemens Aktiengesellschaft | Method and system for hemodynamic assessment of aortic coarctation from medical image data |
CN112001133B (en) | 2020-08-21 | 2023-05-12 | 中国船舶科学研究中心 | Fluid-solid acoustic coupling calculation method based on ship three-dimensional acoustic elastic time domain analysis |
-
2021
- 2021-01-22 JP JP2021008481A patent/JP7496972B2/en active Active
Also Published As
Publication number | Publication date |
---|---|
JP7496972B2 (en) | 2024-06-10 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US6985836B2 (en) | Computer-aided engineering method and apparatus for predicting a quantitative value of a physical property at a point from waves generated by or scattered from a body | |
Cicirello et al. | The vibro-acoustic analysis of built-up systems using a hybrid method with parametric and non-parametric uncertainties | |
Squicciarini et al. | The effect of different combinations of boundary conditions on the average radiation efficiency of rectangular plates | |
Cuenca et al. | The image source method for calculating the vibrations of simply supported convex polygonal plates | |
Karimi et al. | Analytical and numerical prediction of acoustic radiation from a panel under turbulent boundary layer excitation | |
Machado et al. | Spectral element-based method for a one-dimensional damaged structure with distributed random properties | |
Hartmann et al. | High-frequency structure-and air-borne sound transmission for a tractor model using dynamical energy analysis | |
US10319121B2 (en) | Acoustic performance calculation device, acoustic performance calculation method, and acoustic performance calculation program | |
António et al. | A three-dimensional acoustics model using the method of fundamental solutions | |
Ou et al. | The effects of elastic supports on the transient vibroacoustic response of a window caused by sonic booms | |
Miller III et al. | Time-domain impedance boundary conditions for acoustic reduced order finite element simulations | |
JP5024613B2 (en) | Sound field analyzer | |
JP7496973B2 (en) | Sound field analysis device, sound field analysis method and program | |
JP6886890B2 (en) | Decay time analysis methods, instruments, and programs | |
JP7496972B2 (en) | Sound field analysis device, sound field analysis method and program | |
JP7285514B2 (en) | SOUND FIELD ANALYZER, SOUND FIELD ANALYSIS METHOD AND PROGRAM | |
Du et al. | Fast simulation of non-linear pulsed ultrasound fields using an angular spectrum approach | |
Prislan et al. | Ray-trace modeling of acoustic Green's function based on the semiclassical (eikonal) approximation | |
JP7285513B2 (en) | SOUND FIELD ANALYZER, SOUND FIELD ANALYSIS METHOD AND PROGRAM | |
Kelly et al. | Approximate analytical time-domain Green's functions for the Caputo fractional wave equation | |
Yu et al. | Spectral decomposition-based explicit integration method for fully non-stationary seismic responses of large-scale structures | |
JP7113168B2 (en) | Sound pressure calculation method and sound pressure calculation program | |
Casadei et al. | Towards sound absorption in a cylindrical lined duct using CFD with time-domain impedance boundary condition | |
Popescu et al. | A finite volume‐based high‐order, Cartesian cut‐cell method for wave propagation | |
Kurowski et al. | Reinforcement learning algorithm and FDTD-based simulation applied to schroeder diffuser design optimization |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20210409 |
|
A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20231101 |
|
A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20240426 |
|
TRDD | Decision of grant or rejection written | ||
A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20240507 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20240520 |
|
R150 | Certificate of patent or registration of utility model |
Ref document number: 7496972 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |