JP6086235B2

JP6086235B2 - 物理量シミュレーションシステム

Info

Publication number: JP6086235B2
Application number: JP2013127371A
Authority: JP
Inventors: 崇石塚; 寛大久保
Original assignee: Shimizu Corp
Current assignee: Shimizu Corp
Priority date: 2013-06-18
Filing date: 2013-06-18
Publication date: 2017-03-01
Anticipated expiration: 2033-06-18
Also published as: JP2015001923A

Description

本発明は、コンピューター上で音波や電磁波の時間に伴う伝搬をシミュレーションするために好適な物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムに関する。

従来より、コンピューター上で音波や電磁波の時間に伴う伝搬を波動理論に基づきシミュレーションする手法として、ＦＤＴＤ（ＦｉｎｉｔｅＤｉｆｆｅｒｎｃｅＴｉｍｅ
Ｄｏｍａｉｎ）法やＣＩＰ（ＣｏｎｓｔｒａｉｎｅｄＩｎｔｅｒｐｏｌａｔｉｏｎＰｒｏｆｉｌｅ）法が提案されている。

屋外等の開放空間の波動伝搬のシミュレーションを行う場合、ＦＤＴＤ法やＣＩＰ法のような空間離散化型解法では、解析対象領域を有限で打ち切る必要がある。その際、解析対象領域外周の仮想境界面からの反射波が、シミュレーション結果に影響を与えないようにする必要がある。

ＣＩＰ法では、シミュレーションの計算過程において波動を進行波成分と後退波成分（ある軸方向の正方向と負方向に進む波動成分）に分けるため、外周境界面で解析対象領域内に向かう成分を０と置き換えることで反射波を抑えることができる。

ただし、この方法は、境界面に対して正面から波動が入射する場合、すなわち、境界面と波動の進行方向が直交する場合には有効であるが、境界面に対して斜めに入射する場合には反射波を生じる。

一方で、ＰＭＬ（ＰｅｒｆｅｃｔｌｙＭａｔｃｈｅｄＬａｙｅｒ）吸収境界は、あらゆる入射方向に対して有効な手法として、ＦＤＴＤ法では広く用いられている（非特許文献１参照）。また、ＰＭＬ吸収境界のＣＩＰ法への適用手法は、非特許文献２おいて電磁波解析への適用例が開示されている。
「音環境の数値シミュレーション―波動音響解析の技法と応用―」、日本建築学会編、日本建築学会、２０１１年 "ＩｍｐｌｅｍｅｎｔａｔｉｏｎｏｆｔｈｅＰｅｒｆｅｃｔＭａｔｃｈｅｄＬａｙｅｒｔｏｔｈｅＣＩＰＭｅｔｈｏｄ"，ＹｏｓｈｉａｋｉＡｎｄｏａｎｄＭａｓａｓｈｉＨａｙａｋａｗａ，ＩＥＩＣＥＴｒａｎｓ．Ｅｌｅｃｔｒｏｎ．，Ｖｏｌ．Ｅ８９−Ｃ，Ｎｏ．５，２００６．

しかしながら、非特許文献２記載の手法は、Ｍ型ＣＩＰ法に対してＰＭＬ吸収境界を実装する手法であり、より計算精度の高いＣ型ＣＩＰ法に対する実装方法はこれまで提供されていない、という問題があった。

また、非特許文献２記載の手法は、ＰＭＬ吸収境界層内部での波動の減衰を表す微分方程式について、微分係数を差分近似して計算を行っている。この微分の差分近似により計算誤差が生じてしまうこととなる。特に、非特許文献２記載の手法は、時間微分を１次精度の前進差分で近似しており誤差を生じやすい、という問題があった。

この発明は、上記課題を解決するものであって、請求項１に係る発明は、解析対象領域と前記解析対象領域の周縁に配される吸収境界層とに設定されている格子点に対応する物理量をシミュレーションする物理量シミュレーションシステムであって、前記吸収境界層における吸収係数を格子点の位置の関数として定める吸収境界設定ステップと、移流方程式に基づいて、前記吸収境界層における格子点の物理量を、ＣＩＰ法により数値的に演算する第１演算ステップと、前記第１演算ステップにより演算された物理量と、前記関数と、減衰方程式とに基づいて、前記吸収境界層における格子点の物理量を、解析的に演算する第２演算ステップとをコンピューターで実行することを特徴とする。

また、請求項２に係る発明は、請求項１に記載の物理量シミュレーションシステムにおいて、前記第２演算ステップには、２階微分された物理量が含まれることを特徴とする。

また、請求項３に係る発明は、請求項１又は請求項２に記載の物理量シミュレーションシステムにおいて、物理量は、特定方向に移流する波動の成分を表す特性値である。

また、請求項４に係る発明は、請求項３に記載の物理量シミュレーションシステムにおいて、前記特性値から音圧を求めることを特徴とする。

また、請求項５に係る発明は、請求項３又は請求項４に記載の物理量シミュレーションシステムにおいて、前記特性値から粒子速度を求めることを特徴とする。

本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムでは、第１演算ステップにより演算された物理量と、関数と、減衰方程式とに基づいて、吸収境界層における格子点の物理量を演算する第２演算ステップにおいて、これを解析的に行うので、微分係数の差分近似は含まれていないため、近似誤差を生じることなく吸収境界層内部の波動伝搬シミュレーションを高精度で行うことが可能となる。

また、本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムは、Ｍ型ＣＩＰ法に比べ、より計算精度の高いＣ型ＣＩＰ法に対するＰＭＬ吸収境界の実装方法を提供するものである。

また、本発明により、特定方向に移流する波動の成分を表す特性値であるｆ_x±、ｇ_x±及び移流方向と直交する方向への微分値η_x±、ξ_x±などを用いるＣ型ＣＩＰ法による波動伝搬シミュレーションにおいて、仮想境界面からの反射波を低減することができる。

本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムは、ＰＭＬ吸収境界内の波動伝搬のみでなく、一般的な減衰のある場におけるＣＩＰ法を用いた波動伝搬シミュレーションに対しても適用でき、減衰しつつ伝搬する波動の高精度な予測計算が可能になる。

また、本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムは、音波の伝搬解析以外にも、電磁波など他の波動一般の伝搬解析にも適用可能である。

本発明の実施形態に係る物理量シミュレーションシステムを構成するコンピューターの一例を示す図である。解析対象領域の周囲に設定される吸収境界層を説明する図である。本発明に係る物理量シミュレーション方法をコンピューターに実行させるためのフローチャート例を示す図である。

以下、本発明の実施の形態を図面を参照しつつ説明する。図１は本発明の実施形態に係る物理量シミュレーションシステムを構成するコンピューターの一例を示す図である。図１において、１０はシステムバス、１１はＣＰＵ（ＣｅｎｔｒａｌＰｒｏｃｅｓｓｉｎｇＵｎｉｔ）、１２はＲＡＭ（ＲａｎｄｏｍＡｃｃｅｓｓＭｅｍｏｒｙ）、１３はＲＯＭ（ＲｅａｄＯｎｌｙＭｅｍｏｒｙ）、１４は外部情報機器との通信を司る通信制御部、１５はキーボードコントローラなどの入力制御部、１６は出力制御部、１７は外部記憶装置制御部、１８はキーボード、ポインティングデバイス、マウスなどの入力機器からなる入力部、１９は印刷装置などの出力部、２０はＨＤＤ（ＨａｒｄＤｉｓｋＤｒｉｖｅ）等の外部記憶装置、２１はグラフィック制御部、２２はディスプレイ装置をそれぞれ示している。

図１において、ＣＰＵ１１は、ＲＯＭ１３内のプログラム用ＲＯＭ、或いは、大容量の外部記憶装置２０に記憶されたプログラム等に応じて、外部機器と通信することでデータを検索・取得したり、また、図形、イメージ、文字、表等が混在した出力データの処理を実行したり、更に、外部記憶装置２０に格納されているデータベースの管理を実行したり、などといった演算処理を行うものである。

また、ＣＰＵ１１は、システムバス１０に接続される各デバイスを統括的に制御する。ＲＯＭ１３内のプログラム用ＲＯＭあるいは外部記憶装置２０には、ＣＰＵ１１の制御用の基本プログラムであるオペレーティングシステムプログラム（以下ＯＳ）等が記憶されている。また、ＲＯＭ１３あるいは外部記憶装置２０には出力データ処理等を行う際に使用される各種データが記憶されている。メインメモリーであるＲＡＭ１２は、ＣＰＵ１１の主メモリ、ワークエリア等として機能する。

入力制御部１５は、キーボードや不図示のポインティングデバイスからの入力部１８を制御する。また、出力制御部１６は、プリンタなどの出力部１９の出力制御を行う。

外部記憶装置制御部１７は、ブートプログラム、各種のアプリケーション、フォントデータ、ユーザーファイル、編集ファイル、プリンタドライバ等を記憶するＨＤＤ（ＨａｒｄＤｉｓｋＤｒｉｖｅ）や、或いはフロッピーディスク（ＦＤ）等の外部記憶装置２０へのアクセスを制御する。本発明の物理量シミュレーション方法を実現するシステムプログラムは、上記のような外部記憶装置２０に記憶されている。また、グラフィック制御部２１は、ディスプレイ装置２２に表示する情報を描画処理するための構成である。

また、通信制御部１４は、ネットワークを介して、外部機器と通信を制御するものであり、これによりシステムが必要とするデータを、インターネットやイントラネット上の外部機器が保有するデータベースから取得したり、外部機器に情報を送信したりすることができるように構成される。

外部記憶装置２０には、ＣＰＵ１１の制御プログラムであるオペレーティングシステムプログラム（以下ＯＳ）以外に、本発明の物理量シミュレーションシステムをＣＰＵ１１上で動作させるシステムプログラム、及びこのシステムプログラムで用いるデータなどがインストールされ保存・記憶されている。

本発明の物理量シミュレーション方法を実現するシステムプログラムで利用されるデータとしては、基本的には外部記憶装置２０に保存されていることが想定されているが、場合によっては、これらのデータを、通信制御部１４を介してインターネットやイントラネット上の外部機器から取得するように構成することも可能である。また、本発明の物理量シミュレーション方法を実現するシステムプログラムで利用されるデータを、ＵＳＢメモリやＣＤ、ＤＶＤなどの各種メディアから取得するように構成することもできる。

次に、上記のようなシステム構成のコンピューターにより実行可能な本発明に係る物理量シミュレーション方法について、以下説明する。

以下の例では、上記のようなコンピューター上で、解析対象領域における音波の時間に伴う伝搬をシミュレーションする例を説明するが、本発明は電磁波など他の波動一般の伝搬解析にも適用可能である。

図２は解析対象領域の周囲に設定される吸収境界層を説明する図である。解析対象領域及び前記解析対象領域の周縁に配される吸収境界層には、不図示のｘ方向と平行なグリッド及びｙ方向に平行なグリッドが定義されており、本発明に係る物理量シミュレーション方法においては、このグリッドが交差する点として定義される各格子点における物理量（２次元の音波の伝搬の場合は、音圧ｐと粒子速度ｕ_x、ｕ_y）がシミュレーションされる。また、音圧ｐと粒子速度ｕ_x、ｕ_yを算出する上で、ＣＩＰ法では、特定方向に移流する波動の成分を表す特性値ｆやｇが演算される。

解析対象領域におけるＣＩＰ法の適用については、
・「Ｃ型ＣＩＰ法を用いた音場解析に関する検討」、信学技報、Ｖｏｌ．１０６，Ｎｏ．４８１，ｐｐ．１７−２２，ｓ２００６−９８，２００６．
・「ＣＩＰ法による音波伝搬シミュレーションにおけるサブグリッド・テクニックの精度比較」，日本音響学会２０１１年秋季研究発表会講演論文集，ｐｐ．１４９７−１５００，２０１１．
などに記載のＣＩＰ法に関する記載を参照して援用する。

ここでは、吸収境界層における物理量シミュレーションについて、より詳しく説明する。本発明に係る物理量シミュレーション方法は、ＣＩＰ法による波動伝搬シミュレーションにおいてＰＭＬ吸収境界（ＰＭＬ層）を実現する手段として、より精度の高いＣ型ＣＩＰ法への実装方法、および差分近似を用いない高精度な実装方法を提供するものである。

以下では、２次元の音響伝搬シミュレーションを例に記述するが、本発明は振動、電磁波など音響以外の波動伝搬シミュレーション、またこれらの３次元シミュレーションに対しても有効である。

ＰＭＬ吸収境界とは、解析対象領域の外周に設定する、波動が徐々に減衰するように物性値を与えた吸収層である。このＰＭＬ吸収境界は、図２の解析対象領域の外周における（Ａ）、（Ｂ）、（Ｃ）に示す領域である。音響伝搬シミュレーションでは、ＰＭＬ吸収境界として、吸収係数Ｒを徐々に大きくした吸収層を解析対象領域の外周部に設定する。

このとき、吸収係数は音波を吸収したい軸方向にのみ設定する。すなわち、図２において右側の解析対象領域端では、ｘ方向に進行する音波のみ吸収したいので、Ｒ_xを徐々に
大きくなるように設定する一方で、Ｒ_y＝０とする。すなわち、図２の（Ａ）の領域では
、下式（１）、（２）のように設定される。

また、図２において、上側の解析対象領域端ではｙ方向に進行する音波のみ吸収したいので、Ｒ_yを徐々に大きくなるように設定する一方で、Ｒ_x＝０とする。すなわち、図２の（Ｂ）の領域では、下式（３）、（４）のように設定される。

図２右上の隅部においては、ｘ，ｙ両方向に進行する音波を吸収したいので、Ｒ_x及び
Ｒ_yを徐々に大きくなるように設定する。すなわち、図２の（Ｃ）の領域では、下式（５
）、（６）のように設定される。

Ｒ_x及びＲ_yはＰＭＬ層の表面からの距離の関数として式（１）、（４）〜（６）のように与える。

ここで、Ｒ_maxは吸収係数の最大値（定数）、ｘ_o及びｙ_oはＰＭＬ層表面のｘ及びｙ座
標、ｄはＰＭＬ層の厚さ、ｍはＰＭＬ層内での吸収係数の分布を求める次数（ｍ＝１，２，３，・・・）である。

吸収係数Ｒによる減衰のある音場のＣＩＰ法によるシミュレーション方法を以下に示す。ＣＩＰ法では、２次元以上の音場の解析シミュレーションは、支配方程式を各軸方向に分離し、それぞれの軸方向への波動伝搬を互い違いに計算することで行う。以下では、ｘ軸方向への音波の伝搬計算の方法を記述する。ｙ方向への音波の伝搬は、以下の各式においてｘとｙを入れ替えて同じ手順に従えば計算できる。

吸収係数による減衰のある場におけるｘ軸方向の音波伝搬の支配方程式は、以下の式（
７）、（８）、（９）の通りである。

ここで、ｐは音圧、ｕ_xはｘ方向粒子速度、ｃは音速、ρは媒質密度、Ｚは媒質の特性
インピーダンスである。ｆ_x±は、±ｘ方向へ移流する波動の成分を表す特性値である。

ＣＩＰ法では、移流を表す移流方程式（下式の式（１０）に示すもの）と、減衰を表す減衰方程式（下式の式（１1）に示すもの）に分離して計算を行う。本発明においては、
移流方程式を演算する上では数値計算を行い、減衰方程式を演算する上ではこれを解析的に演算する。また、時間ステップｎにおける演算では、まず移流方程式の演算を実行することで、中間ステップの値を得て、次に減衰方程式の演算を実行することで、時間ステップ（ｎ＋１）の値を得る。

まず、移流を表す移流方程式の演算について説明する。

移流方程式をｘに関して偏微分することで次式を得る。

ここで、

である。

式（１０）と式（１２）とは、ｆ_x±とその微分値ｇ_x±に関する移流方程式であるので、ＣＩＰ法による移流計算が可能である。また、ある時間ステップの中間ステップにおける移流計算の結果を、ｆ^n* _x±およびｇ^n* _x±とする。

Ｃ型ＣＩＰ法による２次元音場の解析では、ｆ_x±及びｇ_x±のｙ方向微分値の移流計算を行う必要がある。式（１０）と式（１２）をｙ方向で偏微分すると次式を得る。

ここで、

である。

また、

ここで、

である。

式（１４）と式（１６）とは、η_x±とその微分値ξ_x±に関する移流方程式であるので、ＣＩＰ法による移流計算が可能である。ある時間ステップ（中間ステップ）における移流計算の結果を、η^n* _x±およびξ^n* _x±とする。

次に、減衰方程式の演算について説明する。

式（１１）の減衰に関する偏微分方程式は解析的に解くことが可能である。時間ステップ（ｎ＋１）におけるｆ_x±は次式（１８）で計算できる。

ここで、Δｔは時間ステップ幅である。

また、ｇ_x±に関する減衰は、式（１８）をｘ方向に偏微分した次の式（１９）で解析
的に計算を行うことができる。

また、式（１９）中の∂ｒ／∂ｘは、式（１）、（９）より以下の式で計算することができる。

また、η_x±に関する減衰については、式（１１）をｙ方向に偏微分すると、下式（２
１）が得られる。

これを、ｆ_x±と同様に解析的に解いて、

により計算することができる。
また、ξ_x±に関する減衰は、式（２２）をｘ方向に偏微分した次の式で（２３）計算で
きる。

なお、上記のような演算で得られるｆ_x±から、音圧ｐ及び粒子速度ｕ_xは、それぞれ下式（２４）及び（２５）によって計算することができる。

なお、図２に示されていないが図２において左側の解析領域端にＰＭＬ吸収境界を設定した場合、Ｒ_xは左側の解析領域端に向けて徐々に大きくなるように与えるので、以上の
計算において、式（１）及び（２０）をそれぞれ次のように置き換えて計算する。

ここで、ｘ_O’は左側解析領域端に設定したＰＭＬ層表面のｘ座標である。

次に、以上のような本発明に係る物理量シミュレーション方法を、コンピューターで演算する際のアルゴリズムの一例を説明する。図３は本発明に係る物理量シミュレーション方法をコンピューターに実行させるためのフローチャート例を示す図である。

図３において、ステップＳ１００でシミュレーションの処理が開始されると、続くステップＳ１０１では、ステップを表す変数ｎに値として１がセットされる。

続く、ステップＳ１０２では、前記移流方程式を、ＣＩＰ法を用いて数値的に演算する。特許請求の範囲では、この演算を第１演算と称している。このステップＳ１０２の第１演算が実行されると、（２８）及び（２９）に示すような中間ステップにおける値を得ることができる。

ステップＳ１０３においては、吸収境界層内であるか否かが判定される。ステップＳ１０３における判定がＹＥＳである場合には、ステップＳ１０５で、これまで説明したように、減衰方程式を解析的に演算する。特許請求の範囲では、この演算を第２演算と称している。このステップＳ１０５の第２演算が実行されると、（３０）及び（３１）に示すような中間ステップからステップ（ｎ＋１）における値を得ることができる。

一方、ステップＳ１０３における判定がＮＯである場合には、一般のＣＩＰ法通り、下式（３２）乃至（３５）によりステップ（ｎ＋１）における値を得る。

ステップＳ１０６では、全ての格子点で、ｆⁿ⁺¹、ｇⁿ⁺¹、ηⁿ⁺¹、ξⁿ⁺¹の各値が演算されたか判定される。この判定がＮＯであれば、ステップＳ１０８に進み、次の格子点について演算するためにループする。

一方、判定がＹＥＳであれば、ステップＳ１０７に進み、最終ステップｎ_lastに到達したか否かが判定される。この判定がＮＯであれば、ステップＳ１０９に進み、ｎの値が１インクリメントされ、再び演算ステップをループする。

ステップＳ１０７における判定がＹＥＳとなれば、ステップＳ１１０に進み、処理を終了する。

以上、本発明に係る物理量シミュレーション方法及びそれを用いた物理量シミュレーションシステムでは、第１演算ステップにより演算された物理量と、関数と、減衰方程式とに基づいて、吸収境界における格子点の物理量を演算する第２演算ステップにおいて、これを解析的に行うので、微分係数の差分近似は含まれていないため、近似誤差を生じることなく吸収境界層内部の波動伝搬シミュレーションを高精度で行うことが可能となる。

１０・・・システムバス
１１・・・ＣＰＵ（ＣｅｎｔｒａｌＰｒｏｃｅｓｓｉｎｇＵｎｉｔ）
１２・・・ＲＡＭ（ＲａｎｄｏｍＡｃｃｅｓｓＭｅｍｏｒｙ）
１３・・・ＲＯＭ（ＲｅａｄＯｎｌｙＭｅｍｏｒｙ）
１４・・・通信制御部
１５・・・入力制御部
１６・・・出力制御部
１７・・・外部記憶装置制御部
１８・・・入力部
１９・・・出力部
２０・・・外部記憶装置
２１・・・インターフェイス部
２１・・・グラフィック制御部
２２・・・ディスプレイ装置

Claims

解析対象領域と前記解析対象領域の周縁に配される吸収境界層とに設定されている格子点に対応する物理量をシミュレーションする物理量シミュレーションシステムであって、
前記吸収境界層における吸収係数を格子点の位置の関数として定める吸収境界設定ステップと、
移流方程式に基づいて、前記吸収境界層における格子点の物理量を、ＣＩＰ法により数値的に演算する第１演算ステップと、
前記第１演算ステップにより演算された物理量と、前記関数と、減衰方程式とに基づいて、前記吸収境界層における格子点の物理量を、解析的に演算する第２演算ステップとをコンピューターで実行することを特徴とする物理量シミュレーションシステム。
前記第２演算ステップには、２階微分された物理量が含まれることを特徴とする請求項１に記載の物理量シミュレーションシステム。
物理量は、特定方向に移流する波動の成分を表す特性値である請求項１又は請求項２に記載の物理量シミュレーションシステム。
前記特性値から音圧を求めることを特徴とする請求項３に記載の物理量シミュレーションシステム。
前記特性値から粒子速度を求めることを特徴とする請求項３又は請求項４に記載の物理量シミュレーションシステム。